控制系统仿真及CAD试题研

1、控制系统仿真及CAD 试题（研 2012）1、4-2 某反馈控制系统的开环传递函数为G s Hs试绘制其根轨迹。Ks s4s24s20解：在 MATLAB 命令窗口中输入下列命令：ng=1.0;dg=poly(0,-4,-2+4j,-2-4j);rlocus(ng,dg)运行结果为：4-3 已知某系统传递函数为80( 1s 1)W (s)100( 1 s1)( s ) 22 0.31 s 140200200试绘制其伯德图。解：分子分母同乘100*200 得到80200( s100)W ( s)2s(2.5s100)(20.3s200)200在 Matlab 窗口中输入下列命令：k=80*200

2、;num=1 100;den=conv(2.5 100, (1/200) 2*0.3 200);w=logspace(-1,1,100);m,p=bode(k*num,den,w);subplot(2,1,1);semilogx(w,20*log10(m);grid;xlabel('Frequency(rad/sec)');ylabel('Gain(dB)');subplot(2,1,2);semilogx(w,p);grid;xlabel('Frequency(rad/s)');ylabel('Phase(deg)');可绘制该

3、系统的伯德图如下所示。4-4 设控制系统具有如下的开环传递函数G s HsKs s 1s5试求取当 K=10 时的相角裕度和幅值裕度，并画出其伯德图。解：在 MATLAB 命令窗口中输入下列命令：k=10;num=1;den=poly(0,-1,-5);m,p,w=bode(k*num,den);subplot(2,1,1);semilogx(w,20*log10(m);grid on;ylabel('Gain(dB)');subplot(2,1,2);semilogx(w,p);grid on;xlabel('Frequency(rad/sec)');ylab

4、el('Phase(deg)');gm,pm,wcg,wcp=margin(m,p,w)这里 gm,wcg 为幅值裕度值与相应的频率pm,wcp 为相角裕度值与相应的频率，运行结果为： gm =3.0000， pm = 25.4489， wcg =2.2361， wcp = 1.2241。因此，系统的幅值裕度和相角裕度分别为3dB 和 25.4489°。系统的伯德图如下所示。4-13 对于高阶系统的设计问题， 往往要进行降阶近似处理， 并要验证近似效果。已知某高阶系统模型为194480422964 s511812 s2278376 s382402 s413285 s5

5、 10861s635s7G (s)28880 s37492 s227470 s3 11870 s43017 s5437s633s7s89600经简化处理后，模型等效为284.9835sG (s)12.3110.114ss2试比较两个模型在单位阶跃信号作用下的响应情况，并分析近似效果。解：在 Matlab 命令窗口中输入下列命令：num1=35,10861,13285,82402,278376,511812,422964,194480;den1=1,33,437,3017,11870,27470,37492,28880,9600;sys1=tf(num1,den1);num2=35,284.98

6、;den2=1,10.114,12.31;sys2=tf(num2,den2);step(sys1,'-',sys2,':r',15)legend('原模型 ', '简化后模型 ') ;grid on程序运行结果如下：从曲线中可以看出， 降阶后系统响应无超调， 调整时间缩短，系统响应变快，但是简化前后终值有差异。2-5 用四阶龙格 -库塔法求解题2-3 数值解，并与前两题结果相比较。yk 1ykh (k12k2 2k3k4 )6k1f (tk , yk )解：四阶龙格 -库塔法表达式hh，其截断误差为5k2f (tk2 , yk2

7、 k1 )hk3f (tkh , ykh k2 )22k4f (tkh, ykhk3 )同阶无穷小，当h 步距取得较小时，误差是很小的.(1) 程序如下： h=0.1;disp('四阶龙格 -库塔方法求解函数数值解为 '); disp('y=');y=1;for t=0:h:1disp(y);k1=-y;k2=-(y+k1*h/2);k3=-(y+k2*h/2);k4=-(y+k3*h);y=y+(k1+2*k2+2*k3+k4)*h/6;end得到结果四阶龙格 -库塔方法求解函数数值解为y=10.90480.81870.74080.67030.60650.54

8、880.49660.44930.40660.3679(2) 比较这几种方法：对于四阶龙格 -库塔方法真10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679值龙10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679库误00000000000差显然四阶龙格 -库塔法求解精度很高，基本接近真值。三种方法比较可以得到精度（四阶） 精度（二阶）精度（欧拉）2-9 用题 2-8 仿真程序求解题2-7 系统的闭环输出响应y(t).解：题 2-7 单位反馈系统的开环传递函数

9、已知如下G (s)5s1004.6)( s23.4s16.35)s( s已知开环传递函数，求得闭环传递函数为G(s)5s1004.6)( s23.4s16.35)5s 100s( s在 matlab 命令行里键入 >> a=1 0;>> b=1 4.6;>> c=1 3.4 16.35;>> d=conv(a,b) ;>> e=conv(d,c)e =1.00008.000031.990075.21000>> f=0 0 0 5 100;>> g=e+fg = 1.00008.000031.990080.210

10、0100.0000% 以上是计算闭环传递函数的特征多项式%>> m=5 100;>> z=roots(m)z = -20% 计算零点 %0100000100那么 A=00,B=,C= 100 5 0 0 ,D=001010080.2131.9981m 文件为：function y= hs(A,B,C,D,R,T,h)%T 为观测时间 ,h 为计算步长， R 为输入信号幅值 %disp('数值解为 ');y=0;r=R;x=0;0;0;0;N=T/h;for t=1:N;k1=A*x+B*R;k2=A*(x+h*k1/3)+B*R;k3=A*(x+2*h*

11、k2/3)+B*R;x=x+h*(k1+3*k3)/4;y(t)=C*x+D*R;end在命令行里键入 >> A=0 1 0 000100001-100 -80.21 -31.99 -8;>> B=0 0 0 1'>> C=-100 5 0 0;>> D=0;>> T=1;>> R=1;>> h=0.01;>> y= hs (A,B,C,D,R,T,h)数值解为08.3333e-0075.8659e-0061.8115e-0053.9384e-0057.0346e-005。 %仅取一部分 %

12、3-5 图 3-61 中，若各环节的传递函数已知为：G (s)1; G( s)10.17s ; G (s)1;110.01s20.085s310.01sG4 (s)10.15s ; G5 (s)170; G6( s)10.21 ;0.051s0.0067s0.15sG7 (s)1300.1; G9 ( s)0.0044s;G8 (s)11;0.01s0.01s但 G10（s）=0.212 ；重新列写联接矩阵 W ,W0 , 和非零元素矩阵 WIJ ，将程序 sp4_2.m完善后，应用 sp4_2.m求输出 y7 的响应曲线。y0y7G1(s)G2(s)G3(s)G4(s)G5(s)G6(s)G

13、7(s)-G8(s)G10(s)G9(s)解：根据图中 ui , yi 拓扑连结关系，可写出每个环节输入 ui受哪些环节输出 yi 的影响，现列如入下 :u1y0u2y1y9u3y2u4y3y8u5y4u6y50.212 y7u7y6u8y 6u9y7把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来，U WY W0Y00000000001100000001001000000000010000100（2）W 000100000W000000100.212000000001000000000100000000001000101211291321431WIJ481541651670.2127

14、61861971程序为：P=1,0.01,1,0;0,0.085,1,0.17;1,0.01,1,0;0,0.051,1,0.15;1,0.0067,70,0;1,0.15,0.21,0;0,1,130,0;1,0.01,0.1,0;1,0.01,0.0044,0;WIJ=1,0,1;2,1,1;2,9,-1;3,2,1;4,3,1;4,8,-1;5,4,1;6,5,1;6,7,-0.212;7,6,1;8,6,1;9,7,1;n=9;Y0=10;Yt0=0,0,0,0,0,0,0,0,0;h=0.001;L1=10;T0=0;Tf=5;nout=7;%. 形成闭环各系数阵 .A=diag(P

15、(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:mif (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3);else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;%.数值积分求解 Y=Yt0'y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);R=Y0;for i=1

16、:N;for j=1:L1;k1=Ab*Y+b1*R;k2=Ab*(Y+h*k1/2)+b1*R;k3=Ab*(Y +h*k2/2)+b1*R;k4=Ab*(Y+k3*h)+b1*R;Y=Y+(k1+2*k2+2*k3+k4)*h/6;end;y=y,Y(nout);t=t,t(i)+h*L1;end;t',y'plot(t,y); xlabel('Time t/s')ylabel('幅度 ')3-7 用离散相似法仿真程序sp3_4.m 重求题 3-5 输出 y7 的数据与曲线，并与四阶龙格库塔法比较精度。解：离散相似法 ：P=1,0.01,1,

17、0;0,0.085,1,0.17;1,0.01,1,0;0,0.051,1,0.15;1,0.0067,70,0; 1,0.15,0.21,0;0,1,130,0;1,0.01,0.1,0;1,0.01,0.0044,0;WIJ=1,0,1;2,1,1;2,9,-1;3,2,1;4,3,1;4,8,-1;5,4,1;6,5,1;6,7,-0.212;7,6,1;8,6,1;9,7,1; n=9;Y0=10;Yt0=0,0,0,0,0,0,0,0,0;h=0.001;L1=10;T0=0;Tf=5; nout=7;A=diag(P(:,1); B=diag(P(:,2); C=diag(P(:,

18、3); D=diag(P(:,4); m=length(WIJ(:,1); W0=zeros(n,1);W=zeros(n,n);for k=1:mif (WIJ(k,2)=0); W0(WIJ(k,1)=WIJ(k,3);else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);end;end;for i=1:nif(A(i,i)=0);FI(i,i)=1;FIM(i,i)=h*C(i,i)/B(i,i);FIJ(i,i)=h*h*C(i,i)/B(i,i)/2;FIC(i,i)=1;FID(i,i)=0;if (D(i,i)=0);FID(i,i)=D(i,i)/B(i,i);e

19、lseendelseFI(i,i)=exp(-h*A(i,i)/B(i,i);FIM(i,i)=(1-FI(i,i)*C(i,i)/A(i,i);FIJ(i,i)=h*C(i,i)/A(i,i)-FIM(i,i)*B(i,i)/A(i,i);FIC(i,i)=1;FID(i,i)=0;if (D(i,i)=0);FIC(i,i)=C(i,i)/D(i,i)-A(i,i)/B(i,i);FID(i,i)=D(i,i)/B(i,i);elseendendendY=zeros(n,1);X=Y; y=0; Uk=zeros(n,1); Ub=Uk; t=T0:h*L1:Tf; N=length(t)

20、; for k=1:N-1for i=1:L1Ub=Uk;Uk=W*Y+W0*Y0;Udot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'*X+FIM'*Uk+FIJ'*Udot;Y=FIC'*X+FID'*Uf;endy=y,Y(nout);endplot(t,y,'*r');xlabel('Time t/s')ylabel('幅度 ')分析：用离散相似法与四节龙格库塔比较来看，两者精度很接近，离散相似法更容易推广解决非线性系统仿真问题。46 某过程控制系统如图5-71 所示，试设计 PID 调节

21、器参数，使该系统动态性能达到最佳。RePID0.1s 110sYs25s6e解：本题选用 ITSE 准则为目标函数。首先，编写 M 函数文件，存为optm.m。function ss=optm(x)global kp;global ki;global kd;kp=x(1);ki=x(2);kd=x(3);tt,xx,yy=sim('opt',40,);ss=yy(end);其次，输入恰当的模块参数，使用SIMULINK建立仿真模型文件，存为opt.mdl。11ClocksOut1ProductIntegratoru2FcnScopePID0.1s+1s2 +5s+6StepPI

22、D ControllerTransfer Fcn1TransportScope1Delay在 MATLAB 命令窗口中输入下列指令 : global kp;global ki; global kd; result=fminsearch('optm',2,1,2)运行结 果为： result = 2.6914 0.4610 2.1267,即 kp=2.6914， ki=0.4610；kd=2.1267。用 SIMULINK仿真验证控制器的控制效果，得到阶跃响应曲线如下图所示(蓝色实线所示 )，与初始值kp2，ki=1，kd=2 时的阶跃响应比较（图中红色虚线所示），动态性能得到了

23、很好的改善。仿真结果如下图所示。2优化后参数1.8原 参 数1.61.41.210.80.60.40.200100200300400500600700t/s3、对现代仿真技术、控制系统仿真的认识答：仿真技术是以相似原理、信息技术、系统技术及其应用领域有关的专业技术为基础，以计算机和各种物理效应设备为工具，利用系统模型对实际的或设想的系统进行试验研究的一门综合性技术。仿真技术的应用已不仅仅限于产品或系统生产集成后的性能测试试验，仿真技术已扩大为可应用于产品型号研制的全过程，包括方案论证、战术技术指标论证、设计分析、生产制造、试验、维护、训练等各个阶段。仿真技术不仅仅应用于简单的单个系统，也应用于

24、由多个系统综合构成的复杂系统。CAD技术 , 是将计算机高速而精确的计算能力，大容量存储和处理数据的能力与设计者的综合分析，逻辑判断以及创造性思维结合起来，用以加快设计进程，缩短设计周期，提高设计质量的技术。控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为主要内容的现代控制理论。此外，自适应控制，自校正控制以及最优控制等现代控制测策略都可利用CAD技术实现有效的分析与设计。选 3-8 求图 3-62 非线性系统的输出响应y(t),并与无非线性环节情况进行比较。r (t ) 10e(t )s0.55u(t )20y(t)s0.15s( s2)(s10)解：（ 1）不考虑非线性环

25、节影响时，求解过程如下：1）先将环节编号标入图中。2) 在 MATLAB命令窗口下，按编号依次将环节参数输入 P 阵；>> P=0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0;3) 按各环节相对位置和联接关系，有联接矩阵如下：00011101W1000 ，W00，所以非零元素矩阵141WI J2 1 10 1 00000100321431>> WIJ=1 0 1;1 4 -1;2 1 1;3 2 1;4 3 1 ;4）由于不考虑非线性影响，则非线性标志向量和参数向量均应赋零值；>> Z=0 0 0 0;S=0 0 0 0;5）输入

26、运行参数： 开环截至频率c 约为 1，故计算步长 h 取经验公式值，即10.02，取 h=0.01；每 0.25秒输出一点。故取 L1 =25。h50 c>>h=0.01;>>L1=25;>>n=4;>>T0=0>>Tf=20;>>nout=4;>>Y0=10;>>sp3_4;>> plot(t,y,'r')>> hold on运行结果如图中红色实线所示。(2) 考虑非线性环节 N 影响时，只需将非线性标志向量 Z 和参数向量 S 的相应分量正确输入即可。在

27、MATLAB命令窗口中输入下列语句：>> Z=4 0 0 0;S=5 0 0 0;%第一个线性环节后有饱和非线性，参数值为5。>> sp3_4;>> plot(t,y,'-')运行结果如图中蓝色虚线所示。从图中可以清楚的地看出，饱和非线性环节对线性系统输出响应的影响。附： sp3_4 函数为 :A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:mif (WIJ(k,2)=0); W0(WIJ(k,1)=WIJ(k,3

28、);else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);end;end;for i=1:nif(A(i)=0);FI(i)=1;FIM(i)=h*C(i)/B(i);FIJ(i)=h*h*(C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)=0);FID(i)=D(i)/B(i);elseendelseFI(i)=exp(-h*A(i)/B(i);FIM(i)=(1-FI(i)*C(i)/A(i);FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);FIC(i)=1;FID(i)=0;if(D(i)=0);FIC(i)=C(i)/D(

29、i)-A(i)/B(i);FID(i)=D(i)/B(i);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);for k=1:N-1for i=1:L1Ub=Uk;Uk=W*Y+W0*Y0;for i=1:nif(Z(i)=0)if (Z(i)=1)Uk(i)=satu(Uk(i),S(i);endif(Z(i)=2)Uk(i)=dead(Uk(i),S(i);endif(Z(i)=3)Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i);ende

30、ndendUdot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;Yb=Y;Y=FIC'.*X+FID'.*Uf;for i=1:nif(Z(i)=0)if (Z(i)=4)Y(i)=satu(Y(i),S(i);endif(Z(i)=5)Y(i)=dead(Y(i),S(i);endif(Z(i)=6)Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i);endendendendy=y,Y(nout);end附：饱和非线性函数satu.m为：function Uc=

31、satu(Ur,S1)if(abs(Ur)>=S1)if(Ur>0)Uc=S1;else Uc=-S1;endelse Uc=Ur;end%输入数据 %系统中不能出现纯比例、纯微分环节，含有微分项系数的环节不直接与外加参考输入联接P=input('请按各环节输入参数矩阵P：');Wij=input('请输入非零元素连接阵Wij:');n=input('请输入环节个数（系统阶次）n=');Y0=input('请输入阶跃输入幅值Y0=');Yt0=input('请输入各环节初值Yt0=');h=input(

32、'请输入计算步长h=');L1=input(' 请输入打印间隔点数（每隔点 l1 输出一次）：'); T0=input(' 请输入起始时间 T0=');Tf=input('请输入终止时间Tf=');nout=input('请输入输出环节编号： ');%形成闭环各系数阵 %A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(Wij(:,1);%求非零元素的个数W0=zeros(n,1);W=zeros(n,n);%建立初始 W（n*n 型方

33、阵）、W0阵（ n 维列向量）for k=1:mif (Wij(k,2)=0);W0(Wij(k,1)=Wij(k,3);else W(Wij(k,1),Wij(k,2)=Wij(k,3);endendQ=B-D*W;Qn=inv(Q);%求 Q和 Q的逆R=C*W-A;V1=C*W0;%求 R和 V1阵Ab=Qn*R;b1=Qn*V1;%形成闭环系数阵%数值积分求解 %Y=Yt0'y=Y(nout);t=T0;%置初值，做好求解准备N=round(Tf-T0)/(h*L1);for i=1:N;%每循环一次，输出一点数据for j=1:L1;%每输出点之间计算L1 次K1=Ab*Y+

34、b1*Y0;K2=Ab*(Y+h*K1/2)+b1*Y0;K3=Ab*(Y+h*K2/2)+b1*Y0;K4=Ab*(Y+h*K3)+b1*Y0;Y=Y+h*(K1+2*K2+2*K3+K4)/6;endy=y,Y(nout);t=t,t(i)+h*L1;endshuchu=t',y'plot(t,y)执行后在命令窗口输入：请按各环节输入参数矩阵 P：1 0.01 1 0;0 0.085 1 0.17;1 0.01 1 0;0 0.051 1 0.15;1 0.0067 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.004

35、40请输入非零元素连接阵 Wij:1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1;4 8 -1;5 4 1;6 5 1;6 7 0.212;7 6 1;8 6 1;9 7 1请输入环节个数（系统阶次）n=9请输入阶跃输入幅值Y0=1请输入各环节初值Yt0=0 0 0 0 0 0 0 0 0请输入计算步长 h=0.1请输入打印间隔点数（每隔点l1 输出一次）：3请输入起始时间 T0=0请输入终止时间 Tf=10请输入输出环节编号： 7执行结果为：2-3.用欧拉法求下面系统的输出响应y(t)在 0 t 1 上， h=0.1 时的数值。y 'y, y(0) 1要求保留 4 位小

36、数，并将结果与真解y(t ) e t 比较。yk 1yk h * f (tk , yk )解：欧拉法 y'f (tk , yk )（前向欧拉法 ,可以自启动）其几何意义：把 f(t,y)y(t 0 )y0在 tk , yk 区间内的曲边面积用 矩形面积 近似代替。利用 matlab 提供的 m 文件编程，得到算法公式。如下所示（ 1） m 文件程序为 h=0.1;disp('函数的数值解为 ');%显示 中间的文字 %disp('y=');%同上 %y=1;for t=0:h:1m=y;disp(y);%显示 y 的当前值 %y=m-m*h;end保存文件 q2.m在 matalb 命令行中键入 >> q2得到结果 函数的数值解为y=10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487（ 2 ）另 建一 个 m 文件 求解 ye t 在t0,1 的数 值（% ye t 是y