控制系统数字仿真习题答案.doc

控制系统数字仿真与CAD第二章习题答案2-1 思考题：（1）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，各有什么特点？（2）数学模型各种形式之间为什么要互相转换？（3）控制系统建模的基本方法有哪些？他们的区别和特点是什么？（4）控制系统计算机仿真中的“实现问题”是什么含意？（5）数值积分法的选用应遵循哪几条原则？答：（1）微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础。状态方程能够反映系统内部各状态之间的相互关系，适用于多输入多输出系统。传递函数是零极点形式和部分分式形式的基础。零极点增益形式可用于分析系统的稳定性和快速性。利用部分分式形式可直接分析系统的动态过程。（2）不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。（3）控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系。该方法需要对系统的内部结构和特性完全的了解，精度高。统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型。该方法建立的数学模型受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到更高的要求。混合法是上述两种方法的结合。（4）“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法，将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反映系统各变量动态性能，得到可靠的仿真结果。（5）数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速度和并保证计算结果的稳定。2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：(1) G（s）= (2) =y=0 2 0 2 X（1） 解：（1）状态方程模型参数:编写matlab程序如下 num=1 7 24 24; den=1 10 35 50 24; A B C D=tf2ss(num,den) 得到结果：A=,B=,C=,D=0所以模型为: =X+u,y=X (2) 零极点增益:编写程序 num=1 7 24 24; den=1 10 35 50 24; Z P K=tf2zp(num,den)得到结果Z= -2.7306 + 2.8531 , -2.7306 - 2.8531i ,-1.5388 P= -4, -3 ,-2 ,-1 K=1 (3) 部分分式形式:编写程序 num=1 7 24 24; den=1 10 35 50 24; R P H=residue(num,den) 得到结果R= 4.0000 ,-6.0000, 2.0000, 1.0000 P= -4.0000, -3.0000 , -2.0000 ,-1.0000 H= G(s)=（2）解：（1）传递函数模型参数：编写程序 A=2.25 -5 -1.25 -0.52.25 -4.25 -1.25 -0.250.25 -0.5 -1.25 -11.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) 得到结果num = 0 4.0000 14.0000 22.0000 15.0000den =1.0000 4.0000 6.2500 5.2500 2.2500 (2) 零极点增益模型参数:编写程序 A=2.25 -5 -1.25 -0.52.25 -4.25 -1.25 -0.250.25 -0.5 -1.25 -11.25 -1.75 -0.25 -0.75； B=4 2 2 0; C=0 2 0 2; D=0; Z,P,K=ss2zp(A,B,C,D) 得到结果Z =-1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000P= -0.5000 + 0.8660i -0.5000 - 0.8660i -1.5000 -1.5000 K = 4.0000 表达式 （3）部分分式形式的模型参数：编写程序 A=2.25 -5 -1.25 -0.52.25 -4.25 -1.25 -0.250.25 -0.5 -1.25 -11.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) R,P,H=residue(num,den) 得到结果R = 4.0000 -0.0000 0.0000 - 2.3094i 0.0000 + 2.3094iP = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iH =2-3.用欧拉法求下面系统的输出响应y(t)在0t1上，h=0.1时的数值。 要求保留4位小数，并将结果与真解比较。解：欧拉法（前向欧拉法,可以自启动）其几何意义：把f(t,y)在区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示（1） m文件程序为 h=0.1;disp(函数的数值解为); %显示 中间的文字%disp(y=); %同上%y=1;for t=0:h:1 m=y; disp(y); %显示y的当前值% y=m-m*h;end保存文件q2.m 在matalb命令行中键入 q2 得到结果 函数的数值解为y= 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487（2）另建一个m 文件求解在t0,1的数值 （ %是的真解%）程序为h=0.1;disp(函数的离散时刻解为);disp(y=);for t=0:h:1 y=exp(-t); disp(y);end 保存文件q3.m在matalb命令行中键入 q3 函数的离散时刻解为y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679比较欧拉方法求解与真值的差别欧拉10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679误差0-0.0048-0.00070.01180.01420.01600.01740.01830.0188-0.0192-0.0192显然误差与为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。解：我们经常用到 预报-校正法 的二阶龙-格库塔法， 此方法可以自启动，具有二阶计算精度，几何意义：把f(t,y)在区间内的曲边面积用上下底为和、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示 （1）m文件程序为 h=0.1;disp(函数的数值解为);disp(y=);y=1;for t=0:h:1 disp(y); k1=-y; k2=-(y+k1*h); y=y+(k1+k2)*h/2;end 保存文件q4.m在matlab的命令行中键入 q4 显示结果为 函数的数值解为y= 1 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685（2） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龙库10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685误差00.00020.00030.00040.00050.00060.00060.00060.00070.00060.0006明显误差为得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。2-5用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。解：四阶龙格-库塔法表达式，其截断误差为同阶无穷小，当h步距取得较小时，误差是很小的. (1)编辑m文件程序h=0.1;disp(四阶龙格-库塔方法求解函数数值解为);disp(y=);y=1;for t=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 保存文件q5.m 在matlab命令行里键入 q5得到结果 四阶龙格-库塔方法求解函数数值解为y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.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-6已知二阶系统状态方程为写出取计算步长为h时，该系统状态变量X=的四阶龙格-库塔法递推关系式。解：四阶龙格-库塔法表达式所以状态变量的递推公式可以写作： A=,B=,可以写成则递推形式2-7单位反馈系统的开环传递函数已知如下 用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。解：已知开环传递函数，求得闭环传递函数为 在matlab命令行里键入 a=1 0; b=1 4.6; c=1 3.4 16.35; d=conv(a,b)； e=conv(d,c)e = 1.0000 8.0000 31.9900 75.2100 0 f=0 0 0 5 100; g=e+fg = 1.0000 8.0000 31.9900 80.2100 100.0000%以上是计算闭环传递函数的特征多项式% p=roots(g) %计算特征多项式的根，就是闭环传递函数的极点%p =-0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697i m=5 100; z=roots(m)z = -20 %计算零点% 综上：当闭环传函形如时，可控标准型为：； 所以可控标准型是2-8用matlab语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能接收状态方程各系数阵（A,B,C,D）,和输入阶跃函数r(t)=R*1(t);程序出口应给出输出量y（t）的动态响应数值解序列。解：m文件为：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*k2/3)+B*R;x=x+h*(k1+3*k3)/4;y(t)=C*x+D*R;end在命令行里键入A= B= C= D= R= T= h=y=hs(A,B,C,D,R,T,h) 得到结果。2-9用题2-8仿真程序求解题2-7系统的闭环输出响应y(t).解：A=,B=,C=,D=0在命令行里键入 A=0 1 0 00 0 1 00 0 0 1-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)数值解为 0 8.3333e-007 5.8659e-006 1.8115e-005 3.9384e-005 7.0346e-005。 %仅取一部分%2-10.用式（2-34）梯形法求解试验方程，分析对计算步长h有何限制，说明h对数值稳定性的影响。解：编写梯形法程序为得到 稳定系统最终渐进收敛。系统稳定则 计算得。h的选取不能超出上述范围，否则系统不稳定。2-11如图2-27所示斜梁滚球系统，若要研究滚球在梁上的位置可控性，需首先建立其数学模型，已知力矩电机的输出转矩M与其电流i成正比，横梁为均匀可自平衡梁（即当电机不通电且无滚球时，横梁可处于=0的水平状态），是建立系统的数学模型，并给出简化后系统的动态结构图。解：设球的质心到杆的距离为0，该系统为特殊情况下的球棒系统。另令分别表示棒的惯量、球的质量和球的惯量。则球质心的位置和速度为其中，。因而动能的移动部分为因而动能的移动部分为 球棒系统的旋转动能为 因而，系统总的动能等于其中为常数。此系统的拉格朗日方程组为综合以上公式的系统的方程组为设系统在平衡点附近，则系统方程可化为对上式进行拉普拉斯变换并化简后可得到。参考文献：1 Hauser, S. Sestry, and P. Kokotovi