1、机电系统动态性能的计算机仿真4机电系统动态性能的计算机仿真41 概述 机电系统计算机仿真是目前对复杂机电系统进行分析的重要手段与方法。在进行机电系统分析综合与设计工作过程中，除了需要进行理论分析外，还要对系统的特性进行实验研究。系统性能指标与参数是否达到预期的要求？它的经济性能如何？这些都需要在系统设计中给出明确的结论。对于那些在实际调试过程中存在很大风险或实验费用昂贵的系统，一般不允许对设计好的系统直接进行实验，然而没有经过实验研究是不能将设计好的系统直接放到生产实际中去的，因此就必须对其进行模拟实验研究。当然在有些情况下可以构造一套物理模拟装置来进行实验，但这种方法十分费时而且费用又高，而
2、在有的情况下物理模拟几乎是不可能的。近年来随着计算机的迅速发展，采用计算机对机电系统进行数学仿真的方法已被人们采纳。所谓机电系统计算机仿真就是以机电系统的数学模型为基础，借助计算机对机电系统的动静态过程进行实验研究。这里讲的机电系统计算机仿真是指借助数字计算机实现对机电系统的仿真分析。这种实验研究的特点是：将实际系统的运动规律用数学表达式加以描述，它通常是一组常微分方程或差分方程，然后利用计算机来求解这一数学模型，以达到对系统进行分析研究的目的。 对机电系统进行计算机仿真的基本过程包括：首先建立系统的数学模型，因为数学模型是系统仿真的基本依据，所以数学模型极为重要。然后根据系统的数学模型建立相
3、应的仿真模型，一般需要通过一定的算法或数值积分方法对原系统的数学模型进行离散化处理，从而建立起相应的仿真模型，这是进行机电系统仿真分析的关键步骤；最后根据系统的仿真模型编制相应的仿真程序，在计算机上进行仿真实验研究并对仿真结果加以分析。 机电系统计算机仿真的应用与发展已经过了近40年的历程，进入20世纪80年代以来，随着微型计算机技术以及软件技术的飞速发展与广泛应用，使得机电系统计算机仿真获得了实质性的发展，并使其走进广大的机电系统生产、设计、研究的第一线。目前在该领域应用最为广泛的软件包就是美国Mathworks公司开发的MATLAB语言软件。虽然它最初并非是为机电系统仿真与设计开发的，但是
4、它的强大的矩阵运算功能与图形处理及绘制能力，以及在MATLAB平台上开发出来的各种工具箱，和面向结构图的SIMULINK 系统分析环境，为机电系统计算机仿真提供了强有力的软件工具，从而为机电系统计算机仿真与辅助设计开辟新的方法与途径提供了充分的条件。 考虑到部分读者不具备MATLAB基础知识，因此，在本书的附录一中对MATLAB语言基础知识作了简单的介绍。42面向微分方程的连续系统仿真与MATLAB实现 对于控制系统的最基本数学描述就是微分方程，因此有必要在介绍机电系统计算机仿真时对这种方法加以介绍。实际上在前面章节已经介绍了几种典型数值积分方法以及它们的MATLAB实现。数值积分的数学基础就
5、是微分方程（状态方程），如果一个系统以微分方程的形式加以描述，就可以利用前面的知识对该系统进行仿真研究，有关这部分内容本章就不再加以讨论了。这里将介绍如何利用MATLAB中为用户所提供的常微分方程解函数ode和对微分方程进行求解计算，以及在面向微分方程的系统仿真研究中如何应用M函数。421 基于ode函数的面向微分方程的系统仿真1. 常微分方程解函数ode的基本格式 （l）ode函数的基本分类。在 MATLAB 5.x中 ode函数可以分为两类五种形式。一类是应用于非刚性微分方程求解，如ode45；另一类属于求解刚性微分方程的ode解函数，如ode23s.他们的具体形式如下：1 ode45函数
6、用于求解非刚性微分方程它属于中阶龙格一库塔法；2 ode23函数用于求解非刚性微分方程，它属于底阶龙格-库塔法；3 ode113函数用于求解非刚性微分方程的变阶法；4 ode15s函数用于求解刚性微分方程的变阶法；5 ode23s函数用于求解刚性微分方程的低阶法。 （2）ode函数的基本格式。上还这五种 ode函数的格式基本上一样这里主要以最常用的ode45函数为例加以介绍，其结果多数也可用于其他ode函数。ode45的基本格式主要有以下几种。 t，y=ode45(F，tspan，y0) 其中：F是ode文件名字符串；tspan= T0 TFINAL为一个向量，它指定了积分的启始时刻和结束时刻
7、； y0为积分运算的初始条件。 该函数对微分方程系统从T0到TFINAL时间段进行积分。函数返回一个列向量。在解向量Y中的每一行与列向量T的返回时间相对应。为了在指定的时刻点T0，T1，TFINAL得到解，则使用tspan=T0，T1，TFINAL。 t , yode45（F， tspan， y0，options）该函数的解如同前一格式，只是在这里用选项options中的参数来代替前一格式中的缺省值，options中的参数变量由函数odeset予以建立。最常使用的选项参数是标量相对容差RelTol（缺省值为1e-3）与绝对容差向量AbsTol（缺省的全部向量元素均为1e-6）。 （3）odes
8、et函数的基本格式。odeset函数用来建立或改变options选项结构，其基本格式为 options=odeset（namel，valuel，name2，value2，）通过该odeset函数建立积分器选项options的结构，在options中所命名的属性具有指定值，对于任何未指定的属性均为缺省值。 optionsodeset（oldopts，namel， valuel，）通过该odeset函数改变现存选项结构oldopts。 options=odeset（oldopt， newopts）通过该odeset函数将现有选项结构oldopts与新选项结构newopts结合在一起，新选项中的属性
9、将覆盖对应原有属性。2机电系统仿真的ode函数实现由ode函数可知，当已知系统由微分方程加以描述时，并建立以F作为ODE文件名，则可应用 t，y=ode45(F，tspan，y0)实现对系统的求解分析。这里所讲的ODE文件实际就是M函数文件，因此对于以高阶微分方程加以描述的系统，首先应将其写成一组一阶微分方程的形式，当然如果系统是以状态方程描述的，即可直接对其编写M函数文件。函数文件的特征就是文件第一行一定是以function开始，而且在函数文件中的变量均为局部变量，同时由函数文件构造的功能函数可以像MAT-LAB的其他函数一样进行调用。实际上在ode函数t，yode45（F，tspan，y0
10、）中，对以F作为文件名的ODE文件就是这类函数文件。下面通过实例说明对系统仿真过程中的具体应用。 例4-1 之已知一三阶系统，其状态方程表达式为 其中 ， 输入为u= 100，试应用ode函数求解该系统。 解： 由于该系统是以状态方程形式给出的，因此即可直接对其编写M函数文件。 function xd= rr（t，x） u100； xd=-14*x(1)+9*x(2)+10*x(3)+2.8*u；12*x(1)-9*x(2)+10*x(3)+4*u；24*x(1)-24*x(2)-18*x(3)+12*u；该函数文件对所研究系统的状态方程加以描述并形成rr.m文件以便在ode解函数中使用后即可
11、利用ode解函数对系统进行求解分析。考虑状态的零初值条件，即 图4-1 系统状态曲线 clear x0=0；0；0； t，x=ode45(rr，0 3，x0)； plot(t，x)对其同样可存一M文件，如li32.m文件，这样即可在命令窗口执行该文件。在运行过程中它会自动调用相关 rr函数文件，实现预定的工作，图 4-1给出了相应的运行结果。根据ode函数用途可知，该函数应用于求解常微分方程，因此对于含有输出方程的系统状态空间表达式就不能直接利用ode函数求解系统输出响应。考虑到ode函数返回变量的向量格式，可以根据输出方程的具体要求，由状态变量的线性组合直接得到所要求的系统输出。利用ode函
12、数可以方便地对以微分方程形式描述的系统进行仿真分析，然而ode解函数并非是万能的，在应用ode函数时对其进行干预一般是不方便的，因此在很多情况下，对以微分方程形式给出的系统，还必须采用前面所介绍的各种数值积分方法进行仿真分析。422 基于 M函数的面向微分方程的系统仿真 M函数是MATLAB中所特有的一类函数，前面介绍ode函数时已经使用了M函数。这里将介绍采用M函数建立描述系统的数学模型并在此基础上编制相应的程序。根据需要调用该M文件，实现所要求的运算。下面通过实例加以说明。例4-2已知系统状态方程为 应用数值积分方法（预估校正梯形法）实现对系统状态的仿真研究。 解：首先建立描述该系统数学模
13、型的M函数文件，相应的程序如下。 function df=rrl（X，u） A= -0.5572 -0.7814； 0.7814 0； B=1；0 U= 20； df=A*X+B*u 将该函数文件存盘取作rrlm文件。然后根据预估校正梯形公式编制相应的程序。 h=0.0l；t0=0；tf=15； X=0 0； X1=0 0； T=(tf-t0)/h for i=l：T K1=rrl（X，1）； X= X+h*K1 K2=rrl（X，1）； X= X+h*（K1+K2）/2 X1=X1；X；endplot（X1）在该程序段中两次调用了rr1.m函数文件，在第一次调用rrl.m后，对各状态的计算将
14、产生新值，因此K2就是经过预估计算所得到的相关值,从而即可应用梯形法实现对系统状态的仿真分析。图4-2给出了系统状态响应仿真曲线。 图4-2 系统状态响应43 面向传递函数的连续系统仿真 一个控制系统的动态模型经常是一传递函数的形式给出，特别是当给出开环传递函数， 在进行系统仿真时就要将开环传递函数转换为开环状态方程，再由开环状态方程求出对应的闭环状态方程，然后，再应用数值积分方法进行分解和分析仿真。下面将就这类问题加以介绍并给出MATLAB的实现。431 面向传递函数的系统仿真 考虑一个连续系统，它的开环传递函数的一般描述为 其相应的状态空间表达式为 其中A、B、C为相应的系统系数阵，这些系
15、数阵可以写作能控标准型形式或能观标准型形式，当然也可以写成其他的形式，这些都不会影响系统分析的结果。现如果假定系统结构图如图4-3所示，对于SISO系统，其v、u、y 、f 均为标量。由结构图可知u = v-fy，将其代入状态空间表达式得到 图 4-3 系统结构图式中 （4-1）输出方程为： （4-2）（4-1）式即为系统闭环状态方程。因此可以方便地利用各种数值积分方法对其进行求解，实现对以传递函数形式给出的系统进行仿真分析研究。显然，如果系统的数学模型是以微分方程或状态方程形式给出的，一般就可以直接利用前述数值积分方法或ode函数方法对系统进行仿真研究。432 面向传递函数的系统仿真的MAT
16、LAB实现 下面通过一个例子说明如何利用MATLAB语言进行以传递函数形式给出的系统的仿真分析。 例4-3 设单位反馈系统的开环传递函数为 试给出系统闭环的阶跃响应。 解 根据开环传递函数求出系统开环状态方程，应用MATLAB控制工具箱中的函数，可以方便地求出系统开环状态方程。 num =2*1,2; den = conv(conv(1,0,1,3),1,2,2) G1 = tf(num,den); G2 = ss(G1); a,b,c,d = ssdata (G);这样即可以得到系统系数阵分别为 a = -5.0000 -2.0000 -0.7500 0 4.0000 0 0 0 0 2.0
17、000 0 00 0 1.0000 0 b = 1 0 0 0 c = 0 0 0.2500 0.5000 d = 0进而就可以方便地求出系统闭环状态方程 Ab = a-b*c; Bb = b; Cb = c Db = 0Ab = -5.0000 -2.0000 -1.0000 -0.5000 4.0000 0 0 0 0 2.0000 0 00 0 1.0000 0 Bb = 1 0 0 0 Cb = 0 0 0.2500 0.5000 Db = 0从而可以应用数值积分方法求取系统的阶跃响应。下面给出了应用MATLAB语言编写的四阶龙格-库塔数值积分方法的仿真程序，相应的仿真曲线由图4-4
18、给出。X = 0;0;0;0; Y = 0; t = 0; U = 25; Tf = 15; h=0.01; N = Tf/h; Ab = -5 -2 -1 -0.5; 4 0 0 0; 0 2 0 0; 0 0 1 0;图4-4 四阶龙格-库塔数值积分方法的仿真曲线 Bb = 1; 0; 0;0; Cb = 0 0 0.25 0.5; Db = 0; for i = 1:n K1 = Ab * X + Bb*U; K2 = Ab * (X + h*K1/2)+Bb*U;K3 = Ab * (X + h*K2/2)+Bb*U;K4 = Ab * (X + h*K3)+Bb*U; X = X+
19、h*(K1 + 2*K2+ 2*K3 +K4)/6; Y = Y, Cb*X; t = t, t(i)+h; end plot (t,Y) Grid 当由已知系统开环传递函数求取系统开环状态方程时，除了可以利用ss函数直接求取外，还可以使用状态方程相似变换函数ss2ss进行求取，下面给出上例中应用相似变换函数ss2ss求取系统开环状态方程的基本用法。这里主要利用矩阵水平翻转函数fliplr构造出反对角阵作为变换阵，从而得到相应的状态方程。 num = 2 *1,2; den = conv(conv(1,0,1,3),1,2,2); G1= tf (num, den); G2 = ss(G1);
20、 T = fliplr (eye(4); G = ss2ss(G2,T)运行结果： a = x1 x2 x3 x4 x1 0 1.00000 0 0 x2 0 0 2.00000 0 x3 0 0 0 4.00000 x4 0 -0.75000 -2.00000 -5.00000 b = u1x1 0 x2 0 x3 0 x4 1.00000 c = x1 x2 x3 x4 y1 0.50000 0.25000 0 0 d = u1 y1 0 从计算结果可以看出，该状态方程十分类似通常定义的系统可控标准型形式。 当然如果在研究系统时，需要状态方程以可控标准型形式给出，则可应用MATLAB编程得
21、到。但不论哪种形式的状态方程，最终所要研究的系统阶跃响应都是一样的。44面向结构图的连续系统仿真441 基于典型环节的系统仿真 由于习惯上往往采用结构图形式对系统进行描述，因此在对系统进行仿真时，有必要研究如何以结构图为基础进行。在此将介绍以典型环节为基本环节的面向结构图的系统仿真方法，这种方法的主要优点是： 便于研究各环节参数对系统的影响 可以得到每个环节的动态响应； 对于存在系统非线性环节较容易处理； 不但适合单人单出系统，而且对于多入多出系统也适用。 该方法只要建立起各个典型环节的以及各环节之间的连接关系，就能方便地求取任意环节和系统的动态响应。由于直接输入各环节的参数，因此研究参数变化
22、对系统的影响也十分方便。1典型环节的确定 控制系统常见的环节有：比例环节 K、积分环节、惯性环节、一阶超前（滞后）环节、二阶振荡环节等。为了仿真的方便，必须选择一种能够代表这些环节的典型环节。经过分析可以看出，只要采用下述形式的环节作为典型环节，就可以满足对系统研究的要求，即典型模型可表示为 （4-3）显然，当典型环节中时即为惯性环节；时即是积分环节；至于振荡环节则只要用两个典型环节串联并加上一个负反馈即可得到。其中典型环节中的、均为对角阵，其对角线元素值由各个环节中的相关参数确定。根据典型环节可以得到第i个环节的输出为 因此得到 （i=1，2，n） 将其写为矩阵方程形式，得到 （4-4） 其
23、中 ， B= ，2系统连接矩阵 前面给出了从微观角度对系统环节的描述，然而为了研究一个系统，还必须从宏观角度对系统加以描述，这就需要将各个环节之间的连接关系用一种数学形式加以表示。考察图4-5所示系统，该系统由4个典型环节构成，根据系统结构图可以得到各环节的输入输出关系方程组。 图4-5系统结构图 将其写作矩阵形式 （4-5）其中 ， ， ， 一般称与为系统连接矩阵。为一方阵，元素表示结构图中第j 个环节对第i个环节的作用；若第j个环节对第i个环节没有作用，则。表示外部输入信号对系统的连接情况。3周环系统的状态方程 将(4-5)式代入(4-4)式，消去U得 对上式展开，整理并取拉氏反变换得 考
24、虑在外输入为阶跃信号时，避免对应的初始导数无穷大，则要求。如果令 （4-6） （4-7）可以得到 如果的逆存在，则得到 若令 ， （4-8）则上式可写作 （4-7）显然上式即为一典型的状态方程，从而可利用前面介绍的数值积分方法对其进行求解，实现对系统的仿真研究。下面以一个结构图形式的系统为例，介绍基于MATLAB的面向结构图的系统仿真实现。 例4-4 已知系统结构图如图4-5所示，采用MATLAB实现面向结构图的系统仿真 。 解：根据系统结构图首先利用MATLAB建立各个环节参数构成的系数阵（此处的系数阵均为对角阵）与系统连接阵，在此基础上再依（4-6）式、（4-7）式、（4-8）式建立闭环系
25、统的相应系数阵。下面给出在MATLAB环境下的相应程序。 AA=diag（1 6.5 1 1）； BB=diag（4 1 3 2 ）; CC=diag（4 1 2.4 5）； DD=diag（0 0 0 1）； W=0 0 -4 0;1 0 0 -1;0 1 0 0 ;0 0 1 0; W0=1 0 0 0 ； Q1BB-DD*W； Qinv(Q1）； RCC*WAA； AQ*R； BQ*CC*W0； 一旦得到了闭环系统的系数阵，即可应用数值积分方法实现对系统的分析仿真。应注意，在该状态方程中各个状态变量的实际含义。这里应用四阶龙格一库塔法对其进行仿真分析，下面给出了相应的仿真程序。图4-6给
26、出了采用面向结构图系统阶跃响应曲线。图4-6面向结图仿真（阶跃响应） 图4-7环节参数改变后面向结构图的仿真 X0；0；0；0；Y0；t=0；U25； Tf=15；h001； n=Tfh； Cb=0 01 0； for i 1:n K1A*X+B*U； K2A*（Xh*K1/2） B*U； K3=A*（X h*K2/2） B*U； K4=A*（X h*K3） B*U； XX+h*(K12*K22*K3K4）/6； Y=Y，Cb*X； t=t,t(i)+h; end plot(t,y) grid在这里由于采用的是面向结构图形式的系统仿真，因此研究因各环节参数发生变化对系统的影响就十分方便，仅对环
27、节系数阵加以修改即可。如果将第三个环节修改为2.4/( s +1)时，只需修改B阵的相关元素即可，其他部分无需改动，图4-7给出了改变环节三参数后的仿真结果。由响应曲线看出，系统的响应快速性明显提高了，同时超调量也有所降低。另外采用面向结构图仿真，还可十分方便地研究各个环节输出的响应情况，只要根据需要修改相应系数阵即可。图4-8（a）、（b）和（c） 图4-8 环节 1（a）、环节2（b）、环节4（c）的响应分别给出了环节1、环节2和环节4的输出响应。442 基于Connect连接函数的系统仿真 除了前面介绍的面向结构图的系统仿真方法外，还可利用MATLAB与控制工具箱所提供的相关函数与命令，
28、实现以结构图形式给出的系统仿真。下面通过实例说明如何以这种方法进行仿真工作。 例4-5 已知系统如图4-9所示，研究该系统的阶跃响应。 解： 显然该系统是以结构图形式给出的较复杂系统；利用MATLAB和控制工具箱可以方便地实现系统仿真。就基于Connect连接函数的仿真而言，一般可分为四步进行。1建立系统中间状态模型 首先应对系统各环节进行编号（图4-9），然后即可据此与给定系统参数建立有关中间模型，以便运行时在MATLAB的工作空间建立相应的系统环节模型。其基本语句和格式为 nblocks10； nl=1；dl= 001 1； n2=017 1；d2=0085 0; n31；d3=001 1
29、； n4=015 1；d4=0051 0； n5=70；d5=00067 1； 图4-9系统结构图 n6= 021；d6= 015 1； n7= 130；d7=1 0; n8=-0212；d81； n9= - 01 ；d90011 ； n10=- 00044；d10= 0 01 1； blkbuild； 其中ni、di表示第i个环节传递函数的分子、分母多项式的系数向量，如果所研究的系统是以状态方程形式给出的各个环节，则应采用ai、bi、ci、di来表示环节状态方程的系数阵。对于每个环节只能取二者之一，MATLAB会自动根据所输入的阵列元素数量加以判定并进行相应的运算。nblocks需根据系统环节个数对其进行赋值， blkbuild语句将对nblocks语句后的各环节输入情况进行校验，并自动生成系统的隐含中间状态方程模型（a、b、c、d系数阵）。2建立系统连接矩阵 在MATLAB工作空间建立了系统各环节的模型后，还应该建立反映各环节之间相互连接关系的Q阵，以及输入、输出向量（inputs、outputs)，下面结合实例将Q阵与inputs、outputs列写规则加以介绍。根据已知系统可以得到系统连接阵Q和输入、输出向量： Q1 0 0； 2 1 10； 3 2 0； 4 3 9； 5 4 0； 6 5 8； 7 6 0；
