1 %预测控制书上的P79例5-1 得到的输出曲线趋近于无穷 不对 不知错误在哪里 pid控制器也是趋近于无穷大 2 %不明白采样周期Ts怎么用，什么意思？？？ 3 %将阶跃响应 离散状态空间模型的采样周期都设为Ts=5 预测步长P=50 M=1都有了很好的效果 4 %所以有两个重要的参数：采样周期Ts 预测步长P 还有M参数的作用，要弄清楚 5 clear all 6 %传递函数模型 7 %{ 8 num=[8611.77]; 9 p1=[1,1.1,36.3025]; 10 p2=[1,0.5,237.2225]; 11 den=conv(p1,p2); 12 sys=tf(num,den); 13 %} 14 sys=tf(0.6,[2400 85 1]); 15 16 Ts=5;%Ts为采样周期 17 delay=0;%延迟时间 即纯滞后模块 18 startvalue=0;%系统初始输出值 19 x1=startvalue; 20 x2=0; 21 c=3;%阶跃值 22 pipestartvalue=0;%管温初始值 23 step1=101;%仿真长度 注意变量名字不能与MATLAB中的函数名相同 否则函数不能再调用 24 %P=50;%效果最好 之前一直不稳定 可能是因为P取得太小 或者是采样周期T保持了一致 25 %M=1; 26 P=50; 27 M=1; 28 Q=eye(P);%构造预测输出误差加权矩阵 29 for i=1:1:delay 30 Q(i,i)=0; 31 end %预测输出误差加权阵，对应纯滞后长度的权值为0 32 S=zeros(P);%构造移位矩阵 33 for i=1:1:P 34 if ii&j
M 87 A(i,j)=a11(i-j+1); 88 end 89 end 90 end 91 AT=A'; 92 %计算DT 93 DT=d*inv(AT*Q*A+R)*AT*Q;%计算得到控制向量DT 94 a=a11(1:P);%计算a列向量 95 %a=a';% % 96 qu1=linspace(0,0,P); 97 qu1(1)=1;%构建取1向量 98 for i=1:1:step1 99 Uk(i)=0;%初始化UK，用来记录控制量100 Yk1(i)=0;%初始化Yk1，用来记录实际仿真输出值101 end102 Uk=Uk';103 Yk1=Yk1';104 for i=1:1:P105 Y0(i)=startvalue;%初始化106 Yp0k1(i)=0;107 Ycork1(i)=0;108 end109 Y0=Y0';110 Yp0k1=Yp0k1';111 Ycork1=Ycork1';112 y1k1=0;113 daltauk=0;%初始化控制增量114 uk1=pipestartvalue;% 0115 uk2=0;116 yk1=0;117 yk2=0;118 for n=1:1:step1+90119 %x2=(0.98^(n*1))*1+(1-0.98^(n*1))*c;%参考轨迹参数a=0.992120 %x1=x2;121 Yrk1(n)=3;%计算参考轨迹yrk1，记录到Yrk1(i)122 %参考轨迹设为定值3 可以看出PID控制器输出有超调，而DMC可以快速稳定的达到设定值 无超调123 end124 Yrk1=Yrk1';125 %仿真第一步126 Yp0k1=Y0;127 TempYrk1=Yrk1(1:P);128 daltauk=DT*(TempYrk1-Yp0k1);129 uk2=uk1+daltauk;%计算当前控制量uk130 uk1=uk2;131 Uk(1)=uk1;132 Yk1(1)=Y0(1);%第一步采样值保存到Yk1133 yk1=Y0(1);%第一步不用移位操作，直接取实际系统的输出值作为预测值134 Y1k1=Yp0k1+a*daltauk;%一步预测135 %{136 %Ts=5;137 As=[ -1.6,-17.13,-2.18,-8.41;16,0,0, 0; 0,8,0,0;0,0,8,0];%状态空间方程系数138 As=eye(4)+As*Ts;139 Bs=[4;0;0;0];140 Bs=Bs*Ts;141 Cs=[0,0,0,2.102];142 xs0=[0;0;0;0];143 xs1=[0;0;0;0];144 %}145 146 %Ts=5;147 As=[-0.03542,-0.02667;0.01563,0];148 As=eye(2)+As*Ts;149 Bs=[0.125;0];150 Bs=Bs*Ts;151 Cs=[0,0.128];152 xs0=[0;0];153 xs1=[0;0];154 155 %第二步及其以后的仿真156 for i=2:1:step1157 %前15步，由于纯滞后，所以输出为0158 %if i<=delay159 %采样 yk1160 %yk2=-0.01667*yk1+0.125*0.1333*0;%对象离散模型161 %yk2=-0.2315*yk1+0.6991*0;162 %end163 %if i>delay164 %yk2=-0.01667*yk1+0.125*0.1333*Uk(i-delay);%离散模型的参数165 %yk2=-0.2315*yk1+0.6991*Uk(i-delay);166 %end167 xs1=As*xs0+Bs*uk1;168 yk2=Cs*xs1;169 xs0=xs1;170 171 %if yk2<=startvalue172 % yk2=startvalue;173 %end174 yk1=yk2;175 Yk1(i)=yk1;%采样结束，并保存到Yk1中176 Y0k1=Y1k1;177 y1k1=qu1*Y0k1;%计算y1k1，即Y0k1的第一个元素178 Ycork1=Y0k1+H*(yk2-y1k1);%计算校正预测值179 Yp0k1=S*Ycork1;%移位，计算初始预测值180 TempYrk2=Yrk1(i:i+P-1);181 daltauk=DT*(TempYrk2-Yp0k1);182 uk2=uk1+daltauk;183 %{184 if uk2>4;185 uk2=4;186 end187 %}188 if uk2<0189 uk2=0;190 end191 192 uk1=uk2;193 Uk(i)=uk1;194 Y1k1=Yp0k1+a*daltauk;%一步预测195 end196 Yrklend=Yrk1(1:step1);%整理计时器值，做曲线时使用197 figure198 x=Ts*(1:step1);199 plot(t,Yrklend,t,Yk1);%将实际输出与期望输出两条曲线画在一张图中，要保证二者矢量长度相同200 title(['预测时域P=',num2str(P)]);201 202 %以下为增量式PID控制算法203 y(1)=0; 204 kp=0.35; % 0.4效果会好一些 曲线形式相同205 ki=0.1; % 0.54206 kd=0.62; % 0.2207 actual=0;208 e=0;209 e1=0;210 e2=0;211 uk0_pid=0;212 x0=[0;0];213 x1=[0;0];214 215 for i=1:1:step1-1216 e=Yrk1(i)-actual;217 %e=set-actual;218 increase=kp*(e-e1)+ki*(e)+kd*(e-2*e1+e2);219 uk_pid=uk0_pid+increase;220 %y(i+1)=-0.2315*y(i)+0.6991*uk_pid;%离散模型参数 离散模型参数可由传递函数得到ss(system)221 222 x1=As*x0+Bs*uk_pid;223 y(i+1)=Cs*x1;224 x0=x1;225 226 e1=e;227 e2=e1;228 actual=y(i+1);229 uk0_pid=uk_pid;230 nums(i)=e; 231 end232 Yrklend=Yrk1(1:step1);233 x=1.*(1:step1);234 figure235 plot(x,Yk1,'b',x,y,'r');%将DMC控制与PID控制的输出值，画在一张表上进行比较236 title(['预测时域P=',num2str(P)]);237 figure238 plot(x,y,'r');title(['采样周期Ts=',num2str(Ts)]);%PID控制器输出曲线 注意参数的选择，尤其是采样周期Ts，控制时域P，一般都是先选定M，再调整P