我这人有个例子,你参考一下吧- clc
- clear
- %DMC预测控制在加热系统温度控制中的应用仿真程序
- startvalue=0;%系统初始输出值
- x1=startvalue;
- x2=0;
- c=3;%阶跃值
- pipestartvalue=0;%管温初始值
- step=250;%仿真长度
- P=80;%预测时域长度赋值
- M=2;%控制时域长度赋值
- Q=eye(P);%构造预测输出误差加权阵
- for i=1:1:15
- Q(i,i)=0;
- end%预测输出误差加权阵.对应纯滞后长度的权值取0
- S=zeros(P);%构造移位矩阵
- for i=1:1:P
- if i<P
- S(i,i+1)=1;
- end
- if i==P
- S(P,P)=1;
- end
- end
- Rl=eye(M);%构造控制增量加权矩阵R
- R=0.1*Rl;
- HT=linspace(1,1,P);
- H=HT;%构造误差校正向量
- for i=2:1:P
- H(i)=0.9;
- end
- d1=linspace(0,0,M);%构造向量d
- d1(1)=1;
- d=d1;
- %给阶跃响应序列al赋值
- a1=[2.1,2,2.2,2.19,2.18,2.17,2.16,2.15,2.19,2.20,2.11,2.2,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20,2.11,2.12,2.13,2.14,2.15,2.16,2.17,2.18,2.19,2.20];
- a1=a1*1;
- %计算AT
- for i=1:1:P
- for j=1:1:M
- if j<=i
- A(i,j)=a1(i-j+1);
- end
- if j>i&j<=M
- A(i,j)=0;
- end;
- if i>=M
- A(i,j)=a1(i-j+1);
- end
- end
- end
- AT=A';
- %计算DT
- DT=d*inv(AT*Q*A+R)*AT*Q;%计算得到行向量DT(1xP,1x80)
- a=a1(1:P);%计算a列向量(80x1,Pxl)
- qul=linspace(0,0,P);
- qul(1)=1;%构建取1向量
- for i=1:1:step
- Uk(i)=0;%初始化Uk,用来记录控制量
- Yk1(i)=0;%初始化Ykl,用来记录实际仿真输出值
- t(i)=i;%计时器
- end
- Uk=Uk';
- Yk1=Yk1';
- for i=1:1:P
- Y0(i)=startvalue;%初始化YO(i),YpOkl(i),Ycorkl(i)
- Yp0k1(i)=0;
- Ycork1(i)=0;
- end
- Y0=Y0';
- Yp0k1=Yp0k1';
- Ycork1=Ycork1';
- y1k1=0;
- daltauk=0;%初始化控制增量daltauk
- uk1=pipestartvalue;
- uk2=0;
- yk1=0;
- yk2=0;
- for n=1:1:step+90
- x2=(0.992^n)*x1+(1-0.992^n)*c;%参考轨迹参数a=0.992
- x1=x2;
- Yrk1(n)=1;%x2;%计算参考轨迹yrkl,记录到Yrkl(i)
- end;
- Yrk1=Yrk1;
- %仿真第一步
- Yp0k1=Y0;
- TempYrk1=Yrk1(1:P);
- daltauk=DT*(TempYrk1'-Yp0k1);
- uk2=uk1+daltauk;%计算控制量uk
- uk1=uk2;
- Uk(1)=uk1;
- Yk1(1)=Y0(1);%第一步采样值保存到Ykl;
- %第一步不用移位操作,直接取实际系统输出值作预测值
- yk1=Y0(1);
- Y1k1=Yp0k1+a'*daltauk;%一步预测
- %第二步及其以后步仿真
- for i=2:1:step
- %前15步,由于纯滞后,所以输出值为0
- if i<=15
- %采样ykl
- yk2=0.9689*yk1+0.004644*0;%对象离散模型
- end
- if i>15
- yk2=0.9689*yk1+0.004644*Uk(i-15);
- end
- if yk2<=startvalue
- yk2=startvalue;
- end
- yk1=yk2;
- Yk1(i)=yk1;%采样结束,并保存到Yk1中
- YOk1=Y1k1;%开头大写表示向量,小写表示数值
- y1k1=qul*YOk1;%计算ylkl,既就是YOU的第一个元素
- Ycork1=YOk1+H'*(yk2-y1k1);%计算校正预测值
- YpOk1=S*Ycork1;%移位;计算初始预测值
- TempYrk2=Yrk1(i:i+P-1);
- daltauk=DT*(TempYrk2'-YpOk1);%计算控制增量
- uk2=uk1+daltauk;%计算控制摄uk
- if uk2>75
- uk2=75;
- end
- if uk2<0
- uk2=0;
- end
- uk1=uk2;
- Uk(i)=uk1;
- Y1k1=YpOk1+a'*daltauk;%一步预测
- end
- Yrk1end=Yrk1(1:step);%整理计时器值,做曲线时使用
- %仿真程序结束
- [m_1,n_1]=size(Y1k1);
- ttt=1:1:m_1;
- figure(1);
- plot(ttt,Y1k1)
- grid on
- figure(2);
- plot(Yk1)
- grid on
复制代码 |