连续系统仿真的方法

第3章 连续系统仿真的方法3.1 数值积分法连续系统数值积分法，就是利用数值积分方法对广微分方程建立离散化形式的数学模型——差分方程，并求其数值解。

可以想象在数学计算机上构造若干个数字积分器，利用这些数字积分器进行积分运算。

在数字计算机上构造数字积分器的方法就是数值积分法，因而数字机的硬件特点决定了这种积分运算必须是离散和串行的。

把被仿真系统表示成一阶微分方程组或状态方程的形式。

一阶向量微分方程及初值为()(),00t Y Y t Y ⎧⎫⎪⎪⎨⎬⎪⎪⎩⎭&Y =F =（3-1）其中，Y 为n 维状态向量，F （t ，Y ）为n 维向量函数。

设方程（3-1）在011,,,,n n t t t t t +=…处的形式上的连续解为()()()()n+1n+1t t n+10t t t =Y t +,(),n Y F t Y dt Y t F t Y dt=+⎰⎰（3-2）设 n =()n Y Y t ，令1n n n Y Y Q +=+（3-3）则有：()1n+1t n Y Y +=也就是说，1(,)n nt n t Q F t Y dt +≈⎰（3-4）如果n Y 准确解()n Y t 为近似值，n Q 是准确积分值的近似值，则式（3-4）就是式（3-2）的近似公式。

换句话说，连续系统的数值解就转化为相邻两个时间点上的数值积分问题。

因此，所谓数值解法，就是寻求初值问题（3-1）的真解在一系列离散点12n t t t <…<…上的近似解12,,,n Y Y Y ……，相邻两个时间离散点的间隔1n n n t t +=-h ，称为计算步距或步长，通常取n =h h 为定值。

可见，数值积分法的主要问题归结为对函数(,)F t y 的数值积分问题，即如何求出该函数定积分的近似解。

为此，首先要把连续变量问题用数值积分方法转化成离散的差分方程的初值问题，然后根据已知的初值条件0y ，逐步地递推计算后续时刻的数值解(1,2,)i y i =…。

所以，解初值问题的数值方法的共同特点是步进式的，采用不同的递推算法，就出现各种不同的数值积分方法。

3.2 替换法基于数值积分的连续系统仿真方法具有成熟、计算精度比较高的优点，但算法公式比较复杂、计算量比较大，通常只有在对速度要求不高的纯数字仿真时使用。

当进行实时仿真或在计算机控制系统中实现数字控制器的算法时，要求计算速度快，以便能在一个采样周期内完成全部计算任务，这就需要一些快速计算方法。

用数值积分方法在数字机上对一个连续系统进行仿真时，实际上已经进行了离散化处理，只不过在离散化过程中每一步都用到连续系统的模型，离散一步计算一步。

那么，能否先对连续的模型进行离散化处理，得到一个“等效”的离散化模型，以后的每一步计算都直接在这个离散化模型基础上进行，而原来的连续数学模型不再参与计算呢？回答是肯定的。

这些结构上比较简单的离散化模型，便于在计算机上求解，不仅用于连续系统数字仿真，而且也可用于数字控制器在计算机上实现。

替换法的基本思想是：对于给定的函数G （s ），设法找到s 域到z 域的的某种映射关系，它将S 域的变量s 映射到z 平面上，由此得到与连续系统传递函数G （s ）相对应的离散传函G （z ）。

进而再根据G （z ）由z 反变换求的系统的时域离散模型——差分方程，据此便可以进行快速求解。

根据z 变换理论，s 域到z 域的最基本的映射关系是Ts Z e =或 1ln s z T=如果按这一映射关系直接代入G （s ），得到的G （z ）是相当复杂的，不便于算法实现，所以往往借助于Z 变换的基本映射关系Ts Z e =或1ln s z T=作一些简化和近似处理。

3.3 离散相似法“离散相似法”——将一个连续系统进行离散化处理,然后求得与它等价的离散模型(差分方程)的方法。

获取离散相似模型的两个途径：（1）对传递函数作离散化处理得离散传递函数——称为“频域离散相似模型”；（2）基于状态方程离散化——称为“时域离散相似模型[3]”；对连续系统进行数字仿真可以先在系统加入虚拟的采样器和保持器，如图3-1所示，保持器G （s ）uyTT图3-1 连续系统离散化结构图附注：图3-1所示系统的采样开关和保持器实际上是不存在的，而是为了将(3-5)式离散化而虚构的。

然后利用Z 变换的方法求出系统的脉冲传递函数，再从脉冲传递函数求出对应于系统G(s)的差分方程。

根据图3-1,有脉冲传递函数 :()()()()()Y z G z Z G s G s h U z ⎡⎤==⎣⎦（3-5）其中Gh (s )是保持器的传递函数。

若选择不同的保持器，则可得不同的G (z )，见表3-1。

表3-1 不同保持器的G (z )假设连续系统的状态方程为:x Ax Bu=+&（3-5）若人为地在系统的输入端及输出端加上采样开关，同时为了使输入信号复原为原来的信号，在输入端还要加一个保持器，如图3-2所示。

