控制系统数字仿真与CAD习题.pdf

控制系统计算机仿真及辅助设计控制系统计算机仿真及辅助设计 课程上机作业课程上机作业 1 1 月月 5 5 日上机作业一：日上机作业一： 2-12-1 思考题：思考题： （1 1）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式， 各有什么特点？各有什么特点？ （2 2）数学模型各种形式之间为什么要互相转换？）数学模型各种形式之间为什么要互相转换？ （3 3）控制系统建模的基本方法有哪些？他们的区别和特点是什么？）控制系统建模的基本方法有哪些？他们的区别和特点是什么？ （4 4）控制系统计算机仿真中的）控制系统计算机仿真中的“实现问题实现问题”是什么含意？是什么含意？ （5 5）数值积分法的选用应遵循哪几条原则？）数值积分法的选用应遵循哪几条原则？ （1） 微分方程是直接描述系统输入和输出量之间的制约关系， 是连续控制系统其他 数学模型表达式的基础。状态方程能够反映系统内部各状态之间的相互关系，适用于多 输入多输出系统。传递函数是零极点形式和部分分式形式的基础。零极点增益形式可用 于分析系统的稳定性和快速性。利用部分分式形式可直接分析系统的动态过程。 （2）不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。 （3）控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。机 理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分 析简化建立起来的各物理量间的关系。该方法需要对系统的内部结构和特性完全的了 解，精度高。统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系 统辨识等理论建立的系统模型。该方法建立的数学模型受数据量不充分，数据精度不一 致，数据处理方法的不完善，很难在精度上达到更高的要求。混合法是上述两种方法的 结合。 （4）“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法，将模 型方程转换为适合在计算机上运行的公式和方程， 通过计算来使之正确的反映系统各变 量动态性能，得到可靠的仿真结果。 （5）数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速 度和并保证计算结果的稳定。 2-2.2-2.用用 matlabmatlab 语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式 的模型参数，并分别写出其相应的数学模型表达式：的模型参数，并分别写出其相应的数学模型表达式： (1)G（s）= 32 432 72424 10355024 sss ssss (2)Equation Section (Next) . X= 2.25 -5 -1.25 -0.54 2.25 -4.25 -1.25 -0.252 0.25 -0.5 -1.25 -12 1.25 -1.75 -0.25 -0.75 0 X u y=0 2 0 2 X 解： （1） 程序代码如下： clear; clc; num=1 7 24 24; den=1 10 35 50 24; Z,P,K=tf2zp(num,den); A,B,C,D=zp2ss(Z,P,K); R,P1,H=residue(num,den); G1=tf(num,den) G2=zpk(Z,P,K) G3=ss(A,B,C,D) 程序的输出结果如下： Transfer function: s3 + 7 s2 + 24 s + 24 - s4 + 10 s3 + 35 s2 + 50 s + 24 Zero/pole/gain: (s+1.539) (s2 + 5.461s + 15.6) - (s+4) (s+3) (s+2) (s+1) a = x1x2x3x4 x1-3-1.41400 x21.414000 x311.088-7-3.464 x4003.4640 b = u1 x11 x20 x30 x40 c = x1x2x3x4 y111.088-1.5391.038 d = u1 y10 Continuous-time model. （2） 程序代码如下： clear; clc; A=2.25 -5 -1.25 -0.5 2.25 -4.25 -1.25 -0.25 0.25 -0.5 -1.25 -1 1.25 -1.75 -0.25 -0.75; B=4 2 2 0; C=0 2 0 2; D=0; num,den=ss2tf(A,B,C,D); Z,P,K=ss2zp(A,B,C,D); R,P1,H=residue(num,den); G1=tf(num,den) G2=zpk(Z,P,K) G3=ss(A,B,C,D) 程序输出结果如下： Transfer function: 4 s3 + 14 s2 + 22 s + 15 - s4 + 4 s3 + 6.25 s2 + 5.25 s + 2.25 Zero/pole/gain: 4 (s+1.5) (s2 + 2s + 2.5) - (s+1.5)2 (s2 + s + 1) a = x1x2x3x4 x12.25-5-1.25-0.5 x22.25-4.25-1.25-0.25 x30.25-0.5-1.25-1 x41.25-1.75-0.25-0.75 b = u1 x14 x22 x32 x40 c = x1x2x3x4 y10202 d = u1 y10 Continuous-time model. 2-3.用欧拉法求下面系统的输出响应用欧拉法求下面系统的输出响应 y(t)在在 0t1 上，上，h=0.1 时的数值。时的数值。 , (0)1yy y 要求保留要求保留 4 位小数，并将结果与真解位小数，并将结果与真解( ) t y te比较。比较。 解： 程序代码如下： clear; clc; h=0.1; y=1; disp(欧拉法求得值为：); for i=0:h:1 m=y; disp(y); y=m-m*h; end disp(真解为：); for i=0:h:1 y=exp(-i); disp(y); end 输出结果如下： 欧拉 法值 10.900 0 0.810 0 0.729 0 0.656 1 0.590 5 0.531 4 0.478 3 0.430 5 0.387 4 0.348 7 真值10.904 8 0.818 7 0.740 8 0.670 3 0.606 5 0.548 8 0.496 6 0.449 3 0.406 6 0.367 9 2-4 用二阶龙格库塔法求解用二阶龙格库塔法求解 2-3 的数值解，并于欧拉法求得的结果比较。的数值解，并于欧拉法求得的结果比较。 解： 程序代码如下： clear; clc; h=0.1; y=1; disp(二阶龙格库塔法求得值为：); for i=0:h:1 disp(y); k1=-y; k2=-(y+k1*h); y=y+(k1+k2)*h/2; end disp(欧拉法求得值为：); y=1; for i=0:h:1 m=y; disp(y); y=m-m*h; end 程序输出结果如下： 欧拉 法值 10.900 0 0.810 0 0.729 0 0.656 1 0.590 5 0.531 4 0.478 3 0.430 5 0.387 4 0.348 7 龙格 库塔 法值 10.905 0 0.819 0 0.741 2 0.670 8 0.607 1 0.549 4 0.497 2 0.450 0 0.407 2 0.368 5 2-5用四阶龙格用四阶龙格-库塔法求解题库塔法求解题 2-3 数值解，并与前两题结果相比较。数值解，并与前两题结果相比较。 解： 程序代码如下： clear; clc; h=0.1; y=1; disp(四阶龙格库塔法求得值为：); for i=0:h:1 disp(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 程序结果如下： 欧拉 法值 10.900 0 0.810 0 0.729 0 0.656 1 0.590 5 0.531 4 0.478 3 0.430 5 0.387 4 0.348 7 二阶 龙格 库塔 法值 10.905 0 0.819 0 0.741 2 0.670 8 0.607 1 0.549 4 0.497 2 0.450 0 0.407 2 0.368 5 四阶 龙格 库塔 法值 10.904 8 0.818 7 0.740 8 0.670 3 0.606 5 0.548 8 0.496 6 0.449 3 0.406 6 0.367 9 2-6已知二阶系统状态方程为已知二阶系统状态方程为 . 1 1112 10111 . 22220 2122 2 (0) ; (0) x aaxxbx u xbxxaa x 写出取计算步长为写出取计算步长为 h 时，该系统状态变量时，该系统状态变量 X= 12 ,x x的四阶龙格的四阶龙格-库塔法递推关系式库塔法递推关系式。 2-7 单位反馈系统的开环传递函数已知如下单位反馈系统的开环传递函数已知如下 2 5100 ( ) (4.6)(3.416.35) s G s s sss 用用 matlab 语句语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标 准型实现。准型实现。 解： 程序代码如下： clear; clc; num=5 100; den1=1 0; den2=1 4.6; den3=1 3.4 16.35; den4=0 0 0 5 100; den=conv(den1,den2); den=conv(den,den3); den=den+den4 Z,P,K=tf2zp(num,den); P1=roots(den) Z1=roots(num) Z P 输出结果如下： den = 1.00008.000031.990080.2100100.0000 P1 = -0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697i Z1 = -20 Z = -20 P = -0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697i 2-8 用用 matlab 语言编制单变量系统三阶龙格语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能接收状库塔法求解程序，程序入口要求能接收状 态方程各系数阵态方程各系数阵（A,B,C,D）,和输入阶跃函数和输入阶跃函数 r(t)=R*1(t);程序出口应给出输出量程序出口应给出输出量 y（t） 的动态响应数值解序列的动态响应数值解序列 01 , n yyy。 解： 函数程序代码如下： function y = hs( A,B,C,D,R,T,h ) y=0; r=R; x=0;0;0;0; N=T/h; for t=1:1:N k1=A*x+B*R; k2=A*(x+h*k1/3)+B*R; k3=A*(x+2*h*k2/3)+B*R; x=x+h*(k1+3*k3)/4; y=C*x+D*R; end 将程序保存为名为 hs 的文件，用以求解 2-1 求解出的 A,B,C,D 进行验证如下： hs(A,B,C,D,1,1,1) ans = 2.6667 2-10.用式用式（2-34）梯形法求解试验方程梯形法求解试验方程 1 yy ，分析对计算步长分析对计算步长 h 有何限制有何限制，说说明明 h 对数值稳定性的影响。对数值稳定性的影响。 1 1 月月 6 6 日上机作业二：日上机作业二： 3-1 求解下列线性方程，并进行解得验证：求解下列线性方程，并进行解得验证： (1) 7 2 1 -24 9 15 3 -27 -2 -2 11 51 1 3 2 130 x ,(2) 5 7 6 5 124 7 10 8 7 234 6 8 10 9 336 5 7 9 10 435 1 2 3 4 515 x 96 136 144 140 60 由由 A*X=B 得：得：X=AB 解： （1） 程序代码如下： clear; clc; A=7 2 1 -2 9 15 3 -2 -2 -2 11 5 1 3 2 13; B=4 7 -1 0; x=inv(A)*B 输出结果如下： x = 0.4979 0.1445 0.0629 -0.0813 （2） 程序代码如下： clear; clc; A= 5 7 6 5 1 7 10 8 7 2 6 8 10 9 3 5 7 9 10 4 1 2 3 4 5; B=24 96 34 136 36 144 35 140 15 60; x=inv(A)*B 输出结果如下： x = 1.00004.0000 1.00004.0000 1.00004.0000 1.00004.0000 1.00004.0000 3-2进行下列计算，给出不使用进行下列计算，给出不使用 for 和和 while 等循环语句的计算方法。等循环语句的计算方法。 （1） 63 0 2i i k (2)求出求出 y=x*sin(x) 在在 0x100 条件下的每个峰值条件下的每个峰值 解： （1） 程序代码如下： clear; clc; a1=1; q=2; n=64; sum=a1*(1-qn)/(1-q) 输出结果如下： sum = 1.8447e+019 （2） 先画出 0x100 条件下 y 的图形 代码如下： x=0:0.1:100; y=x.*sin(x); plot(x,y) xlabel(x); ylabel(y); 图形如下： 可以看出 sin(x)在峰值时 y 也是在峰值位置，即每个点的峰值为)31,1 , 0( 2 kk 3-33-3绘制下面的图形。绘制下面的图形。 （1 1）sin(1/t),-1t1(2)(2) 3 1 cos (7 ) t-1t1 解： （1） 代码如下： t=-1:0.01:1; y=sin(1./t); plot(t,y); xlabel(t); ylabel(y); 输出图形如下： （2） 代码如下： clear; clc; t=-1:0.01:1; y=1-(cos(7*t).3; plot(t,y); xlabel(t); ylabel(y); 输出图形如下： 3-4已知元件的实验数据如下，拟合这一数据，并尝试给出其特性方程。已知元件的实验数据如下，拟合这一数据，并尝试给出其特性方程。 X0.01001.01002.01003.01004.0100 Y2.54377.88849.624211.607111.9727 X5.01006.01007.01008.01009.0100 y13.218914.267914.613415.404515.0805 解： 程序代码如下： clear; clc; x=0.01001.01002.01003.01004.01005.01006.0100 7.01008.01009.0100; y=2.54377.88849.624211.607111.972713.218914.2679 14.613415.404515.0805 ; P=polyfit(x,y,3); xi=0:0.01:9.01; yi=polyval(P,xi); plot(x,y,o,xi,yi,*); plot(x,y,xi,yi); xlabel(x); ylabel(y); 输出结果如下： 3-5分别使用解微分方程方法分别使用解微分方程方法、控制工具箱控制工具箱、simulink 求解具有如下闭环传递函数的系统的阶跃响求解具有如下闭环传递函数的系统的阶跃响 应。应。 432 10 ( ) 8364010 s ssss 解： Simulink 中建立如下系统 阶跃响应图像如下： 输入以下代码： num=6 26 6 20; den=1 3 4 2 2; step(num,den,10); 输出波形如下 3-6已知系统的闭环传递函数已知系统的闭环传递函数 32 432 626620 ( ) 3422 sss s ssss ，试分析该系统的稳定性。，试分析该系统的稳定性。 解： 画出系统零极点图代码如下： num=6 26 6 20; den=1 8 36 40 10; pzmap(num,den) 输出结果为： 可以看出系统的几点都在左半平面中，所以该系统为稳定系统。 3-8某小功率随动系统动态结构如下图所示，已知：某小功率随动系统动态结构如下图所示，已知： 12012 0.01,0.05,1,300, 1,0.08. c TTKKK K 若系统输入分别为若系统输入分别为( )1( ),1( ) 1(1.5) srsrsr tttt，适用，适用 simulink 分析系统的输出分析系统的输出 ( )sc t 分别如何？分别如何？ （1）输入为）输入为 1（t）: 输出为： （2）输入为）输入为 t 时：时： 输出： （3）输入为输入为 l(t)-l(t-1.5)时：时： 1 1 月月 7 7 日上机作业三：日上机作业三： 4-2 设典型闭环结构控制系统如图设典型闭环结构控制系统如图 4-47 所示所示， 当阶跃输入幅值当阶跃输入幅值20R 时时， 用用 sp4_1.m 求取输出求取输出( )y t 的响应。的响应。 ( )y t ( )r t 2 432 3025 0.0160.8643.273.421 s ssss 解： 代码为： a=0.016 0.864 3.27 3.42 1; b=30 25; X0=0 0 0 0;%系统状态向量初值为零 V=2;%反馈系数 n=4; T0=0;Tf=10; h=0.01;R=20 ; sp4_1; plot(t,y) 输出结果为： 4-8 求图求图 4-49 非线性系统的输出响应非线性系统的输出响应 y(t),并与无非线性环节情况进行比较。并与无非线性环节情况进行比较。 0.5 0.1 s s 20 (2)(10)s ss 5 5 ( ) 10r t ( )e t ( )u t ( )y t 解： 程序代码为： P=0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0; WIJ=1 0 1;1 4 -1;2 1 1;3 2 1;4 3 1 ; Z=0