動態矩陣控制 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