《计算机辅助设计——机电系统CAD》教学讲义教材：第四章 机电系统动态性能的计算机仿真

4.机电系统动态性能的计算机仿真 4.1概述 机电系统计算机仿真是目前对复杂机电系统进行分析的重要手段与方法。在进 行机电系统分析综合与设计工作过程中,除了需要进行理论分析外,还要对系统的 特性进行实验研究。系统性能指标与参数是否达到预期的要求?它的经济性能如 何?这些都需要在系统设计中给出明确的结论。对于那些在实际调试过程中存在很 大风险或实验费用昂贵的系统,一般不允许对设计好的系统直接进行实验,然而没 有经过实验研究是不能将设计好的系统直接放到生产实际中去的,因此就必须对其 进行模拟实验硏究。当然在有些情况下可以构造一套物理模拟装置来进行实验,但 这种方法十分费时而且费用又高,而在有的情况下物理模拟几乎是不可能的。近年 来随着计算机的迅速发展,采用计算机对机电系统进行数学仿真的方法已被人们采 纳。所谓机电系统计算机仿真就是以机电系统的数学模型为基础,借助计算机对机 电系统的动静态过程进行实验研究。这里讲的机电系统计算机仿真是指借助数字计 算机实现对机电系统的仿真分析。这种实验研究的特点是:将实际系统的运动规律 用数学表达式加以描述,它通常是一组常微分方程或差分方程,然后利用计算机来 求解这一数学模型,以达到对系统进行分析研究的目的 对机电系统进行计算机仿真的基本过程包括:首先建立系统的数学模型,因为 数学模型是系统仿真的基本依据,所以数学模型极为重要。然后根据系统的数学模 型建立相应的仿真模型,一般需要通过一定的算法或数值积分方法对原系统的数学 模型进行离散化处理,从而建立起相应的仿真模型,这是进行机电系统仿真分析的 关键步骤;最后根据系统的仿真模型编制相应的仿真程序,在计算机上进行仿真实 验研究并对仿真结果加以分析。 机电系统计算机仿真的应用与发展已经过了近40年的历程,进入20世纪80 年代以来,随着微型计算机技术以及软件技术的飞速发展与广泛应用,使得机电系 统计算机仿真获得了实质性的发展,并使其走进广大的机电系统生产、设计、研究 的第一线。目前在该领域应用最为广泛的软件包就是美国 Mathworks公司开发的 MATLAB语言软件。虽然它最初并非是为机电系统仿真与设计开发的,但是它的强 大的矩阵运算功能与图形处理及绘制能力,以及在 MATLAB平台上开发出来的各种 工具箱,和面向结构图的 SIMULINK系统分析环境,为机电系统计算机仿真提供了 强有力的软件工具,从而为机电系统计算机仿真与辅助设计开辟新的方法与途径提 供了充分的条件。 考虑到部分读者不具备 MATLAB基础知识,因此,在本书的附录一中对 119

119 4．机电系统动态性能的计算机仿真 4．1 概述 机电系统计算机仿真是目前对复杂机电系统进行分析的重要手段与方法。在进 行机电系统分析综合与设计工作过程中，除了需要进行理论分析外，还要对系统的 特性进行实验研究。系统性能指标与参数是否达到预期的要求？它的经济性能如 何？这些都需要在系统设计中给出明确的结论。对于那些在实际调试过程中存在很 大风险或实验费用昂贵的系统，一般不允许对设计好的系统直接进行实验，然而没 有经过实验研究是不能将设计好的系统直接放到生产实际中去的，因此就必须对其 进行模拟实验研究。当然在有些情况下可以构造一套物理模拟装置来进行实验，但 这种方法十分费时而且费用又高，而在有的情况下物理模拟几乎是不可能的。近年 来随着计算机的迅速发展，采用计算机对机电系统进行数学仿真的方法已被人们采 纳。所谓机电系统计算机仿真就是以机电系统的数学模型为基础，借助计算机对机 电系统的动静态过程进行实验研究。这里讲的机电系统计算机仿真是指借助数字计 算机实现对机电系统的仿真分析。这种实验研究的特点是：将实际系统的运动规律 用数学表达式加以描述，它通常是一组常微分方程或差分方程，然后利用计算机来 求解这一数学模型，以达到对系统进行分析研究的目的。 对机电系统进行计算机仿真的基本过程包括：首先建立系统的数学模型，因为 数学模型是系统仿真的基本依据，所以数学模型极为重要。然后根据系统的数学模 型建立相应的仿真模型，一般需要通过一定的算法或数值积分方法对原系统的数学 模型进行离散化处理，从而建立起相应的仿真模型，这是进行机电系统仿真分析的 关键步骤；最后根据系统的仿真模型编制相应的仿真程序，在计算机上进行仿真实 验研究并对仿真结果加以分析。 机电系统计算机仿真的应用与发展已经过了近 40 年的历程，进入 20 世纪 80 年代以来，随着微型计算机技术以及软件技术的飞速发展与广泛应用，使得机电系 统计算机仿真获得了实质性的发展，并使其走进广大的机电系统生产、设计、研究 的第一线。目前在该领域应用最为广泛的软件包就是美国 Mathworks 公司开发的 MATLAB 语言软件。虽然它最初并非是为机电系统仿真与设计开发的，但是它的强 大的矩阵运算功能与图形处理及绘制能力，以及在 MATLAB 平台上开发出来的各种 工具箱，和面向结构图的 SIMULINK 系统分析环境，为机电系统计算机仿真提供了 强有力的软件工具，从而为机电系统计算机仿真与辅助设计开辟新的方法与途径提 供了充分的条件。 考虑到部分读者不具备 MATLAB 基础知识，因此，在本书的附录一中对

