adaptivegpc.m

自适应广义预测控制S函数程序

function [sys,x0,str,ts] = gpc_la(t,x,u,flag,N1,N2,Nu,r,rho,k_d,P_pocz,alfa,ts)
num=1;
den=[2 8 1];
[num1,den1]=c2dm(num,den,ts);
A_pocz=num1;%zeros(1,size(den)); A_pocz(1,1)=1;
B_pocz=den1;%zeros(1,size(num));

nA = length(A_pocz);nB = length(B_pocz)-1; k =nA+nB+1;
kp = k*k; kt = kp+1; ktend = kp+k; kf = ktend+1;
kfend = ktend+k; ky = kfend+1; ku = ky+1; kend = ky+k_d;
P = zeros(k,k); x= x(:)';
switch flag
    case 0
        sizes = simsizes;%Read the default values of the system varibles
        sizes.NumContStates = 0; sizes.NumDiscStates = kend;
        sizes.NumOutputs = 1; sizes.NumInputs =2;
        sizes.DirFeedthrough = 0; sizes.NumSampleTimes = 1;
        sys = simsizes(sizes); str = []; ts = [ts 0];
        x0 = zeros(1,kend); x0(1:k+1:kp) = P_pocz*ones(k,1);
        x0(kt:kt+nA-1) = A_pocz; x0(kt+nA:ktend) = B_pocz;
    case 2
        Phi = [x(ky),x(kf:kf+nA-2),x(kend),x(kf+nA:kfend-1)];
        P(:) = x(1:kp); P = (1/alfa)*(P-(P*Phi'*Phi*P)/(alfa+Phi*P*Phi'));
        Theta = x(kt:ktend)+Phi*P*(u(2)-Phi*x(kt:ktend)');
        k_M = max([nA+1,nB+k_d]);
        num = [zeros(1,k_d-1),Theta(nA+1:k),zeros(1,k_M-nB-k_d)];
        den = [1,Theta(1:nA),zeros(1,k_M-nA-1)]; h = dstep(num,den,N2);
        for i = 1:Nu, Qt(1:N2,i) = [zeros(i-1,1); h(1:N2-i+1)]; end
        Q = Qt(N1:N2,:); q = [1,zeros(1,Nu-1)]*inv(Q'*Q+r*eye(Nu))*Q';
        [w,xw] = dlsim(rho,1-rho,1,0,u(1)*ones(N2+1,1),u(2));
        A = [1,Theta(1:nA)]; B = Theta(nA+1:k); Bm = [B,0]; Bm = Bm-[0,B];
        Am = [A,0]; Am = Am-[0,A]; Ared = Am(2:nA+2);
        Bred = [zeros(1,k_d-1),Bm]; Y = [u(2),-x(ky),-x(kf:kf+nA-2)];
        U = [x(ku),x(ku:kend),x(kf+nA:kf+nA+nB-1)];
        for i = 1:N2
            yp(i) = -Ared*Y'+Bred*U'; Y = [yp(i),Y(1:nA)];
            U = [U(1),U(1:nB+k_d)];
        end
        nu = x(ku)+q*(w(N1+1:N2+1)-yp(N1:N2)');
        sys = [P(:)',Theta,x(ky),x(kf:kf+nA-2),x(kend),x(kf+nA:kfend-1),-u(2),nu,x(ku:kend-1)];
    case 3, sys = x(ku);
    otherwise, sys = [];
end