mgpc 多变量预测控制

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（2001）;
t3=zeros（3001）;
d1=input（‘（2*n+1）阶方阵P的形式:0-对角阵1-方阵:‘）;
d2=input（‘（2*n+1）阶方阵P的初始值:0-自动赋值1e+51-键盘输入:‘）;
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（n1）;u1=zeros（m1）;uu2=zeros（n1）;u2=zeros（m1）;
yy1=zeros（n1）;y11=zeros（n1）;yy2=zeros（n1）;y12=zeros（n1）;
Q1=zeros（2*n+11）;Q1（11）=1;Q1（n+11）=1;Q1（2*n+11）=1;
Q2=Q1;
%产生时间为T，幅值为1的阶跃信号的给定值
T=300;
t=1:T;
yr0（t1）=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（n1）-0.6787*yy1（n-11）+0.0202*yy1（n-21）-0.0012*uu1（n1）+0.2167*uu1（n-11）-0.4678*uu1（n-21）+0.2543*uu1（n-31）+0.0003*uu2（n-11）+0.0396*uu2（n-21）+0.0301*uu2（n-31）-0.0025*uu2（n-41）;
    y2=2.208*yy2（n1）-1.5457*yy2（n-11）+0.3339*yy2（n-21）-0.0008*uu1（n-11）-0.0061*uu1（n-21）+0.0091*uu1（n-31）-0.0027*uu1（n-41）+0.0052*uu2（n-11）-0.0007*uu2（n-21）-0.0033*uu2（n-31）;
    %（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（i1）=yy1（i+11）;yy2（i1）=yy2（i+11）;
    end
    yy1（n1）=y1;yy2（n1）=y2;
    yyy1=[yyy1;（y1-1.2970）*210];yyy2=[yyy2;（y2+0.7014）*18];%保存各k时刻的nm个输出量以便绘图
    %根据最小二乘公式，由y1（k）和y2（k）计算G11，G12，G21和G22阵的各元素值g0g1…gn
    for i=1:n
        X（1i）=uu1（i1）;X（1i+n）=uu2（i1）;
    end
    X（12*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）;
    %根据元素值g0g1…gn求G11G12G21和G22阵
    for j=1:m
        for i=n:-1:j
            i1=n-i+j;
            G11（i1j）=Q1（i1）;G12（i1j）=Q1（i+n1）;
            G21（i1j）=Q2（i1）;G22（i1j）=Q2（i+n1）;
        end
    end
    %for i=1:nG12（i2）=0;G21（i2）=0;end
    %求nn维y01y02向量（y11和y12为上一时刻的y01和y02向量）
    e1=y1-y11（11）;e2