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

轉載于:https://www.cnblogs.com/1987-05-04/p/6957020.html

總結

以上是生活随笔為你收集整理的动态矩阵控制 MATLAB代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。