fangzhen.m

广义预测控制源程序(单输入单输出)

clear;
mode=input('请输入模式（1为定常，2为自校正）：');

yr=1;
arfa=0.2;
k=1;



if mode==1
na=input('请输入na:');
nb=input('请输入nb:');
nb=nb+1;
A(1)=input('请输入A:')
for i=2:1:na+1
    A(i)=input('     ');
end
A=fliplr(A);
B(1)=input('请输入B:');
for i=2:1:nb
    B(i)=input('    ');
end
B=fliplr(B);

N=input('请输入N:');
Nu=input('请输入Nu:');
lamda=input('请输入lamda:');

y=[zeros(1,na),0]';u=zeros(1,nb)';

ta=[-1,1];                    %求G
A1=conv(A,ta);
E=[1];
EE=[zeros(1,N-1),E];
F=-A1(1:end-1);
FF=F;
EB=conv(E,B);
G=[EB(end),zeros(1,Nu-1)];
H=EB(1:end-1);
for j=2:1:N
    zj=[1,zeros(1,j-1)];
    e=F(end);
    E=polyadd(E,e*zj);
    f=polyadd(F,-e*A1);
    F=f(1:end-1);
    EE(j,:)=[zeros(1,N-j),E];
    FF(j,:)=F;
    EB=conv(E,B);
    if j<=Nu
        G(j,:)=[EB(end-j+1:end),zeros(1,Nu-j)];  
    else
        G(j,:)=EB(end-j+1:end-j+Nu);
    end
    H(j,:)=EB(1:end-j);
end
G=inv(G'*G+lamda*eye(Nu,Nu))*G';%%%%%%%%%%%%%



for i=1:k:100                             %求U
    %kesi=-0.05+0.1*rand(1);
    yk=(-A)*y+B*u;%+kesi;
    yy=[y(1:na)',yk]';
    for l=1:1:(na-1)
        y(l)=y(l+1);
    end
    y(na)=yk;
    yd=yk;
    yK(i)=yk;
    for j=1:1:N
        yd=arfa*yd+(1-arfa)*yr;
        ydk(j)=yd;
        y0=FF(j,:)*yy+H(j,:)*(u(2)-u(1));
        y0k(j)=y0;
    end
    drtu=G*(ydk-y0k)';
    for m=1:1:nb-1
        u(m)=u(m+1);
    end
    u(nb)=drtu(1)+u(nb);%%%%%%%%%%%%%%%%%%%%%%%%
end
i=1:k:100;                     %画图
 plot(i,yr,'r',i,yK,'g');%%%%%%%%%%%%%%%%%%%
 
 
elseif mode==2
na=input('请输入na:');
nb=input('请输入nb:');
nb=nb+1;

N=input('请输入N:');
Nu=input('请输入Nu:');
lamda=input('请输入lamda:');

A=[ones(1,na)/1000,1];
B=[ones(1,nb)/1000];
c0=[A(1:na),B]';
p0=10^6*eye(na+nb,na+nb);
EEE=0.000000005;%相对误差E=0.000000005
e2=ones(na+nb,1);
y=[zeros(1,na)/100,0]';u=zeros(1,nb)';

for i=1:k:100
     %kesi=-0.05+0.1*rand(1);
     yk=(-[0.7,-1.5,1])*y+[0.5,1]*u;%+kesi;采样
     yy=[y(1:na)',yk]';
     for l=1:1:(na-1)
         y(l)=y(l+1);
     end
     y(na)=yk;
     yd=yk;
     yK(i)=yk;
     

     %if e2>EEE
        h1=[-yy(1:2)',u']'; %辨识
        x=h1'*p0*h1+1; 
        x1=inv(x);
        k1=p0*h1*x1;
        d1=yk-h1'*c0; 
        c1=c0+k1*d1;
        e1=c1-c0;%求参数当前值与上一次的值的差值
        e2=e1./c0;%求参数的相对变化
        c0=c1;
        p1=p0-k1*k1'*[h1'*p0*h1+1];
        p0=p1;
        A=[c1(1:na)',1];B=c1(na+1:na+nb)';%%%%%%%%%%%%%%%%%%%%%%%%%%
      
        
        ta=[-1,1];                        %求G
        A1=conv(A,ta);
        E=[1];
        EE=[zeros(1,N-1),E];
        F=-A1(1:end-1);
        FF=F;
        EB=conv(E,B);
        G=[EB(end),zeros(1,Nu-1)];
        H=EB(1:end-1);
        for j=2:1:N
            zj=[1,zeros(1,j-1)];
            e=F(end);
            E=polyadd(E,e*zj);
            f=polyadd(F,-e*A1);
            F=f(1:end-1);
            EE(j,:)=[zeros(1,N-j),E];
            FF(j,:)=F;
            EB=conv(E,B);
            if j<=Nu
               G(j,:)=[EB(end-j+1:end),zeros(1,Nu-j)];  
            else
               G(j,:)=EB(end-j+1:end-j+Nu);
            end
            H(j,:)=EB(1:end-j);
        end
        G=inv(G'*G+lamda*eye(Nu,Nu))*G';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % end
 
     
     for j=1:1:N                              %求U
         yd=arfa*yd+(1-arfa)*yr;
         ydk(j)=yd;
         y0=FF(j,:)*yy+H(j,:)*(u(2)-u(1));
         y0k(j)=y0;
     end
     drtu=G*(ydk-y0k)';
     for m=1:1:nb-1
         u(m)=u(m+1);
     end
     u(nb)=drtu(1)+u(nb);
end                                         %%%%%%%%%%%%%%%%%%%%%%


i=1:k:100;                         %画图
plot(i,yr,'r',i,yK,'g');%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end