控制系统数字仿真

现代工程控制理论实验报告实验名称：控制系统数字仿真技术实验时间： 2015/5/3目录一、实验目的 (2)二、实验容 (3)三、实验原理 (3)四、实验方案 (6)1、分别离散法； (6)2、整体离散法； (7)3、欧拉法 (9)4、梯形法 (9)5、龙格——库塔法 (10)五、实验结论 (11)小结： (14)一、实验目的1、探究多阶系统状态空间方程的求解；2、探究多种控制系统数字仿真方法并对之进行精度比较；二、 实验容1、 对上面的系统进行仿真，运用分别离散法进行分析；2、 对上面的系统进行仿真，运用整体离散法进行分析；3、 对上面的系统进行仿真，运用欧拉法进行分析；4、 对上面的系统进行仿真，运用梯形法进行分析；5、 对上面的系统进行仿真，运用龙泽——库塔法进行分析；6、 对上面的几种方法进行总计比较，对他们的控制精度分别进行分析比较； 三、 实验原理1、 控制系统状态空间方程整体离散法的求解；控制系统的传递函数一般为x Ax Bu Y Cx Du⋅=+=+有两种控制框图简化形式如下： KI 控制器可以用框图表示如下：惯性环节表示如下：高阶系统(s)(1)nKG T =+的框图如下对于上面的框图可以简写传递函数x Ax Bu Y Cx Du⋅=+=+根据各环节间的关系可以列写出式子中出现的系数A 、B 、C 和D ，下面进行整体离散法求传递函数的推导00()0...*()...()(t)(0)...*(t)(0)(t)(0)()(0)At At At At At tt Attt A AT tATA At ttAtAtA AtA t x Ax Bu e ex e Ax e Bud e x dt Bue dt dt e x Bue dtex x Bue d e x x e e Bued x xe Bue d t KT x kT x e ττττττττ⋅-⋅-----------=+=+=⇒=⇒=+=+⇒=+==⎰⎰⎰⎰⎰⎰①①得②③③得令()0(1)(1)[(1)]0(1)[(1)]0...(1)[(1)](0)...*(1)()(1)T (1)()()(),kTA kTA kT k TA k TA k T AT k TATA k T kTTTAT At AT At ATBue d t K T x k T x e Bue d e x k ex k Bue d k tx k e x k e Budt e x k e Bdt u k eτττττττ⋅-+⋅++-++-+=++=+-+-=+-=+=+=+⋅Φ=⎰⎰⎰⎰⎰④令⑤⑤④得令令0(1)()(1)TAt m m e Bdtx k x k x k Φ=+=Φ⋅+Φ⋅+⎰得这样，如果知道系数，就可以知道高阶系统的传递函数和状态空间方程。

