循环流化床床温控制系统(GPC)-单模型

%一次风量调节床温closeall;clc;clear;%对象本身num0=[-7.8-1.3];den0=conv([1151],[1151]);sys0=tf(num0,den0);z=tf('z',-1);sys00=c2d(sys0,10,'tustin');T0=840;V10=71.3;B0=sys00.num{1};A0=sys00.den{1};G00=filt(B0,A0);Bb=G00.num{1};Aa=G00.den{1};hf=60;%典型工况170%num1=[-1.2-0.3];den1=conv([1751],[1701]);sys1=tf(num1,den1);z=tf('z',-1);sys11=c2d(sys1,10,'tustin');B11=sys11.num{1};A11=sys11.den{1};G01=filt(B11,A11);B1=G01.num{1};A1=G01.den{1};%实际的初值p1=8;N11=1;N12=41;y1(1:200)=0;u1(1:200)=0;dU1=zeros(p1,1);dUK1(1:1200)=0;du=zeros(length(B1),1)';dy=zeros(length(A1)-1,1)';fy=zeros(length(du)+length(dy),1)';zeta0=[A1(2:end)B1(1:end)]';zeta1=zeros(length(du)+length(dy),1);P0=10000*eye(length(du)+length(dy));t=[1:1200];%仿真时间fort=1:1200;ift=200;yr(t)=0;elseif(t700&t=1200)yr(t)=2;elseyr(t)=4;%设定值扰动endendfort=200:1199;%解丢潘图方程[E1,F1]=diupantu(A1,N11,N12);fori=1:N12-N11;CF1(i,:)=conv(B1,E1(:,i));%多项式相乘end%将G0的每行的前j列取出放入G中fori=1:N12-N11;ifi=p1;forj=1:i;G11(i,i-j+1)=CF1(i,j);endelseforj=1:p1;G11(i,p1-j+1)=CF1(i,j);endendendG12(1,:)=G11(1,:);fori=2:N12-N11;G12(i,:)=G11(i,:)-G11(i-1,:);endG13(1,:)=G12(1,:);fori=2:N12-N11;G13(i,:)=G12(i,:)-G12(i-1,:);end%将G0的每行的第j列之后的元素取出放入H中fori=1:N12-N11;forj=i+1:i+length(B1)-1;%从G0的第i行第i+1列开始取，只取非零项H1(i,j-i)=CF1(i,j);%在H中的列序号往前移动iendendK11=12;K12=15;K13=18;%设定PID参数%设定权值Q1=eye(p1)*2;D11=inv(K11*G11.'*G11+K12*G12.'*G12+K13*G13.'*G13+Q1)*G11.'*K11;D12=inv(K11*G11.'*G11+K12*G12.'*G12+K13*G13.'*G13+Q1)*G12.'*K12;D13=inv(K11*G12.'*G12+K12*G11.'*G11+K13*G13.'*G13+Q1)*G13.'*K13;v=0.9;w1=cankaoguiji(v,yr,t,N11,N12,y1);%参考轨迹%求最优控制f1=H1*dUK1(t:-1:t-length(B1)+2).'+F1.'*y1(t:-1:t-length(A1)+1).';d11w(1)=w1(1);d11f(1)=f1(1);fori=2:N12-N11;d11w(i)=w1(i)-w1(i-1);d11f(i)=f1(i)-f1(i-1);endd12w(1)=d11w(1);d12f(1)=d11f(1);fori=2:N12-N11;d12w(i)=d11w(i)-d11w(i-1);d12f(i)=d11f(i)-d11f(i-1);enddU1=D11*(d11w.'-d11f.')+D12*(w1-f1)+D13*(d12w.'-d12f.');u1(t+1)=u1(t)+dU1(1)+0*wgn(1,1,1);%将控制施加于实际对象y1(t+1)=-1*Aa(2:end)*y1(t:-1:t+2-length(Aa)).'+Bb*u1(t+1:-1:t+2-length(Bb)).';dUK1(t+1)=u1(t+1)-u1(t);%用递推最小二乘法在线辨识du=[u1(t:-1:t-length(du)+1)-u1(t-1:-1:t-length(du))];dy=[y1(t:-1:t-length(dy)+1)-y1(t-1:-1:t-length(dy))];fy=[-1*dy(1:end)du(1:end)];K=P0*fy'*inv([fy*P0*fy'+0.9999]);P1=1/0.98*[eye(length(P0))-K*fy]*P0;zeta1=zeta0+K*[y1(t+1)-y1(t)-fy*zeta0];A1=[1;zeta1(1:length(dy))]';B1=[zeta1(length(dy)+1:end)]';zeta0=zeta1;endfort=1:1200;Tb1(t)=y1(t)+T0;V1(t)=u1(t)+V10;Tbr(t)=yr(t)+T0;time(t)=t*10/60;end%绘图figuresubplot(211)plot(time,Tbr,'r--');holdonplot(time,Tb1,'b');gridtitle('循环流化床床温控制系统（GPC）');xlabel('仿真时间time(min)');ylabel('Tb(℃)');legend('设定温度','床温')subplot(212)plot(time,V1,'b');%控制量gridtitle('一次风量');xlabel('仿真时间time(min)');ylabel('V1(kg/s)');