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控制器输出曲线