2、 在控制系统的每一个环节都加一个采样开关，构成分别离散法求解系统的状态空间方程； 采样开关其实是一个零阶保持器(t)(kT)(1)u u kT t k T =≤<+比例环节：(s)(1)*(1)(s)x kp x k kp u k u =⇒+=+ 积分环节：(s)(1)()**()(s)x kix k x k ki dt u k u s=⇒+=+ 惯性环节：(s)(1)exp(/)()[1exp(/)]()(s)1x k x k d x k k d u k u s τττττ=⇒+=-⋅+--+ 四、 实验方案 1、 分别离散法； 系统框图根据上面提到的分别离散法得到仿真的公式 已知系数：K1=0.93; K2=2.086; T1=73.3; T2=96.1; n1=2; n2=4; kp1=0.32; ki1=0.0018; kp2=2; ki2=0.00008;惯性环节的系数：fai1=exp(-dt/T1);faiM1=1-fai1;fai2=exp(-dt/T2);faiM2=1-fai2;PID 控制环节：up1=e*kp1;x(1)=x(1)+ki1*dt*e；up2=e1*kp2;x(2)=x(2)+ki2*dt*e1;惯性环节：x(3)=fai1*x(3)+K1*faiM1*u1;x(4)=fai1*x(4)+faiM1*x(3);x(5)=fai2*x(5)+K2*faiM2*x(4); x(6)=fai2*x(6)+faiM2*x(5);x(7)=fai2*x(7)+faiM2*x(6);x(8)=fai2*x(8)+faiM2*x(7); 2、整体离散法；将系统框图拆开系统的状态空间方程为：'1181'22124128121212112112'1312348111111'43411'254522'65622'76722'878221111111111i i i i p i p i p p p p p p x k x k Rx k x k x k k x k k R K k K k K k k K k k K x x x x x x R T T T T T T x x x T T K x x x T T x x x T T x x x T T x x x T T =-+=--+=----+=-=-=-=-=-此时可以得到此时状态方程的系数12221221211111111111222222222000000000000100011000000100000011000000110000001100i i i i p p p p p K K K K K K K K K K K K K T T T T T T T K A T T T T T T T T -⎧⎫⎪⎪--⎪⎪⎪⎪⎪⎪---⎪⎪⎪⎪⎪⎪-⎪⎪⎪⎪⎪⎪=⎨⎬-⎪⎪⎪⎪⎪⎪-⎪⎪⎪⎪⎪⎪-⎪⎪⎪⎪⎪⎪-⎪⎪⎩⎭由上面的推导可知(1)()(1)m x k x k x k +=Φ⋅+Φ⋅+求出m ΦΦ和就可以得到系统的状态空间方程22312231111(......)2!3!!111(......)B2!3!!ATk k k k m e I A T AT A T A T k T AT A T A T k --Φ==++++++Φ=+++++ 在Matlab 中仿真时为for i=1:n1*n2faiM=faiM+(dt^i)*(a^(i-1))/factorial(i); endfai=faiM*a+eye(n1*n2); faiM=faiM*b; for j=1:lpx=fai*x+faiM*r; y=c*x+d*r; y1=[y1 y]; t=[t j*dt]; end3、 欧拉法由上面已经求出系统的状态空间方程，所以这里直接引用，欧拉法的求解过程如下：在Matlab 中的仿真程序如下：for i=1:lp xk=a*x+b*r; x=x+xk*dt; y=c*x+d*r; y1=[y1 y]; t=[t dt*i]; end4、 梯形法类似于欧拉法，梯形法的推导如下111122()()(),(1)()()(1)(1)(1)()(1)()2(1)()()()()()()(1)2222x Ax Bu x k Ax k Bu k x k x k dt x k x k Ax k Bu k x k x k x k x k x k dt x k dt dt dt dtI dtA A x k A Bu k Bu k ⋅⋅⋅⋅⋅⋅⋅⋅=+⇒=+⇒+=++=+++++=+=+=++++++平均一下得在Matlab 中仿真的程序如下：for i=1:lp xk=a*x+b*r; xk1=x+dt*xk; xk2=a*xk1+b*r; E=(xk+xk2)/2; x=x+dt*E; y=c*x+d*r; y1=[y1 y]; t=[t dt*i]; end5、 龙格——库塔法 推导如下：120031141234234234234234()()11()(),2211()()22(1)(1)22()6(1)()()()()2!3!4!()()23624x Ax Bu e Ax k Bu k e Ax k Bu k e Ax k Bu k e Ax k Bu k e e e e E k x k x k dtE k dt dt dt I dtA A A A x k dt dt dt dt A A A A Bu k ⋅=+⇒=+=+++=+++=++++++=+=+=++++++++在Matlab 中的仿真程序如下：for i=1:lp e1=a*x+b*r; xk1=x+dt*e1/2; e2=a*xk1+b*r; xk2=x+dt*e2/2; e3=a*xk2+b*r; xk3=x+dt*e3/2; e4=a*xk3+b*r; E=(e1+e2+e3+e4)/6; x=x+dt*E; y=c*x+d*r; y1=[y1 y]; t=[t dt*i]; end五、 实验结论 5种方法仿真图形0500100015002000250030003500400045005000放大后的图像此时，可以看出，分别离散已经开始远离其他的线继续放大此时分别离散已经明显远离其他，并且欧拉法也开始远离其他的线最终可以看出，龙格——库塔法与整体离散法得到的仿真曲线最接近。

小结：利用不同的方法对多阶系统的状态方程进行求解，分别离散法，因为零阶保持器的缘故，所以误差比较大；欧拉法通过简单的取切线的端点作为下一步的起点，提升了精确性，但是本身也存在缺点，当步数增加时，误差在逐渐累积；详细实例见附件；梯形法是欧拉法的升级版，首先可以由欧拉法求得下一时刻的值，再代入校正得到一个更精准的值，这样，可以较欧拉法得到更精准的值；龙格库塔法是至尊版，比梯形法更精准，运用不同阶数的龙格库塔法可以得到更精准的值，他运用不同预估值的斜率求取平均值，并赋予不同的权重，提高精度；六、实验中存在的问题没有明显的问题。