MATLAB语言基础知识作了简单的介绍 4.2面向微分方程的连续系统仿真与 MATLAB实现 对于控制系统的最基本数学描述就是微分方程,因此有必要在介绍机电系统计 算机仿真时对这种方法加以介绍。实际上在前面章节已经介绍了几种典型数值积分 方法以及它们的 MATLAB实现。数值积分的数学基础就是微分方程(状态方程), 如果一个系统以微分方程的形式加以描述,就可以利用前面的知识对该系统进行仿 真研究,有关这部分内容本章就不再加以讨论了。这里将介绍如何利用 MATLAB中 为用户所提供的常微分方程解函数ode和对微分方程进行求解计算,以及在面向微 分方程的系统仿真研究中如何应用M函数 4.2.1基于ode函数的面向微分方程的系统仿真 1.常微分方程解函数ode的基本格式 (1)ode函数的基本分类。在 MATLAB5x中ode函数可以分为两类五种形式。 类是应用于非刚性微分方程求解,如ode45:另一类属于求解刚性微分方程的ode 解函数,如ode23s他们的具体形式如下 lodε45函数用于求解非刚性微分方程.它属于中阶龙格一库塔法 2ode23函数用于求解非刚性微分方程,它属于底阶龙格-库塔法 3 odell3函数用于求解非刚性微分方程的变阶法 4odel5s函数用于求解刚性微分方程的变阶法 5ode23s函数用于求解刚性微分方程的低阶法 (2)ode函数的基本格式。上还这五种ode函数的格式基本上一样.这里主要 以最常用的ode45函数为例加以介绍,其结果多数也可用于其他ode函数。ode45的 基本格式主要有以下几种 O [t, yI=ode45('F,tspan, yO 其中:F是ode文件名字符串 tspan=[10 TFINAL]为一个向量,它指定了积分的启始时刻和结束时刻 y0为积分运算的初始条件。 该函数对微分方程系统y=F(y从m0到TNAL时间段进行积分。函数F(y) 返回一个列向量。在解向量Y中的每一行与列向量T的返回时间相对应。为了在指 定的时刻点T0,T1,…, TFINAL得到解,则使用 tspan=[T0,Tl,…, TFINAL] 2 t, y]=ode45('F, tspan, yO, options) 120

120 MATLAB 语言基础知识作了简单的介绍。 4．2 面向微分方程的连续系统仿真与 MATLAB 实现 对于控制系统的最基本数学描述就是微分方程，因此有必要在介绍机电系统计 算机仿真时对这种方法加以介绍。实际上在前面章节已经介绍了几种典型数值积分 方法以及它们的 MATLAB 实现。数值积分的数学基础就是微分方程（状态方程）， 如果一个系统以微分方程的形式加以描述，就可以利用前面的知识对该系统进行仿 真研究，有关这部分内容本章就不再加以讨论了。这里将介绍如何利用 MATLAB 中 为用户所提供的常微分方程解函数 ode 和对微分方程进行求解计算，以及在面向微 分方程的系统仿真研究中如何应用 M 函数。 4．2．1 基于 ode 函数的面向微分方程的系统仿真 1. 常微分方程解函数 ode 的基本格式 （l）ode 函数的基本分类。在 MATLAB 5.x 中 ode 函数可以分为两类五种形式。 一类是应用于非刚性微分方程求解，如 ode45；另一类属于求解刚性微分方程的 ode 解函数，如 ode23s.他们的具体形式如下： 1 ode45 函数用于求解非刚性微分方程．它属于中阶龙格一库塔法； 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］为一个向量，它指定了积分的启始时刻和结束时刻； y0 为积分运算的初始条件。 该函数对微分方程系统 y  = F(t, y) 从 T0 到 TFINAL 时间段进行积分。函数 F(t, y) 返回一个列向量。在解向量 Y 中的每一行与列向量 T 的返回时间相对应。为了在指 定的时刻点 T0，T1，…，TFINAL 得到解，则使用 tspan=[T0，T1，…，TFINAL]。 ② [t , y]＝ode45（'F'， tspan， y0，options）

该函数的解如同前一格式,只是在这里用选项 options中的参数来代替前一格式中的 缺省值, options中的参数变量由函数 odeset予以建立。最常使用的选项参数是标量 相对容差 RelTo"(缺省值为le-3)与绝对容差向量' Abs To'(缺省的全部向量元素均 为1e-6) (3) odeset函数的基本格式。 odeset函数用来建立或改变 options选项结构,其 基本格式为 ptions=odeset (,, valuel, 'name2,, value 通过该 odeset函数建立积分器选项 options的结构,在 options中所命名的属性具有 指定值,对于任何未指定的属性均为缺省值。 options=odeset (oldopts, 'namel,, valuel,.) 通过该 odeset函数改变现存选项结构 adopts 通过该 odeset函数将现有选项结构 adopts与新选项结构 newts结合在一起,新选 项中的属性将覆盖对应原有属性 2.机电系统仿真的ode函数实现 由od函数可知,当已知系统由微分方程y=F(y)加以描述时,并建立以F 作为ODE文件名,则可应用 It, yI=ode45(F, tspan, yo) 实现对系统的求解分析。这里所讲的ODE文件实际就是M函数文件,因此对于以 高阶微分方程加以描述的系统,首先应将其写成一组一阶微分方程的形式,当然如 果系统是以状态方程描述的,即可直接对其编写M函数文件。函数文件的特征就是 文件第一行一定是以 function开始,而且在函数文件中的变量均为局部变量,同时 由函数文件构造的功能函数可以像 MAT-LAB的其他函数一样进行调用。实际上在 ode函数[t,y]=ode45(F, tspan,y0)中,对以F作为文件名的ODE文件就是这 类函数文件。下面通过实例说明对系统仿真过程中的具体应用。 例4-1之已知一三阶系统,其状态方程表达式为 X=AX+Bl 其中 4910 A=12-910B=4 24-24-18 12 输入为u=100,试应用ode函数求解该系统

121 该函数的解如同前一格式，只是在这里用选项 options 中的参数来代替前一格式中的 缺省值，options 中的参数变量由函数 odeset 予以建立。最常使用的选项参数是标量 相对容差'RelTol'（缺省值为 1e-3）与绝对容差向量'AbsTol'（缺省的全部向量元素均 为 1e-6）。 （3）odeset 函数的基本格式。odeset 函数用来建立或改变 options 选项结构，其 基本格式为 options=odeset（'namel'，valuel，'name2'，value2，…） 通过该 odeset 函数建立积分器选项 options 的结构，在 options 中所命名的属性具有 指定值，对于任何未指定的属性均为缺省值。 options＝odeset（oldopts，'namel'， valuel，…） 通过该 odeset 函数改变现存选项结构 oldopts。 options=odeset（oldopt， newopts） 通过该 odeset 函数将现有选项结构 oldopts 与新选项结构 newopts 结合在一起，新选 项中的属性将覆盖对应原有属性。 2．机电系统仿真的 ode 函数实现 由 ode 函数可知，当已知系统由微分方程 y  = F(t, y) 加以描述时，并建立以'F' 作为 ODE 文件名，则可应用 [t，y]=ode45('F'，tspan，y0) 实现对系统的求解分析。这里所讲的 ODE 文件实际就是 M 函数文件，因此对于以 高阶微分方程加以描述的系统，首先应将其写成一组一阶微分方程的形式，当然如 果系统是以状态方程描述的，即可直接对其编写 M 函数文件。函数文件的特征就是 文件第一行一定是以 function 开始，而且在函数文件中的变量均为局部变量，同时 由函数文件构造的功能函数可以像 MAT-LAB 的其他函数一样进行调用。实际上在 ode 函数[t，y]＝ode45（'F'，tspan，y0）中，对以'F'作为文件名的 ODE 文件就是这 类函数文件。下面通过实例说明对系统仿真过程中的具体应用。 例 4-1 之已知一三阶系统，其状态方程表达式为 X  = AX + BU 其中           − − − − = 24 24 18 12 9 10 4 9 10 A ，           = 12 4 2.8 B 输入为 u= 100，试应用 ode 函数求解该系统

解:由于该系统是以状态方程形式给出的,因此即可直接对其编写M函数文 件 +10*x(3)+4*u;24*x(1)-24*x(2)-18*x(3)+12*山 该函数文件对所研究系统的状态方程加以描述并形成rm文件以便在ode解函数中 使用后即可利用ode解函数对系统进行求解分析。考虑状态的零初值条件,即 0.5 1.5 2.5 图4-1系统状态曲线 [t,x]=ode45(r,[03],xO) plot(t, x) 对其同样可存一M文件,如li32m文件,这样即可在命令窗口执行该文件。在运行 过程中它会自动调用相关m函数文件,实现预定的工作,图41给出了相应的运行 结果 根据ode函数用途可知,该函数应用于求解常微分方程,因此对于含有输出方 程的系统状态空间表达式就不能直接利用ode函数求解系统输出响应。考虑到ode 函数返回变量的向量格式,可以根据输出方程的具体要求,由状态变量的线性组合 直接得到所要求的系统输出。利用ode函数可以方便地对以微分方程形式描述的系 统进行仿真分析,然而ode解函数并非是万能的,在应用ode函数时对其进行干预 般是不方便的,因此在很多情况下,对以微分方程形式给出的系统,还必须采用 前面所介绍的各种数值积分方法进行仿真分析。 122

122 解： 由于该系统是以状态方程形式给出的，因此即可直接对其编写 M 函数文 件。 function xd= rr（t，x） u＝100； 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 解函数中 使用后即可利用 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 函数可以方便地对以微分方程形式描述的系 统进行仿真分析，然而 ode 解函数并非是万能的，在应用 ode 函数时对其进行干预 一般是不方便的，因此在很多情况下，对以微分方程形式给出的系统，还必须采用 前面所介绍的各种数值积分方法进行仿真分析

4.2.2基于M函数的面向微分方程的系统仿真 M函数是 MATLAB中所特有的一类函数,前面介绍ode函数时已经使用了M 函数。这里将介绍采用M函数建立描述系统的数学模型并在此基础上编制相应的程 序。根据需要调用该M文件,实现所要求的运算。下面通过实例加以说明 例42已知系统状态方程为 05572-07814x 0.7814 0 应用数值积分方法(预估校正梯形法)实现对系统状态的仿真研究 解:首先建立描述该系统数学模型的M函数文件,相应的程序如下 function dfrrl (X, u) A=[0.5572-0.7814:0.78140]: df=A*X+B*u 将该函数文件存盘取作rl.m文件。然后根据预估校正梯形公式编制相应的程序。 h=00l;t0=0;tf=15 X=[00];X1=00 T=(tf-tO)h for i=l: T KI=r(X,1); X= X+h*KI K2=rrl (X, 1): X=X+h*(K1+K2)/2 XI=[XI: XI plot(XI) 在该程序段中两次调用了rr1m函数文件,在第一次调用rlm后,对各状态的 计算将产生新值,因此K2就是经过预估计算所得到的相关值,从而即可应用梯形法 实现对系统状态的仿真分析。图4-2给出了系统状态响应仿真曲线。 123

123 4．2．2 基于 M 函数的面向微分方程的系统仿真 M 函数是 MATLAB 中所特有的一类函数，前面介绍 ode 函数时已经使用了 M 函数。这里将介绍采用 M 函数建立描述系统的数学模型并在此基础上编制相应的程 序。根据需要调用该 M 文件，实现所要求的运算。下面通过实例加以说明。 例 4-2 已知系统状态方程为 u x x x x       +           − − =      0 1 0 7814 0 0 5572 0 7814 2 1 2 1 . . .   应用数值积分方法（预估校正梯形法）实现对系统状态的仿真研究。 解：首先建立描述该系统数学模型的 M 函数文件，相应的程序如下。 function df=rrl（X，u） A= ［-0.5572 -0.7814； 0.7814 0］； B=[1；0] U= 20； df=A*X+B*u 将该函数文件存盘取作 rrl．m 文件。然后根据预估校正梯形公式编制相应的程序。 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']； end plot（X1） 在该程序段中两次调用了 rr1.m 函数文件，在第一次调用 rrl.m 后，对各状态的 计算将产生新值，因此 K2 就是经过预估计算所得到的相关值,从而即可应用梯形法 实现对系统状态的仿真分析。图 4-2 给出了系统状态响应仿真曲线

10 00 1000 图4-2系统状态响应 4.3面向传递函数的连续系统仿真 一个控制系统的动态模型经常是一传递函数的形式给出,特别是当给出开环传 递函数,在进行系统仿真时就要将开环传递函数转换为开环状态方程,再由开环状 态方程求出对应的闭环状态方程,然后,再应用数值积分方法进行分解和分析仿真 下面将就这类问题加以介绍并给出 MATLAB的实现。 4.3.1面向传递函数的系统仿真 考虑一个连续系统,它的开环传递函数的一般描述为 CrS+Cs-2 +cn-Istc 其相应的状态空间表达式为 X=AX+Bl Y=CX 其中A、B、C为相应的系统系数阵,这些系数阵可以写作能控标准型形式或能观标 准型形式,当然也可以写成其他的形式,这些都不会影响系统分析的结果。现如果 假定系统结构图如图4-3所示,对于SISO系统,其v、u、y、f均为标量。由结构 图可知u=v-fy,将其代入状态空间表达式得到 X=AX+B(u-fy AX+B

124 图 4-2 系统状态响应 4．3 面向传递函数的连续系统仿真 一个控制系统的动态模型经常是一传递函数的形式给出，特别是当给出开环传 递函数， 在进行系统仿真时就要将开环传递函数转换为开环状态方程，再由开环状 态方程求出对应的闭环状态方程，然后，再应用数值积分方法进行分解和分析仿真。 下面将就这类问题加以介绍并给出 MATLAB 的实现。 4．3．1 面向传递函数的系统仿真 考虑一个连续系统，它的开环传递函数的一般描述为 ( ) ( ) ( ) n n n n n n n n n s a s a s a s a c s c s c s c U s Y s G s + + + + + + + + + = = − − − − − − 1 2 2 1 1 1 2 2 1 1   其相应的状态空间表达式为 Y CX X AX BU =  = + 其中 A、B、C 为相应的系统系数阵，这些系数阵可以写作能控标准型形式或能观标 准型形式，当然也可以写成其他的形式，这些都不会影响系统分析的结果。现如果 假定系统结构图如图 4-3 所示，对于 SISO 系统，其 v、u、y 、f 均为标量。由结构 图可知 u = v-fy，将其代入状态空间表达式得到 A X Bv X AX B v fy = b + = + ( − ) 

t= Ax Bu C 图4-3系统结构图 式中A=A-BC 输出方程为: Y=CX (4-1)式即为系统闭环状态方程。因此可以方便地利用各种数值积分方法对其进行 求解,实现对以传递函数形式给出的系统进行仿真分析研究。显然,如果系统的数 学模型是以微分方程或状态方程形式给出的,一般就可以直接利用前述数值积分方 法或ode函数方法对系统进行仿真研究。 4.3.2面向传递函数的系统仿真的 MATLAB实现 下面通过一个例子说明如何利用 MATLA B语言进行以传递函数形式给出的系 统的仿真分析。 例4-3设单位反馈系统的开环传递函数为 G(S) 2(S+2) (S+3)s2+2s+2 试给出系统闭环的阶跃响应。 解根据开环传递函数求出系统开环状态方程,应用 MATLAB控制工具箱中的 函数,可以方便地求出系统开环状态方程。 num=2*[12];den=conv(conv([1,0][1,3]1,22]) la,b, c, d =ssdata(G)

125 图 4-3 系统结构图 式中 Ab = A− BfC （4-1） 输出方程为： Y =CX （4-2） （4-1）式即为系统闭环状态方程。因此可以方便地利用各种数值积分方法对其进行 求解，实现对以传递函数形式给出的系统进行仿真分析研究。显然，如果系统的数 学模型是以微分方程或状态方程形式给出的，一般就可以直接利用前述数值积分方 法或 ode 函数方法对系统进行仿真研究。 4．3．2 面向传递函数的系统仿真的 MATLAB 实现 下面通过一个例子说明如何利用 MATLAB 语言进行以传递函数形式给出的系 统的仿真分析。 例 4-3 设单位反馈系统的开环传递函数为 ( 3)( 2 2) 2( 2) ( ) 2 + + + + = s s s s s G s 试给出系统闭环的阶跃响应。 解 根据开环传递函数求出系统开环状态方程，应用 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);

这样即可以得到系统系数阵分别为 5.0000-2.0000-0.75000 02.0000 01.00000 00000 进而就可以方便地求出系统闭环状态方程 Ab=a-b*c: bb=b: Cb=c Db=0 Ab 50000-2.0000-1.0000-0.5000 4.0000 02.0000 0 0 01.0000 0 Bb 000 00.25000.5000 从而可以应用数值积分方法求取系统的阶跃响应。下面给出了应用 MATLAB语言编 写的四阶龙格-库塔数值积分方法的仿真程序,相应的仿真曲线由图4-4给出 X=[00,0,0]Y=0;t=0,U=25,Tf=15;h=0.01; N=Tf/h Ab=[-5-2-1-0.5;4000,0200,00101 126

126 这样即可以得到系统系数阵分别为 a = -5.0000 -2.0000 -0.7500 0 4.0000 0 0 0 0 2.0000 0 0 0 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 = 0 Ab = -5.0000 -2.0000 -1.0000 -0.5000 4.0000 0 0 0 0 2.0000 0 0 0 0 1.0000 0 Bb = 1 0 0 0 Cb = 0 0 0.2500 0.5000 Db = 0 从而可以应用数值积分方法求取系统的阶跃响应。下面给出了应用 MATLAB 语言编 写的四阶龙格-库塔数值积分方法的仿真程序，相应的仿真曲线由图 4-4 给出。 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];

图44四阶龙格-库塔数值积分方法的仿真曲线 Cb=[000.250.5 KI=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+h*(K1+2*K2+2*K3+K4)/6 Y=[Y, Cb*X t=[t, t(i)+h]: plot(t,Y) 当由已知系统开环传递函数求取系统开环状态方程时,除了可以利用ss函数直 接求取外,还可以使用状态方程相似变换函数sss进行求取,下面给出上例中应用 相似变换函数ss求取系统开环状态方程的基本用法。这里主要利用矩阵水平翻转 函数 flipper构造出反对角阵作为变换阵,从而得到相应的状态方程 num=2*1, 2: den= conv(conv(1, 0],[1, 3 ).1, 2, 2D) 1= tf(num, den); G2=SS(G1) plr(eye(4)); G=sS2ss(G2, T) 12

127 图 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+ 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); T = fliplr (eye(4)); G = ss2ss(G2,T)

运行结果: 2.00000 0 -0.75000 2.00000 500000 b uI 0 41.00000 y10.5000 0.25000 d 0 从计算结果可以看出,该状态方程十分类似通常定义的系统可控标准型形式。 当然如果在研究系统时,需要状态方程以可控标准型形式给出,则可应用MA∏LAB 编程得到。但不论哪种形式的状态方程,最终所要研究的系统阶跃响应都是一样的 4.4面向结构图的连续系统仿真 4.4.1基于典型环节的系统仿真 由于习惯上往往采用结构图形式对系统进行描述,因此在对系统进行仿真时, 有必要研究如何以结构图为基础进行。在此将介绍以典型环节为基本环节的面向结 构图的系统仿真方法,这种方法的主要优点是 ①便于研究各环节参数对系统的影响 ②可以得到每个环节的动态响应; ③对于存在系统非线性环节较容易处理; ④不但适合单人单出系统,而且对于多入多出系统也适用 该方法只要建立起各个典型环节的以及各环节之间的连接关系,就能方便地求 取任意环节和系统的动态响应。由于直接输入各环节的参数,因此研究参数变化对 12

128 运行结果： 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 = u1 x1 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 编程得到。但不论哪种形式的状态方程，最终所要研究的系统阶跃响应都是一样的。 4．4 面向结构图的连续系统仿真 4．4．1 基于典型环节的系统仿真 由于习惯上往往采用结构图形式对系统进行描述，因此在对系统进行仿真时， 有必要研究如何以结构图为基础进行。在此将介绍以典型环节为基本环节的面向结 构图的系统仿真方法，这种方法的主要优点是： ① 便于研究各环节参数对系统的影响 ② 可以得到每个环节的动态响应； ③ 对于存在系统非线性环节较容易处理； ④ 不但适合单人单出系统，而且对于多入多出系统也适用。 该方法只要建立起各个典型环节的以及各环节之间的连接关系，就能方便地求 取任意环节和系统的动态响应。由于直接输入各环节的参数，因此研究参数变化对