广义预测控制（GPC）GPC算法仿真被控对象模型动态矩阵控制算法的编程原理(1)设置GPC参数,例如采样周期，预测时域，控制时域，截断步长等。(2)建立系统阶跃响应模型(3)设置初始时刻参数，例如系统的初始时刻值，柔化系数等。(4)计算参考轨迹(5)计算控制作用增量(6)实施GPC控制(7)输出结果，绘制曲线GPC算法：1.初选控制参数：Q、R、P、M、ysp、?、?(z-1)2.采集输入、输出样本{?u(k),?y(k)}3.用RLS算法估计参数4.递推求解Diophantine方程，得到5.计算F(k)6.在线计算控制器参数dT7.得到控制增量?u(k)和控制输入u(k)＝u(k－1)＋?u(k)8.k+1?k，进入下一周期预测计算和滚动优化GPC程序：%Clarke广义预测控制（C=1）（对象参数已知）%N1=d、N、Nu取不同的值clearall；closeall;a=cell(1,2);b=cell(1,2);c=cell(1,1);d=cell(1,1);%对象参数symsk;k=length(k);if(0<=k<=150)a=[10.9234];elseifb=[7.24020.9485];c=1;d=1;(150<k<=300)a=[10.8981];elseifb=[9.99010.14142];c=1;d=1;(300<k<=450)b=[9.60410.34067];c=1;d=1;a=[10.8838];else(450<k<=600)a=[10.9234];b=[7.24020.9485];c=1;d=1;endna=length(a)-1;b=[zeros(1,d-1)b];nb=length(b)-1;%na、nb为多项式A、B阶次（因d！=1,对b添0）aa=conv(a,[1-1]);naa=na+1;%aa的阶次N1=d;N=15;Nu=5;%最小输出长度、预测长度、控制长度gamma=1*eye(Nu);alpha=0.11;%控制加权矩阵、输出柔化系数L=600;%控制步数uk=zeros(d+nb,1);%输入初值：uk(i)表示u(k-i)duk=zeros(d+nb,1);%控制增量初值yk=zeros(naa,1);%输出初值w=10*[ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4+d,1)];%设定值xi=sqrt(0.01)*randn(L,1);%白噪声序列%求解多步Diophantine方程并构建F1、F2、G[E,F,G]=multidiophantine(aa,b,c,N);G=G(N1:N,:);F1=zeros(N-N1+1,Nu);F2=zeros(N-N1+1,nb);fori=1:N-N1+1forj=1:min(i,Nu);F1(i,j)=F(i+N1-1,i+N1-1-j+1);endforj=1:nb;endF2(i,j)=F(i+N1-1,i+N1-1+j);endfork=1:Lif(1<=k<=150)time(k)=k;a=[10.9234];b=[7.24020.9485];c=1;d=1;y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)dUk=duk(1:nb);%构建向量△U(k-j)elseif(150<k<=300)time(k)=k;a=[10.8981];b=[9.99010.14142];c=1;d=1;y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)dUk=duk(1:nb);%构建向量△U(k-j)elseif(300<k<=450)time(k)=k;a=[10.8838];b=[9.60410.34067];c=1;d=1;y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k);yk(1:na)];%构建向量Y(k)dUk=duk(1:nb);%构建向量△U(k-j)else(450<k<=L)time(k)=k;a=[10.9234];b=[7.24020.9485];c=1;d=1;y(k)=-aa(2:naa+1)*yk+b*duk(1:nb+1)+xi(k);%采集输出数据Yk=[y(k