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

第一章 习题 1-1什么是仿真？它所遵循的基本原则是什么？ 答：仿真是建立在控制理论，相似理论，信息处理技术和计算技术等理论基础之上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假想的系统进行试验，并借助专家经验知识，统计数据和信息资料对试验结果进行分析和研究，进而做出决策的一门综合性的试验性科学。 它所遵循的基本原则是相似原理。 1-2在系统分析与设计中仿真法与解析法有何区别？各有什么特点？ 答：解析法就是运用已掌握的理论知识对控制系统进行理论上的分析，计算。它是一种纯物理意义上的实验分析方法，在对系统的认识过程中具有普遍意义。由于受到理论的不完善性以及对事物认识的不全面性等因素的影响，其应用往往有很大局限性。 仿真法基于相似原理，是在模型上所进行的系统性能分析与研究的实验方法。 1-3数字仿真包括那几个要素？其关系如何？ 答: 通常情况下，数字仿真实验包括三个基本要素，即实际系统，数学模型与计算机。由图可见，将实际系统抽象为数学模型，称之为一次模型化，它还涉及到系统辨识技术问题，统称为建模问题；将数学模型转化为可在计算机上运行的仿真模型，称之为二次模型化，这涉及到仿真技术问题，统称为仿真实验。 1-4为什么说模拟仿真较数字仿真精度低？其优点如何？。 答：由于受到电路元件精度的制约和容易受到外界的干扰，模拟仿真较数字仿真精度低 但模拟仿真具有如下优点： （1） 描述连续的物理系统的动态过程比较自然和逼真。 （2） 仿真速度极快，失真小，结果可信度高。 （3） 能快速求解微分方程。模拟计算机运行时各运算器是并行工作的，模拟机的解题速度与原系统的复杂程度无关。 （4） 可以灵活设置仿真试验的时间标尺，既可以进行实时仿真，也可以进行非实时仿真。 （5） 易于和实物相连。 1-5什么是CAD技术？控制系统CAD可解决那些问题？ 答：CAD技术，即计算机辅助设计（Computer Aided Design）,是将计算机高速而精确的计算能力，大容量存储和数据的能力与设计者的综合分析，逻辑判断以及创造性思维结合起来，用以快速设计进程，缩短设计周期，提高设计质量的技术。 控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为主要内容的现代控制理论。此外，自适应控制，自校正控制以及最优控制等现代控制测略都可利用CAD技术实现有效的分析与设计。 1-6什么是虚拟现实技术？它与仿真技术的关系如何？ 答：虚拟现实技术是一种综合了计算机图形技术，多媒体技术，传感器技术，显示技术以及仿真技术等多种学科而发展起来的高新技术。 1-7什么是离散系统？什么是离散事件系统？如何用数学的方法描述它们？ 答：本所讲的“离散系统”指的是离散时间系统，即系统中状态变量的变化仅发生在一组离散时刻上的系统。它一般采用差分方程，离散状态方程和脉冲传递函数来描述。 离散事件系统是系统中状态变量的改变是由离散时刻上所发生的事件所驱动的系统。这种系统的输入输出是随机发生的，一般采用概率模型来描述。 1-8如图1-16所示某卫星姿态控制仿真实验系统，试说明： （1） 若按模型分类，该系统属于那一类仿真系统？ （2） 图中“混合计算机”部分在系统中起什么作用？ （3） 与数字仿真相比该系统有什么优缺点？ 答：（1）按模型分类，该系统属于物理仿真系统。 （2）混合计算机集中了模拟仿真和数字仿真的优点，它既可以与实物连接进行实时仿真，计算一些复杂函数，又可以对控制系统进行反复迭代计算。其数字部分用来模拟系统中的控制器，而模拟部分用于模拟控制对象。 （4） 与数字仿真相比，物理仿真总是有实物介入，效果逼真，精度高，具有实时性与在线性的特点，但其构成复杂，造价较高，耗时过长，通用性不强。 第二章习题 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.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) 得到结果 num =  0    4.0000  14.0000  22.0000  15.0000 den =1.0000    4.0000    6.2500    5.2500    2.2500 (2) 零极点增益模型参数:编写程序>> 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]; >> [Z,P,K]=ss2zp(A,B,C,D) 得到结果Z =-1.0000 + 1.2247i -1.0000 - 1.2247i  -1.5000 P= -0.5000 + 0.8660i  -0.5000 - 0.8660i  -1.5000  -1.5000          K = 4.0000 表达式 （3）部分分式形式的模型参数：编写程序>> 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) >> [R,P,H]=residue(num,den) 得到结果R =  4.0000  -0.0000  0.0000 - 2.3094i  0.0000 + 2.3094i P =  -1.5000  -1.5000  -0.5000 + 0.8660i  -0.5000 - 0.8660i H =[] 2-3.用欧拉法求下面系统的输出响应y(t)在0≤t≤1上，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 文件求解 在t [0,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 比较欧拉方法求解与真值的差别 欧拉 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487 真值 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 误差 0 -0.0048 -0.0007 –0.0118 –0.0142 –0.0160 –0.0174 –0.0183 –0.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） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值） 真值 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 龙库 1 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685 误差 0 0.0002 0.0003 0.0004 0.0005 0.0006 0.0006 0.0006 0.0007 0.0006 0.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)比较这几种方法： 对于四阶龙格-库塔方法 真值 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 龙库 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 误差 0 0 0 0 0 0 0 0 0 0 0                         显然四阶龙格-库塔法求解精度很高，基本接近真值。三种方法比较可以得到 精度（四阶 ） 〉精度（二阶） 〉精度（欧拉） 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+f g = 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 0 0 0 1 0 0 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. Kokotovic. “Nonlinear control via approximate input-output linearization”. IEEE Trans. on Automatic Control, vol.37:pp.392-398, 1992. [2] R. Sepulchre. “Slow peaking and low-gain designs for global stabilization of nonlinear systems”. submitted for IEEE TAC 1999. [3] R. Sepulchre, M. Jankovic, and P. Kokotovic Constructive Nonlinear Control. Springer-Verlag, 1997. [4] R. Teel. “Using Saturation to stabilize a class of single-input partially linear composite systems”. IFAC NOLCOS'92 Symposium, pages 369-374, June 1992. 2-12如图2-28所示双水箱系统中， 为流入水箱1的液体流量， 为流出水箱2的液体流量，试依据液容与液阻的概念，建立 的系统动态结构图。 解：根据液容和液阻的概念，可分别列出两个水箱的数学模型 对上式进行在零初始条件下进行拉普拉斯变换得 化简后可得 第三章 习题 4-2设典型闭环结构控制系统如图4-47所示，当阶跃输入幅值 时，用sp4_1.m求取输出 的响应。 解：用sp4_1.m求解过程如下： 在MATLAB语言环境下，输入以下命令语句 >> 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 ;                            %仿真步长h=0.01，阶跃输入幅值   >> sp4_1                                  %调用sp4_1.m函数 >> plot(t,y) 运行结果为： 附：sp4_1.m函数为 b=b/a(1);a=a/a(1);A=a(2:n+1); A=[rot90(rot90(eye(n-1,n)));-fliplr(A)]; B=[zeros(1,n-1),1]'; m1=length(b); C=[fliplr(b),zeros(1,n-m1)]; Ab=A-B*C*V; X=X0'; y=0;t=T0; N=round((Tf-T0)/h); for i=1:N K1=Ab*X+B*R; K2=Ab*(X+h*K1/2)+B*R; K3=Ab*(X+h*K2/2)+B*R; K4=Ab*(X+h*K3)+B*R; X=X+h*(K1+2*K2+2*K3+K4)/6; y=[y,C*X]; t=[t,t(i)+h]; end 4-4系统结构图如图4-48，写出该系统的联结矩阵 和 ，并写出联结矩阵非零元素阵 。 解：根据图4-48中 , 拓扑连结关系，可写出每个环节输入 受哪些环节输出 的影响， 现列如入下: 把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来， 即 = ， = ， 4-6若系统为图4-5b 双输入-双输出结构，试写出该系统的联接矩阵 ， ，说明应注意什么？ 解：根据图4-5b中 , 拓扑连结关系，可列写如下关系式: 转换成矩阵形式为 所以联接矩阵 = ， = 此时应注意输入联接矩阵 变为 型。 4-8求图4-49非线性系统的输出响应y(t),并与无非线性环节情况进行比较。 解：（1）不考虑非线性环节影响时，求解过程如下： 1） 先将环节编号标入图中。 2) 在MATLAB命令窗口下，按编号依次将环节参数输入P阵； >> P=[0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0]; 3) 按各环节相对位置和联接关系，有联接矩阵如下： ，  ，所以非零元素矩阵 >> 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）输入运行参数：开环截至频率 约为1，故计算步长h取经验公式值，即 ，取h=0.01；每0.25秒输出一点。故取 =25。 >>h=0.01; >>L1=25; >>n=4; >>T0=0 >>Tf=20; >>nout=4; >>Y0=10; >>sp4_4; >> plot(t,y,'r') >> hold on 运行结果如图中红色实线所示。 (2)考虑非线性环节N影响时，只需将非线性标志向量Z和参数向量S的相应分量正确输入即可。 在MATLAB命令窗口中输入下列语句： >> Z=[4 0 0 0];S=[5 0 0 0];          %第一个线性环节后有饱和非线性，参数值为5。 >> sp4_4; >> plot(t,y,'--') 运行结果如图中蓝色虚线所示。 从图中可以清楚的地看出，饱和非线性环节对线性系统输出响应的影响。 附：sp4_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:m if  (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:n if(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); else end else FI(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(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end end end Y=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-1 for i=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; for i=1:n if(Z(i)~=0) if (Z(i)==1) Uk(i)=satu(Uk(i),S(i)); end if(Z(i)==2) Uk(i)=dead(Uk(i),S(i));  end if(Z(i)==3) [Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));  end end end Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; for i=1:n if(Z(i)~=0) if (Z(i)==4) Y(i)=satu(Y(i),S(i)); end if(Z(i)==5) Y(i)=dead(Y(i),S(i));  end if(Z(i)==6) [Y(i),Ubb(i)]=backlash(Ubb(i),Y(i),Yb(i),S(i));  end end end end y=[y,Y(nout)]; end 附：饱和非线性函数satu.m 为： function Uc=satu(Ur,S1) if(abs(Ur)>=S1) if(Ur>0) Uc=S1; else Uc=-S1; end else Uc=Ur; end 4-10采样控制系统如图4-50所示，编写程序实现对该系统的仿真分析。（提示：连续部分按环节离散化考虑）图中， 为典型数字PID控制器； =0.65为比例系数； =0.7为积分时间常数； =0.2为微分时间常数； 为具有纯滞后特性的典型二阶控制对象； ； ； 。 解：在控制对象前引入零阶保持器，将连续环节部分按环节离散化： = =         设 ，为简化运算及编程，取 为T的整数倍 = 对上式进行Z逆变换，得到 由此可编写仿真程序。 在MATLAB命令窗口中输入下列语句： >> KP=0.65;TI=0.7;TD=0.2; >> T1=0.3;a=1/T1;T3=0.4; >> T=0.1;h=0.001;Tf=10; >>hh 编写M脚本文件，存为hh.m。 %离散化后各参数为： A=1-a*h*exp(-a*h)-exp(-a*h);B=exp(-2*a*h)-exp(-a*h)+a*h*exp(-a*h);C=2*exp(-a*h);D=exp(-2*a*h); P=KP*(1+T/TI+TD/T);H=KP*(1+2*TD/T);M=KP*TD/T; %系统初始值为： E=zeros(1,3);U=zeros(1,2+T3/T+1);Y=zeros(1,2+T3/h+1);R=1; yk=0;yt=0;t=0; %仿真迭代运算： for K1=1:Tf/T ek=R-Y(1); E=[ek,E(1:2)]; uk=P*E(1)-H*E(2)+M*E(3)+U(1); U=[uk,U(1:(2+T3/T))]; for K2=1:T/h yk=A*U(T3/T+1+1)+B*U(T3/T+2+1)+C*Y(T3/h+1)-D*Y(T3/h+2); Y=[yk,Y(1:(2+T3/h))]; end yt=[yt,yk]; t=[t,K1*T]; end %输出波形： plot(t,yt) 运行结果为： 此题可以用SIMULINK仿真进行验证： 建立SIMULINK仿真模型： 运行结果为： 第四章习题 5-1 设控制系统的开环传递函数为 试画出该系统的根轨迹。 解: 在Matlab窗口中输入下列命令： num=[1 1]; a=[1 0]; b=[1 -1]; c=[1 4 16]; d=conv(a,b); den=conv(d,c); rlocus(num,den) grid on 可得到系统的根轨迹如下图所示： 5-2 某反馈控制系统的开环传递函数为 试绘制其根轨迹。 解：在MATLAB命令窗口中输入下列命令： num=1; den=conv(conv([1,0],[1,4]),[1,4,20]); rlocus(num,den) grid on 运行结果为： 5-3.已知某系统传递函数为 试绘制其伯德图。 解：分子分母同乘100*200得到 在Matlab窗口中输入下列命令： k=80*200; num=[1 100]; a=[2.5 100]; b=[(1/200) 2*0.3 200]; den=conv(a,b); 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)'); 可绘制该系统的伯德图如下所示。 5-4设控制系统具有如下的开环传递函数   试求取当K=10时的相角裕度和幅值裕度，并画出其伯德图。 解：在MATLAB命令窗口中输入下列命令： k=10; num=1; den=poly([0,-1,-5]); [m,p,w]=bode(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)'); ylabel('Phase(deg)'); [gm,pm,wcg,wcp]=margin(m,p,w) 这里gm,wcg为幅值裕度值与相应的频率 pm,wcp为相角裕度值与相应的频率，运行结果为：gm =30.0000，pm =76.5911，wcg =2.2361，wcp =0.1960。因此，系统的幅值裕度和相角裕度分别为30dB和 。系统的伯德图如下所示。 5-5已知某单位反馈系统开环传递函数为 若性能指标要求如下： ，试确定校正装置。 解：由速度误差常数的要求可求得K=200,则可绘制 的波德图。在Matlab窗口中输入下列命令： k=200*10*50*100*200; num=1; den=poly([0 -10 -50 -100 -200]); w=logspace(-1,2,200); [m,p]=bode(k*num,den,w); subplot(2,1,1); semilogx(w,20*log10(m)); grid; ylabel('Gain(dB)'); subplot(2,1,2); semilogx(w,p); grid; xlabel('Frequency(rad/s)'); ylabel('Phase(deg)'); 从图中可以看出，相角裕度为-53度, 剪切频率为37.7rad/s，这显然不能满足性能指标要求。由于补偿角为 ，单独使用超前校正装置或是滞后校正装置都不能很好的达到设计要求，因此这里考虑使用串联相位滞后超前校正。设校正装置的传递函数为 ， ， 。 按照给定的要求选定系统剪切频率为 ，为使系统稳定，应保证开环对数幅频特性以-20dB/dec的斜率穿越频率轴，即中频段的斜率为-20dB/dec。为使校正后系统的开环增益不低于200，校正前后系统的低频段特性应保持一致。中频段与低频段之间用斜率为-40dB/dec的直线连接，连接线与中频段相交的交接频率 与 不宜离的太近，取 ；同时取 ，则 , 。 最后来确定 的取值。由于校正后的剪切频率为 ，则对应 处的增益为 ;而未校正系统在 时的增益为 。两者相减就得到串联校正装置在 区间的增益，即 ，解得 。 因此串联相位滞后超前校正装置的传递函数为 校正后系统的开环传递函数为 可通过Matlab仿真验证设计结果是否满足性能指标要求。 在Matlab命令窗口中输入下列命令，绘制校正后系统的伯德图以及闭环系统的阶跃响应： z=[-1/0.6667]; p=[0,-1/8.8671,-1/0.02,-1/0.01,-1/0.0075,-1/0.005]; k=200*0.6667/(8.8671*0.02*0.01*0.0075*0.005); sys=zpk(z,p,k)                        %求出系统的开环传递函数 [gm,pm,wcg,wcp]=margin(sys)            %求出系统的幅值裕度、相位裕度及其所对应的频率 sys_closeloop=feedback(sys,1)            %求出系统的闭环传递函数 bode(sys)                            %绘制系统的伯德图 grid on figure(2) step(sys_closeloop,5)                    %绘制闭环系统的阶跃响应 grid on 运行程序后得到校正后系统的剪切频率 ， ，这些都满足期望的性能指标，同时系统伯德图如下所示。 校正后闭环系统的阶跃响应如下图所示，可见系统具有较好的动态性能与稳态性能。 5—6某过程控制系统如图5-71所示，试设计PID调节器参数，使该系统动态性能达到最佳。 解：本题选用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。 在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仿真验证控制器的控制效果，得到阶跃响应曲线如下图所示(蓝色实线所示)，与初始值 kp＝2，ki=1，kd=2时的阶跃响应比较（图中红色虚线所示），动态性能得到了很好的改善。仿真结果如下图所示。 5-7 试采用smith预估控制方法对题5-6所述系统进行重新设计，并用仿真的方法分析滞后参数变化对系统动态性能的影响。 解：设 为过程控制通道传递函数，其中 为过程中不包含纯滞后部分的传递函数， 为控制器的传递函数。在本题中， ， 。不使用Smith预估控制方法时，系统闭环传递函数为 由于在特征方程中引入了 项，使得闭环系统的品质大大恶化。Smith预估补偿控制方法的实质是通过预估补偿装置的引入，将实际工业过程中不可分割的 与 在传函形式上分开，并以 为过程控制通道的传递函数，以其输出信号作为反馈信号，以此来改善控制质量。预估补偿装置的传递函数为 ，使用后控制系统的框图如下， 系统闭环传递函数为 可见，经过Smith预估补偿，其特征方程中已经消去了 ，纯滞后环节对系统控制品质的不利影响已经消除。因此我们可以对过程中不包含纯滞后部分 进行控制器的设计，纯滞后部分 通过预估补偿器转换到闭环以外。这里选用PI控制器进行系统设计，按照前一题中介绍的方法，可确定PI控制器的传递函数为 在Matlab/Simulink环境下建立过程控制系统的仿真模型如下： 在给定输入为阶跃输入 时，系统的输出响应如下图所示，可见系统具有非常好的动态性能。 在理论上，Smith预估补偿可以克服大滞后的不利影响，但是Smith预估器需要知道被控对象精确的数学模型，并且对模型的误差非常敏感，因而在实际应用中仍面临很多问题。当被控对象的滞后时间常数变为 时，仿真可以得到系统此时的输出响应如下，可见系统的动态性能已经变的很差了。因此使用Smith 预估补偿的控制系统，其鲁棒性较差。 5-8如图5-72所示一带有库仑摩擦的二阶随动系统，试优化设计 参数，并分析非线性环节对系统动态性能的影响。 解：本题选用ITSE准则为目标函数。 首先，编写M函数文件，存为optm2.m。 function ss=optm2(x) global k; k=x; [tt,xx,yy]=sim('opt2',10,[]); ss=yy(end); 其次，输入恰当的模块参数，使用SIMULINK建立仿真模型文件，存为opt2.mdl。 在MATLAB命令窗口中输入下列指令: global k; result=fminsearch('optm2',[1]) 运行结果为： result =2.0182 用SIMULINK仿真验证所的结果，阶跃响应曲线如下： 与初始值 k＝1时的阶跃响应比较（下图），显然优化设计后，动态性能得到了很好的改善。 5-9 试分析图5-73所示系统中死区非线性对系统动态性能的影响。 解： 因为阶跃给定是10，如果c大于10，系统输出响应就为零，研究它已经没有意义，因此死区特性参数c应满足 。在Simulink环境下搭建系统模型如下： 分别在c=0(此时为倍数为一的比例增益)，2，4，6，8的情况下，仿真得到系统的输出曲线如下： 由此可以看出，引入非线性环节后随着c值的增加，系统动态响应速度变慢，系统超调减少，阻尼系数变大，当c=8时系统的输出响应已经表现出过阻尼特性。 5-10如图5-37所示计算机控制系统，试设计一最小拍控制器 ,并用仿真的方法分析最小拍控制器对系统输入信号和对象参数变化的适应性。 解：在此题中，为表述方便，设 表示系统闭环脉冲传递函数： 表示系统误差脉冲传递函数： 表示数字控制器脉冲传递函数； 表示控制对象的z变换； 表示输入信号的z变换； 表示数字控制器输入信号的z变换。        将被控对象离散化后可以得到 本题中，T=0.2，那么 也可用下列MATLAB命令语句求出： sysc=zpk([],[0 -1],[10]): c2d(sysc,0.2) 得到结果为 ，可以看出，与理论计算结果相同，但显然用MATLAB命令语句求解更加简单方便。 （1） 首先按照快速有纹波最小拍控制系统设计。 控制对象的传递函数为 ，由此可见，控制对象中含有单位圆上的极点p=1,为满足稳定性要求， 必须包含z=1的零点；因 滞后一拍，闭环脉冲传递函数 至少也要滞后一拍。系统针对等速输入设计，综上所述，设 = 。 为满足稳态精度，列写下列关系： 解得， ,即 由此可得， = 建立SIMULINK仿真模型如下： 用SIMULINK仿真观察设计结果，得到系统的输出响应曲线如下图所示， 数字控制器的输出曲线如下图所示， 由上面两个图可见，快速有纹波系统其输出值跟随输入值后，在非采样时刻有明显的纹波，而数字控制器振荡收敛，且振荡持续时间很长。这是由于数字控制器 = 有一个p=-0.9356的极点存在，其值与p=-1接近，由此造成数字控制器输出振荡收敛。这种情况下，不仅造成非采样时刻有误差，而且浪费执行机构的功率，增加机械磨损。 因此，我们可以考虑按照无纹波最小拍控制系统设计。 （2） 按照无纹波最小拍控制系统设计。 控制对象 ，由此可见，控制对象中含有单位圆上的极点p=1,以及零点z=-0.9356。若想令输出无纹波，则 必须包含 的全部零点（不管其在z平面的位置）。同时，为满足稳定性要求， 必须包含z=1的零点；因 滞后一拍，闭环脉冲传递函数 至少也要滞后一拍。系统针对等速输入设计，综上所述，设 = 为满足稳态精度，列写下列关系： 解得， , 即 由此可得， = 建立SIMULINK仿真模型如下： 用SIMULINK仿真观察设计结果，得到系统的输出响应曲线如下图所示， 数字控制器的输出曲线如下图所示， 从上两图可见，和有纹波系统相比，虽然无纹波系统的输出跟随输入函数所需时间要多一个采样周期，但采样点之间不存在纹波，且数字控制器输出波形较好，控制质量较高。 (3)按照无纹波最小拍系统设计方法，分析最小拍控制器对系统输入信号和对象参数变化的适应性。 a． 对系统输入信号的适应性。 当输入为阶跃信号时，在Matlab/Simulink环境下建立下图所示的控制系统结构图。 仿真得到系统输出曲线如下： 仿真得到数字控制器输出曲线如下： 由此可见，当输入为阶跃输入时，控制器输出正负跳动很大，而且系统输出超调严重，达到了100%。 当输入为加速度信号时，在Matlab/Simulink环境下建立下图所示的控制系统结构图， 仿真得到系统输出曲线如下图所示,为检验控制系统的性能，将输入输出曲线绘制在一张图上。 可以看出，系统对于加速度输入时，由明显的稳态误差。 综上所述，针对一种典型的输入函数进行设计，得到的最小拍系统闭环脉冲传递函数，用于次数较低的输入函数时，系统将出现超调，响应时间也将增加；而用于次数较高的输入函数时，输出将产生稳态误差。最小拍控制器对系统输入信号适应性不好。 b． 对对象参数变化的适应性。 设控制对象的传递函数为 ，本题中 。为检验最小拍控制器对于对象参数变化的适应性，令 ，得到此时的控制系统输出曲线如下图所示， 可以看出，所设计的最小拍控制器对于对象参数的变化，表现出了一定的适应性。将参数b做较大范围的修改后，输出仍然能够跟踪给定，但是恢复时间略有增加。 5-11 为使图5-75所示系统不产生自激振荡，试分析a、b取值。 解：线性部分传递函数为 因此 令线性部分频率响应 的虚部为零，即 求得 曲线与负实轴交点处的频率为 rad/s，将其代入到 中，得到 因此 曲线与负实轴的交点坐标为（-1.333，j0）。 这可以通过Matlab相关命令加以验证，在Matlab窗口中输入下列命令： num=3; den=conv([1,0],conv([0.8,1],[1,1])); nyquist(num,den) 得到线性部分的奈奎斯特图如下图所示， 将图形局部放大后，同样可得出其与负实轴交点坐标为（-1.333，j0），与理论分析结果一致。 图中所示系统为具有死区的继电器非线性系统，这一非线性特性的描述函数为 据此可画出 轨迹，图中的拐点为 。根据非线性控制理论，若使系统不产生自激振荡，必须保证 轨迹与 轨迹不相交，因此需令 可求得 ，可取 ，此时满足不产生自激振荡条件。 5-12 已知某一地区在有病菌传染下的描述三种类型人数变化的动态模型为 式中，X1表示可能传染的人数；X2表示已经得病的人数；X3表示已经治愈的人数； 。试用仿真方法求未来20年内三种人人数的动态变化情况。 解：编写M函数，并且存为rs.m文件： function dx=rs(t,x) dx=[-0.001*x(1)*x(2);0.001*x(1)*x(2)-0.072*x(2);0.072*x(2)]; 在MATLAB命令窗口中输入下列指令： t_end=20;                                    %设定仿真终止时间为20年  x0=[620;10;70];                              %三种类型人数初值 [t,x]=ode45('rs',[0,t_end],x0);              plot(t,x(:,1),'b-',t,x(:,2),'r--',t,x(:,3),'k.') xlabel('时间（年）'); ylabel('三种类型人数（人）'); legend('X1可能传染的人数','X2已经得病的人数','X3已经治愈的人数'); title('20年内三种人人数动态变化曲线图') x;                                          %查看三种类型人数 floor(x)                                      %人数应取整数      运行程序后可得到20年内三种人人数动态变化曲线，表格中列出了每年的三种 类型人数。 年数 0 1 2 3 4 5 6 7 8 9 10 X1 620 611 598 575 539 487 418 337 254 182 124 X2 10 17 29 49 80 124 182 248 311 359 389 X3 70 70 72 75 79 87 98 113 133 158 185                         年数 11 12 13 14 15 16 17 18 19 20 X1 83 56 37 25 17 12 9 6 5 3 X2 402 400 390 375 356 336 316 297 278 260 X3 213 242 271 299 325 350 373 396 416 436                       5-13.对于高阶系统的设计问题，往往要进行降阶近似处理，并要验证近似效果。已知某高阶系统模型为 经简化处理后，模型等效为 试比较两个模型在单位阶跃信号作用下的响应情况，并分析近似效果。 解：在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]; den2=[1,10.114,12.31]; sys2=tf(num2,den2); step(sys1,'-',sys2,':r',15) legend('原模型', '简化后模型') ; grid on 程序运行结果如下： 从曲线中可以看出，降阶后系统响应无超调，调整时间缩短，系统响应变快，但是简化前后终值有差异。 5-14 针对如图5-76所示的一阶双摆系统，试讨论如下问题（两杆所用相同）： （1）控制量F（施加的作用力）能否在保持双摆不倒的前提下，实现小车的位置伺服控制？ （2）试给出你的具体实现方案。 （提示：说明系统的实物制作原理与控制策略） 解题提示：首先需要对一阶双摆系统进行数学建模，在此基础上验证系统的可控性，如果系统是完全可控的，那么就可以在保持双摆不倒的前提下，实现小车的位置伺服控制。有资料证明，只要两个摆质量不相等（如果材料相同，则长度不相等），那么就可在保持双摆不倒的前提下，实现小车的伺服控制。 参考文献： 1、井上昭，邓明聪，逸见知弘，植木信幸.小车二级摆摆起控制器设计.控制理论与应用,2004,21(5):709-716. 2、易建强.平行二级倒立摆的稳定控制.控制与决策，2003,18(4):504:506. 第五章习题 1．(1) 可以实现 ，即实现自由旋转摆的垂直倒立。 方案： 1) 首先，要选取输入输出量、状态变量，这里给出：输入：M(力矩)，输出： 通过拉格朗日方程建立旋转倒立摆的数学模型为： 其中： 2) 在 位置附近10度范围内对模型进行线性化处理 即 ，所有二次小项都为0； 3) 选取状态变量为 4) 线性化处理后得到线性模型 5) 选取闭环极点P，使系统稳定，利用线性全状态反馈实现稳定控制 利用matlab提供的place(A,B,P)即可求出反馈矩阵K。满足输入 6) 带入此控制矩阵，即可实现倒立摆直立。 (2) 要实现自由旋转摆，即 ，必须满足如下所示的轨迹图。 通过建模可以得到系统的精确数学模型，对模型的可控性进行分析，考察其可控性随运行位置的变化属性，我们可以得到系统的可控性范围，可以采用非线性控制方法(如能量控制、模糊控制等)来实现摆杆在可控性范围内的任意位置稳定。 再者，可以把 作为扰动信号来看待，仍然采用前面的设计控制器方案，只不过保留较大的稳定裕度，使得加入扰动后仍然能实现上面杆的竖直稳定状态。 请参考文献： [1]傅雪冬. Pendubot 的基于能量的非线性控制. 机器人. Vol(22), No6. [2]赖旭芝. 欠驱动系统Pendubot的模糊控制策略. 电路与系统学报. Vol(9), no(6). 2．实现电动机调速范围(1r/10h---1r/min)，设计电机驱动器、数字式给定电路及数字控制器：请参考      [1] 冯绍娟 步进电机驱动器的设计与实现； [2] 张文增 机器人微型电机驱动器的设计与实现； 3.这是“吊车防摆问题”，在驱动力F作用下，小车运动，带动货物运动，绳索方向与垂直方向夹角 ，小车位置 ，这就构成了控制问题：通过F来控制 ，使得 ， 逐渐靠近目标位置。 需要建立系统模型，以 为变量，可以采用拉格朗日方法(基于能量)建立其数学模型；在建模过程中，可以分成两种：一种是摆长 不变化，即定摆长消摆问题，这可以采用全状态反馈方法(参照第一题)来实现控制目标；另一种是摆长变化 ，这就属于非线性控制问题(当然也可以在工作点附近线性化处理后，采用线性控制方法)，但控制水平、质量不能满足要求，这时可以采用滑模变结构、模糊控制等一些非线性控制方法达到控制要求。 请参考文献： [1]王伟, 易建强. 基于滑模方法的桥式吊车系统的抗摆控制. 控制与决策. Vol(19), No(9). 4．这个问题可以抽象为弹性材料的防摆问题，如图所示，要使得小车行走过程中，弹性杆与竖直方向的夹角 ，同时还要实现位置x最快达到目标地点，这是一个最优控制问题。 首先要建立系统模型：需要参考材料力学(包含应力、扭矩、截面、材料等关系)，广义力学(关于能量、广义力、广义坐标)及基本动力学知识(转动惯量、力矩、结构分解等)。 其次，选择合适的控制方法。可以在工作点附近线性化系统，得到线性定常模型 利用最优控制方法(参考线性系统理论或现代控制理论)，设计控制器结构。可以利用matlab工具箱提供的lqr()函数来方便的求取控制器参数K，即 来实现控制目标。 请参考文献： [1]刘才山, 王建明. 柔性机械臂的动力学模型及PD控制。 [2]于景元 柔性臂振动系统角速度反馈控制的最优指数衰减问题. 控制与决策. Vol(15)，NO(2) [3]李晓平 一种机器人柔性臂的振动控制新方法. 5．这道题属于最优控制问题，或者称为电梯的最优调度问题。 32层楼3部电梯，控制目标：实现电梯以最快速度、最短距离达到乘梯者位置、电梯运行方向(上升、下降)的合理变更、乘梯者等候时间最短、同运行方向乘梯者的综合搭载问题等。 设计过程中，可以考虑把每层楼的梯外信号(上升申请、下降申请)、梯内信号(到达某层)、以及三部电梯的当前位置信号以及当前运行方向信号作为考查对象，来设计“电梯最佳调度系统”。 请参考文献： [1]唐桂忠, 张广明. 基于混合遗传算法层间交通模式下的电梯群控派梯新方法. 电气自动化 2005.4 [2]宗群. 一种电梯群控多目标调度方法. 系统工程理论与实践. 2001.11. [3]应柏安. 局内电梯调度问题与竞争算法. 航空计算技术. Vol(31), NO(2).