动态矩阵控制 MATLAB代码
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 i<P 35 S(i,i+1)=1; 36 end 37 if i==P 38 S(P,P)=1; 39 end 40 end 41 R1=eye(M);%构造控制增量加权矩阵R 42 %R=0.1*R1; 43 R=0.0*R1; 44 HT=linspace(1,1,P); 45 H=HT';%构造误差校正向量 46 for i=2:1:P 47 H(i)=0.9; 48 end 49 d1=linspace(0,0,M);%构造向量d 50 d1(1)=1; 51 d=d1; 52 53 %给单位阶跃响应序列a1赋值 54 %{ 55 for i=1:1:step1 56 a11(i)=1*(1-exp(-i*T/60));%采样周期T 57 end 58 figure 59 plot(a11);title('阶跃响应1'); 60 %} 61 t=0:5:500; 62 %t=0:0.5:100; 63 [a11,t0]=step(sys,t);%横轴时间与书上的数量级相差太大，不知原因 64 figure 65 plot(a11);title('阶跃响应step'); 66 %{ 67 for i=1:1:40 68 a11(i)=a1(10*i);%为了与书上的数量级保持一致 乘以系数 求预测模型参数 书上的例5-1 69 end 70 %} 71 %{ 72 a1=1-1.1835*exp(-0.55*t).*sin(6*t+1.4973)-0.18038*exp(-0.25*t).*sin(15.4*t-1.541); 73 figure 74 plot(a1);title('表达式求得的阶跃响应'); 75 %} 76 77 %计算AT A为动态矩阵，AT为其转置 78 for i=1:1:P 79 for j=1:1:M 80 if j<=i 81 A(i,j)=a11(i-j+1); 82 end 83 if j>i&j<M 84 A(i,j)=0; 85 end 86 if i>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 end 102 Uk=Uk'; 103 Yk1=Yk1'; 104 for i=1:1:P 105 Y0(i)=startvalue;%初始化 106 Yp0k1(i)=0; 107 Ycork1(i)=0; 108 end 109 Y0=Y0'; 110 Yp0k1=Yp0k1'; 111 Ycork1=Ycork1'; 112 y1k1=0; 113 daltauk=0;%初始化控制增量 114 uk1=pipestartvalue;% 0 115 uk2=0; 116 yk1=0; 117 yk2=0; 118 for n=1:1:step1+90 119 %x2=(0.98^(n*1))*1+(1-0.98^(n*1))*c;%参考轨迹参数a=0.992 120 %x1=x2; 121 Yrk1(n)=3;%计算参考轨迹yrk1，记录到Yrk1(i) 122 %参考轨迹设为定值3 可以看出PID控制器输出有超调，而DMC可以快速稳定的达到设定值 无超调 123 end 124 Yrk1=Yrk1'; 125 %仿真第一步 126 Yp0k1=Y0; 127 TempYrk1=Yrk1(1:P); 128 daltauk=DT*(TempYrk1-Yp0k1); 129 uk2=uk1+daltauk;%计算当前控制量uk 130 uk1=uk2; 131 Uk(1)=uk1; 132 Yk1(1)=Y0(1);%第一步采样值保存到Yk1 133 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:step1 157 %前15步，由于纯滞后，所以输出为0 158 %if i<=delay 159 %采样 yk1 160 %yk2=-0.01667*yk1+0.125*0.1333*0;%对象离散模型 161 %yk2=-0.2315*yk1+0.6991*0; 162 %end 163 %if i>delay 164 %yk2=-0.01667*yk1+0.125*0.1333*Uk(i-delay);%离散模型的参数 165 %yk2=-0.2315*yk1+0.6991*Uk(i-delay); 166 %end 167 xs1=As*xs0+Bs*uk1; 168 yk2=Cs*xs1; 169 xs0=xs1; 170 171 %if yk2<=startvalue 172 % yk2=startvalue; 173 %end 174 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 end 187 %} 188 if uk2<0 189 uk2=0; 190 end 191 192 uk1=uk2; 193 Uk(i)=uk1; 194 Y1k1=Yp0k1+a*daltauk;%一步预测 195 end 196 Yrklend=Yrk1(1:step1);%整理计时器值，做曲线时使用 197 figure 198 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.54 206 kd=0.62; % 0.2 207 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-1 216 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 end 232 Yrklend=Yrk1(1:step1); 233 x=1.*(1:step1); 234 figure 235 plot(x,Yk1,'b',x,y,'r');%将DMC控制与PID控制的输出值，画在一张表上进行比较 236 title(['预测时域P=',num2str(P)]); 237 figure 238 plot(x,y,'r');title(['采样周期Ts=',num2str(Ts)]);%PID控制器输出曲线
注意参数的选择，尤其是采样周期Ts，控制时域P，一般都是先选定M，再调整P
