gpc_1a.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function [sys, x0,str,ts]=gpc_1a(t,x,u,flag,N1,N2,Nu,r,rho,k_delay,...
             B_pocz,A_pocz,P_pocz,alfa,ts)
nA=length(A_pocz); nB=length(B_pocz)-1; k=nA+nB+1;
kpend=k*k; kt=kpend+1; ktend=kpend+k; kf=ktend+1;
kfend=ktend+k; ky=kfend+1; ku=ky+1; kend=ky+k_delay;
P=zeros(k,k); x=x(:)';
switch flag
   case 0
      sizes = simsizes; % 读入系统变量的默认值
      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:kpend)=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:kpend); 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_max=max([nA+1,nB+k_delay]);
      num=[zeros(1,k_delay-1),Theta(nA+1:k),zeros(1,k_max-nB-k_delay)];
      den=[1, Theta(1:nA),zeros(1,k_max-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_delay-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_delay)];
      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