gpc.m

多变量广义预测的控制算法MATLAB代码

clear
disp('多变量系统的隐式广义预测控制算法的研究(pt70)')
disp('广义预测控制算法初始值')
nn=input('时域长度nn=');
n=input('预测长度n=');
m=input('控制长度m=');
t0=input('控制加权系数t0=');
a=input('柔化系数a=');
disp('最小二乘公式初始值')
t1=1;%input('遗忘因子t1=');
t2=zeros(200,1);
t3=zeros(300,1);
d1=input('(2*n+1)阶方阵P的形式:0-对角阵,1-方阵:');
d2=input('(2*n+1)阶方阵P的初始值:0-自动赋值1e+5,1-键盘输入:');
if(d1==1)
    if(d2==1)
        P=input('在方括号[]中，输入(2*n+1)阶方阵P的值:');
    else
        P=(1e+5)*ones(2*n+1);
    end
else
    if(d2==1)
        PP=input('在方括号[]中，输入(2*n+1)阶对角方阵P对角线上的值:');
        P=diag(PP);
    else
        P=(1e+5)*eye(2*n+1);
    end
end
%参数初始值
uuu1=0;yyy1=0;uuu2=0;yyy2=0;
uu1=zeros(n,1);u1=zeros(m,1);uu2=zeros(n,1);u2=zeros(m,1);
yy1=zeros(n,1);y11=zeros(n,1);yy2=zeros(n,1);y12=zeros(n,1);
Q1=zeros(2*n+1,1);Q1(1,1)=1;Q1(n+1,1)=1;Q1(2*n+1,1)=1;
Q2=Q1;
%产生时间为T，幅值为1的阶跃信号的给定值
T=300;
t=1:T;
yr0(t,1)=heaviside(t);
d3=input('输出曲线是否去掉前100步:0-不,1-去掉:');
nm=length(t);%确定循环次数
for ij=2:nm
    yr1=yr0(ij);
    yr2=yr1;
    %根据系统模型，计算K时刻的输出值y1(k)和y2(k)
    y1=1.629*yy1(n,1)-0.6787*yy1(n-1,1)+0.0202*yy1(n-2,1)-0.0012*uu1(n,1)+0.2167*uu1(n-1,1)-0.4678*uu1(n-2,1)+0.2543*uu1(n-3,1)+0.0003*uu2(n-1,1)+0.0396*uu2(n-2,1)+0.0301*uu2(n-3,1)-0.0025*uu2(n-4,1);
    y2=2.208*yy2(n,1)-1.5457*yy2(n-1,1)+0.3339*yy2(n-2,1)-0.0008*uu1(n-1,1)-0.0061*uu1(n-2,1)+0.0091*uu1(n-3,1)-0.0027*uu1(n-4,1)+0.0052*uu2(n-1,1)-0.0007*uu2(n-2,1)-0.0033*uu2(n-3,1);
    %(nn=n=6;m=2;t=0.8;a=0.3;t1=1)
    %保存k时刻及以前的n个输出值y(k),y(k-1),…,y(k-n),以供模型运算
    for i=1:n-1
        yy1(i,1)=yy1(i+1,1);yy2(i,1)=yy2(i+1,1);
    end
    yy1(n,1)=y1;yy2(n,1)=y2;
    yyy1=[yyy1;(y1-1.2970)*210];yyy2=[yyy2;(y2+0.7014)*18];%保存各k时刻的nm个输出量以便绘图
    %根据最小二乘公式，由y1(k)和y2(k)计算G11，G12，G21和G22阵的各元素值g0,g1,…,gn
    for i=1:n
        X(1,i)=uu1(i,1);X(1,i+n)=uu2(i,1);
    end
    X(1,2*n+1)=1;
    K=P*X'*inv(t1+X*P*X');
    P=(eye(2*n+1)-K*X)*P/t1;
    Q1=Q1+K*(y1-X*Q1);Q2=Q2+K*(y2-X*Q2);
    %根据元素值g0,g1,…,gn,求G11,G12,G21和G22阵
    for j=1:m
        for i=n:-1:j
            i1=n-i+j;
            G11(i1,j)=Q1(i,1);G12(i1,j)=Q1(i+n,1);
            G21(i1,j)=Q2(i,1);G22(i1,j)=Q2(i+n,1);
        end
    end
    %for i=1:n,G12(i,2)=0;G21(i,2)=0;end
    %求nn维y01,y02向量(y11和y12为上一时刻的y01和y02向量)
    e1=y1-y11(1,1);e2=y2-y12(1,1);
    y01(1:nn-1,1)=y11(2:nn,1);y01(nn,1)=y11(nn,1);y01=y01+e1;
    y02(1:nn-1,1)=y12(2:nn,1);y02(nn,1)=y12(nn,1);y02=y02+e2;
    for i=1:n
        y11(i,1)=y01(i,1)+u1(1,1)*Q1(n+1-i,1)+u2(1,1)*Q1(2*n+1-i,1);
        y12(i,1)=y02(i,1)+u2(1,1)*Q2(n+1-i,1)+u1(1,1)*Q2(2*n+1-i,1);
    end
    for i=n+1:nn
        y11(i,1)=y11(n,1);y12(i,1)=y12(n,1);
    end
    %根据y01,y02求n维f1,f2向量
    f1(1:n,1)=y01(1:n,1);
    f2(1:n,1)=y02(1:n,1);
    %保存k时刻及以后的n个控制增量,以供模型运算
    for i=1:n-1
        uu1(i,1)=uu1(i+1,1);uu2(i,1)=uu2(i+1,1);
    end
    uu1(n,1)=u1(1,1);uu2(n,1)=u2(1,1);
    uuu1=[uuu1;u1(1,1)];uuu2=[uuu2;u2(1,1)];%保存各k时刻的nm个控制增量以便绘图
    %由当前k时刻的输出值y(k)和给定值yr,求k时刻以后的n个参考轨迹w(k+1),…,w(k+n)
    w1=a*y1+(1-a)*yr1;w2=a*y2+(1-a)*yr2;
    for i=2:n
        w1=[w1;a^i*y1+(1-a^i)*yr1];
        w2=[w2;a^i*y2+(1-a^i)*yr2];
    end
    %计算k时刻以及以后的m个控制增量Du(k),…,Du(k+m)
    u1=inv(G11'*G11+t0*eye(m))*G11'*(w1-G12*u2-f1);
    u2=inv(G22'*G22+t0*eye(m))*G22'*(w2-G21*u1-f2);
    %控制量限幅
    if(u1(1,1)>1)u1(1,1)=1;end
    if(u1(1,1)<-1)u1(1,1)=-1;end
    if(u2(1,1)>1)u2(1,1)=1;end
    if(u2(1,1)<-1)u2(1,1)=-1;end
end
%绘制给定值，输出值和控制增量曲线
if(d3==1)
    %绘制去掉前100步的给定值,输出值和控制增量曲线
    yyyy1(1:(T-100),1)=yyy1(101:T,1);uuuu1(1:(T-100),1)=uuu1(101:T,1);
    yyyy2(1:(T-100),1)=yyy2(101:T,1);uuuu2(1:(T-100),1)=uuu2(101:T,1);
    for i=1:T
        t3(i,1)=i;
    end
    t2(1:(T-100),1)=t3(101:T,1)-100;yr01(1:(T-100),1)=yr0(101:T,1);
    subplot(2,2,1):plot(t2,(yr01*210),t2,yyyy1);
    axis([0,nm-100,0,300]);xlabel('t');ylabel('yr1,y1')
    subplot(2,2,3):plot(t2,uuuu1);
    axis([0,nm-100,-1.5,1.5]);xlabel('t');ylabel('u1')
    subplot(2,2,2):plot(t2,(yr01*18),t2,yyyy2);
    axis([0,nm-100,0,20]);xlabel('t');ylabel('yr2,y2')
    subplot(2,2,4):plot(t2,uuuu2);
    axis([0,nm-100,-1.5,1.5]);xlabel('t');ylabel('u2')
else
    %绘制完整的给定值,输出值和控制增量曲线
    subplot(2,2,1):plot(t,(yr0*210),t,yyy1);
    axis([0,nm,0,300]);xlabel('t');ylabel('yr1,y1')
    subplot(2,2,3):plot(t,uuu1);
    axis([0,nm,-1.5,1.5]);xlabel('t');ylabel('u1')
    subplot(2,2,2):plot(t,(yr0*18),t,yyy2);
    axis([0,nm,0,20]);xlabel('t');ylabel('yr2,y2')
    subplot(2,2,4):plot(t,uuu2);
    axis([0,nm,-1.5,1.5]);xlabel('t');ylabel('u2')
end