零阶保持器x=Ax+BuuyTT图3-2 采样控制系统结构图若对方程(3-5)式两边进行拉普拉斯变换，得：即： ()()(0)()sI A X s X BU s -=+ 以1()sI A --左乘上式的两边可得 ：11()()(0)()()X s sI A X sI A BU s --=-+-（3-6）考虑到状态转移矩阵：11()At t e L sI A --⎡⎤Φ()==-⎢⎥⎣⎦（3-7）故对(3-6)式反变换可得：保持器的传递函数Gh (s ) 脉冲传递函数G (z )零阶：1Tse s --1()z G s Z z s -⎡⎤⎢⎥⎣⎦一阶： 1Tse s--1()(1)2()2z G s Ts Z z Ts ⎡⎤-+⎢⎥⎣⎦三角形：2(1)2Ts Tse e Ts --2(1)()2z G s Z z Ts ⎡⎤-⎢⎥⎣⎦)()()()(s BU s AX X s sX +=-0()()(0)()0A t At t X t e x e BU d τττ-=+⎰（3-8）此为(3-5)式的连续解，由此可推导出系统的离散解。

根据上式，n 及n+1两个相连的采样瞬间，有：()()(0)()0A kT AkT kT X kT e x e BU d τττ-=+⎰ （3-9） [][](1)(1)(1)(1)(0)()0A k T A k T k T X k T e x e BU d τττ+-+++=+⎰（3-10）将(3-10)式减去(3-9)式后乘以AT e ，得：将(3-11)式右边积分进行变量代换，即令：kT tτ=+（3-12） 则得：[]()(1)()()0A T t AT T X k T e X kT e BU kT t dt -+=++⎰ （3-13）但由图3-2可知：若系统采用零阶保持器时，则两个采样点之间输入量可看做常数，即u(nT+t)=u(nT)，这样(3-13)式可写为：[]()0(1)()()()()()()TATA T t X k T e X kT e Bdt U kT G T X kT H T U kT -⎡⎤+=+⎢⎥⎣⎦=+⎰式中 ：()()()0ATG T e A T t T H T e Bdt =-=⎰第6章 计算机仿真实例6.1 连续系统仿真的离散相似法在研究对象的数学模型时,通过模拟研究可以预测这一对象在不同的输入向量的作用下的行为,可为模型的简化提供数据。

通常通过计算机仿真技[][](1)(1)(1)()()A k T k T AT X k T e X kT e BU d kT τττ+-++=+⎰（3-11）术可以估计各种不同的控制系统,在各种干扰作用下的过渡过程,进行方案的分析比较,为选择最好的方案提供依据[7]。

例如,对于一个复杂组分的控制系统,采用数字计算机进行模拟,可以得到各种工况下的控制系统仿真分布曲线,为正确选择仿真方法及路线提供可靠依据,并可以预测控制系统的动态响应效果,所以在自动控制系统的设计、分析和研究中,计算机仿真技术是一有效的手段。

控制系统方框图如图6-1所示，分析k=1,2时的系统的动态响应，（饱和非线性环节斜率为1），)(1.25)(t t R =。

图6-1 控制系统方框图用离散相似法分析计算如下：第一步 引入采样开关的零阶保持器，变成离散控制系统,如图6-2 所示。

图6-2 控制系统方框图第二步 求对象和调节器的状态方程 ，传递函数为)1()(+=s s ns G由控制系统方框图中的传递函数，给出状态方程1s -1s -∑yx2n x1-1u图6-3 系统状态图取：n=1，1011110022xx u x x ⎡⎤⎡⎤-⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦&& 2y x = 即：1010A -⎡⎤=⎢⎥⎣⎦， 10B ⎡⎤=⎢⎥⎣⎦由图中比例限幅调节器的特性可列出：10[()()]10()[()()],10[()()]1010[()()]10R k y k u k R k y k R k y k R k y k ->⎧⎫⎪⎪=--≤-≤⎨⎬⎪⎪--<-⎩⎭其中y （k ）=[0 1]x(k)第三步，离散状态方程)()()()()1(k U T H K X T G k x +=+第四步，求取)(T G 和)(T H ，{}111101()()11(1)ATs G T eL sI A L s s s ---⎧⎫⎡⎤⎪⎪⎢⎥+⎪⎪⎢⎥==-=⎨⎬⎢⎥⎪⎪⎢⎥+⎪⎪⎣⎦⎩⎭=011TT e e --⎡⎤⎢⎥-⎣⎦00010()0111TTTT T At TT ee H T e Bdt dt dt e e----⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎰⎰⎰ =11T T e T e --⎡⎤-⎢⎥-+⎣⎦设 ()s φ=1()sI A -- 称为预解矩阵。

det （sI-A ）：为其的行列式，adj （sI-A ）：为其的伴随矩阵，预解矩阵：()()det()adj sI A s sI A φ-=-若取45.0=T 则：0.637630()0.362371G T ⎡⎤=⎢⎥⎣⎦0.36237()0.08763H T ⎡⎤=⎢⎥⎣⎦ 第五步，由图6-1构成的离散系统，用递推法求解， 已知：0=k 时，1.25)0(=R ，0)0(=y ，0)0(=u ，0)0(=x1.25)0()0()0(=-=+y R E由调节器方程可知： 10)0(=+u可求出1,k =时0.63763000.36237 3.6237()100.36237100.087630.8763x k ⎡⎤⎡⎤⎡⎤⎡⎤=+=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦2k =，时0.637630 3.62370.3623 5.9336()100.3623710.87630.08763 3.0657x k ⎡⎤⎡⎤⎡⎤⎡⎤=+=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦这样一直运算下去，通过编程，计算机仿真，应用离散相似法控制系统计算机仿真如图6-4示，图（a ）1k =，图（b ）2k =051015202530354045102030（a）k=151015202530354045010203040（b）k=2图6-4 仿真结果图小结：1. 由于各个环节的输入量U（I）及输出量Y（I）每一步的数值都可求出，所以这个程序很容易被推广到包含有非线性环节的系统仿真中去。