shushanglizi.m

此为广义预测控制的一个程序

N1=1;N2=10;Nu=3;P=3;          % N1为最小预测长度，N2为最大预测时域，控制时域Nu，P为预测时域吧
A=[1 -0.8 0 0];
A1=[1 -1.8 0.8 0];     % A1=A*△                         
B=[0.4 0.6 0 0];
na=2;nb=3;                    % na为A的次数，好像也是F的z次数，
S=zeros(P,na+1);              % 
S(1,:)=[1.8 -0.8 0];  % 怎么来的？？
for j=1:P-1
    for i=1:na
S(j+1,i)=S(j,i+1)-A1(i+1)*S(j,1)
S(j+1,na+1)=-S(j,1)*A1(na+2)     
    end
end                           % 递推求解书上的F（即S）
R=zeros(P,P);
r0=1;r=[];r(1)=r0;            % 定义E（即R）
for j=1:P-1
R(j+1,j+1)=S(j,1);            
r(j+1)=R(j+1,j+1);            
end
for j=1:2
    for i=1:j
R(j+1,i)=r(i);
R(j,i)=R(j+1,i);
    end
end                           % 递推求解E （即R）
G=zeros(P,P+nb);              % 书上定义Ｇ为ＮｘＭ。其中Ｎ为建模长度，Ｍ为控制长度？
for j=1:nb+1                  % G=E*B
G(1,j)=B(j);
end
for j=1:P-1
    for i=1:j+nb+1
        if i<=j
            bij=0;
        else 
            bij=B(i-j);
        end
        G(j+1,i)=G(j,i)+S(j,1)*bij;
    end
end
G1=zeros(P,Nu);G2=zeros(P,nb);
for i=1:P
    for j=1:Nu
        if i==j
            G1(i,j)=G(1,1);
        else if i==j+1
                G1(i,j)=G(2,2);
            else if i==j+2
                    G1(i,j)=G(3,3);
                end
            end
        end
    end
end
for i=1:P
    for j=1:nb
G2(i,j)=G(i,i+nb-j+1);                % 从后面看，求解Ｇ1即为书上的G，而G2还？
    end
end
u0=zeros(nb,1);                        
u1=zeros(P,1);      
y0=zeros(P,1);
y=[];                                 % 什么含义y=[] 
u01=zeros(na+2,1);                    % u01是u0的？
u=[];                                 % 
for j=1:na+1
u(j)=0;
end
u=u';
yr=[];yr(1)=0;                        % yr是参考输入
yr1=zeros(P,1);
sy=zeros(P,1);                        % sy代表？
d1=[1,0,0];                           % 取第一行
f=zeros(P,1);                         % f代表？
% w=2;r=2.4;a=0.8;
% w=ones(1,nn+10+1)*5;
% w(101+10:200+10)=-5*ones(1,100);
% w(301+10:400+10)=-5*ones(1,100);
% w(1:10)=zeros(1,10);
w=1;                                  % 定义，为求解参考输入yr
r=2;
a=0.84;
I=eye(P,P);                           
for t=1:400 
    for j=1:P
yr(t+j)=a*yr(t+j-1)+(1-a)*w;          % 对应书上yr的参考轨迹
yr1(j)=yr(t+j);                       % yr1是书上的yd（k+1)
    end
for j=1:P
        sy(j)=S(j,1)*y0(1)+S(j,2)*y0(2)+S(j,3)*y0(3);   % sy（i）代表什么--见求f的公式
    end
f=G2*u0+sy                            % 书上f有具体求解公式，p60
u1=inv(G1'*G1+r*I)*G1'*(yr1-f);       % u1即为△u（k）
ut1=d1*u1;                            
ut=ut1+u01(na+2);                     % ut即为控制增量
u01(1)=u01(2);                        % u0的初始化见63行，代表什么
u01(2)=u01(3);                        
u01(3)=u01(4);                        
u01(4)=ut;                            % 移位，为什么向左移                            
u=[u;ut];                             % u=[u;ut]形式上是什么意思                            
u0(1)=u0(2);                          % u0的初始化见59
u0(2)=u0(3);                             
u0(3)=ut1;                            % u0与u01的区别
y1=B(1)*u(nb+t)+B(2)*u(nb+t-1)+B(3)*u(nb+t-2)+B(4)*u(nb+t-3)-A(2)*y0(2)-A(3)*y0(3); % 次方程又何处的来？
y=[y;y1];                             % 问题同99
y0(3)=y0(2);
y0(2)=y0(1);
y0(1)=y1; 
end 
t=1:400;
plot(t,y);
% grid on 
hold on
t1=410
plot(t1,w);