控制系统数字仿真理论

1、第 7 章控制系统数字仿真理论7.1 弓 I 言仿真主要采用相似性原理。因实际系统是连续的，而计算机系统是离散的（尽管计算机的主频目前可达1 GHz以上，但仍然是断续的），故用计算机进行仿真有两种关键技术：1）建立实际系统的数学模型。2）实际系统的离散方法。系统的离散化方法主要分为两大类，即数值积分方法和直接离散化方法。常用数值积分方法按递推时所需数据步数分为单步法、多步法和预估一校正法。（1）数字仿真的特点连续系统的数学模型一般是微分方程或偏微分方程，因此数字仿真中的主要数值计算工作是微分方程（或偏微分方程）数值解的问题。数字仿真的整个过程是由事先编好的仿真程序 来控制。在大系统实时或超实时

2、仿真中，仿真速度成为一个十分突出的问题。（2）系统仿真技术新动向一个实际的系统可分为连续系统、离散系统、混合系统和定性系统（模糊理论）。而仿真根据其采用的对象可分为计算机仿真、半实物仿真、比例模型仿真和人在回路中仿真。 根据信号的类别可分为数字仿真、模拟仿真、混合仿真。根据仿真时间可分为实时仿真、 超实时仿真（n:1）和欠实时仿真（1:n）。根据应用情况可分为工程系统仿真和非工程系统仿真。还可根据分布情况分为集中式和分布式仿真。系统仿真技术的新动向是：采用分布式、开放式、交互式构架体系，面向对象、网络和 数据库的标准化的应用多媒体和虚拟现实技术进行系统仿真。其发展目标是构成可操作性、 可移植性

3、、交互性强，开放式的仿真体系构架。（3）仿真的可信度仿真的可信度取决于模型的准确性、环境模拟的准确性和干扰处理等3个因素。虚拟现实（virtual reality或灵境，缩写为VR）1989年，美国计算机科学家Jaron Lanier赋VR以现在的含义。虚拟现实综合运用了计算 机图形学、仿真技术、人机接口技术、多媒体技术、传感器技术等，能感知方向、听觉、视 觉、触觉、嗅觉、 味觉，使人有身临其境的感觉。传感器主要有：头盔显示器、数据手套、触觉与力度传感器；跟踪球；空间探针等。1）VR的4个重要特征：multi-se nsory多传感器；immersio n （prese nee）临场感；in

4、teraction交互；aut onomy自主性。2）VR的5个研究内容及关键技术：动态环境建模技术； 实时3D图形生成技术，最好30帧/s以上；立体显示和传感器技术；VR环境的开发平台（VRT，WTK）;系统集成技术，包括信息同步技术、模型标定技术等。7.2 仿真理论要在数字计算机上进行连续系统的仿真，必须先将连续模型变换为离散化的模型，然 后迭代递推出要仿真的变量结果。系统的离散化方法主要分为两大类，即数值积分方法和直接离散化方法。7.2.1数值积分法常用的数值积分方法按递推时所需数据步数分为单步法、多步法和预估 衣正方法3种。(1)单步法属于单步法的主要有欧拉(Euler)法和龙格库塔(

5、Runge-Kutta)法。其中欧拉法最简 单，但由于它有明显的几何意义，可以比较清楚地看出其数值解是如何逼近微分方程精 确解的。1)欧拉法设有一微分方程y(t) =f(t,y(t),且 y(0) =y(7-2-1)若把式(7-2-1)在某一区间(tn,tn 1)上积分，则可得yn i yn二tnn 1f(t, y(t)dt上式右端积分若以一近似公式代之，即tnn1f(t, y(t)dt =hfn其中，h=tn1-tn，即步长。令fn二f(tn，y(tn)，Yn 1二y(tn 1)，Yn二y(tn)，只要h取值比较小，就可以认为：在该步长内的导数近似保持前一时刻tn时的导数值fn。这样用欧拉法

6、离散式(7-2-1)后的递推公式为yn 1=ynhfn(7-2-2)因已知y(O)=y,所以由式(7-2-2)可以求出 力，然后求出y，以此类推。其一般规律即是： 由前点tn上的数值yn就可以求得后一点tn 1上的数值yn 1。这种方法称为单步法。由于它可以直接由微分方程已知的初始值y0作为它递推计算时的初值，而不需其他信息，因此它是一种自启动的算式。下面用一简单例子说明欧拉法的应用及其数值解与精确解的误差。【例7.2.1】 设一微分方程为y y2=0, y(0)=1,试用欧拉法求其数值解。解：因欧拉法递推公式为yn= ynhfn，现y=- y2，所以f (y - y2。若取步长h=0.1，由

7、t=0开始积分，则可得2y1=1+(0.1)( -1 )=0.9y2=0.9+(0.1)-(0.9)2=0.819yn 1 =ynhbjkj2y3=0.819+(0.1)_ (0.819)2=0.752yio=0.482该例的精确解为：y =1/(1 t)。以上结果与精确解比较如下：t00.10.20.31.0精确解y(t)10.909 090 90.833 333 30.769 230 70.5数值解yn10.90.8190.7520.482欧拉法的几何意义为：对微分方程求解的几何意义为求微分函数曲线下面的面积，而欧拉离散法是用采样区间的矩形面积代替其实际的曲线下的面积，然后把所有矩形面积相

8、加得到微分方程的近似解。见图7.2.1所示。欧拉法的代数意义为：微分函数的泰勒级数在tn处展开并保留到一阶。欧拉法可由给定的初值一直递推，而不需其他信息，直到递推出满足精度要求的结yn 1二ynhfn(h2/2)(fnfynfn) (7-2-4)由本例已可看出欧拉法的误差是比较大的,其误差的数量级在10工左右。果为止。其缺点是精度较低，其误差的数量级在10 左右。2)龙格-库塔法为得到精度较高的数值积分方法，龙格和库塔两人先后提出了用函数值f的线性组合来代替f的高阶导数项，则既可避免计算高阶导数，又可提高数值积分的精度。其方 法如下：先将精确解y(t)在y(tntn附近用泰勒级数展开成- h2

9、 h3-h) =y(tn) hy(tn) y(tn) y(tn) 23!(7-2-3)y(tn)y(tn)二fn- fynfn为避免计算fn,fyn等导数项，可以令yn出由以下算式表示:(725)id其中r即阶数，bi是待定系数，ki二f(tn Cjh,yn hljkjhs =0。j=i当r=1时，yn+=yn+hfn即欧拉法。当r=2(7-2-6)即2阶龙格一库塔法。将f (tnC2h,yn- aikih)在(tn，yn)点附近用泰勒级数展开可得f (tn+C2h, yn+aikih)f (tn, ynC2hfn+fynh(7-2-7)将式(7-2-6)和式(7-2-7)代入式(7-2-5)

10、则得yn#=yn+bihki+b2hk2=yn+bihfn+b2h( fn+C2hfn+aifnfynh)(7-2-8)式(7-2-4)与式(7-2-8)右端对应项系数相等，则可得到以下关系式：b| b2= 1b2c2=1/2b2a1=1/2因上述方程组中有4个未知数a1,b1,b2,c2,为求解方程组，可先设定一未知数，常用的有以下几种：设a1=1/2,则c2=1/2,6 = 0,b2=1;设a1=2/3,则c2=2/3, d =1/4,b2=3/4；设a1=1,则c2=1,切=1/2,b2=1/2。相应地，2阶龙格-库法有3个常用的递推公式，在实际应用时可取其中任意一个，即ynynhf(t

11、nh/2,ynhfn/2)yni二yn(h/4)fn3f(tn2h/3,yn2hfn/3)(7-2-9)yni二yn(h/2)fnf (tn- hy hfn)2阶龙格-库塔法的几何意义为：在采样区间内插一个值，然后用插值处的函数值为高的小矩形面积或在插值处分开的2个小矩形的面积代替其实际的曲线下的面积，然后把所有矩形面积相加得到微分方程的近似解，如图7.2.2所示。|_ki= f (tn, yn) = fnk2二f(tnC2h, yn8我山)Otnh tn it图 722 龙格佯塔法几何意义示意图2阶龙格-库塔法的代数意义为：微分函数的泰勒级数在tn处展开并保留到2阶。下面给出3阶和4阶龙格

12、弄塔法的递推公式，其中高且编程实现容易而使用最为广泛。4阶龙格库塔法由于其精度较3阶：yn 1二yn(h/4)(k13k3)(7-2-10)其中：& 二f(tn,yn)k2= f (tnh/3, ynk1h/3) k3=f (tn2h/3, yn2k2h/3)k2=f(tnh/2,ynkih/2) kf (tnh/2,ynk2h/2) k4=f (tnh, ynhk3)对于大部分实际问题，4阶龙格-库塔法已可满足精度要求，它的截断误差正比于h5。龙格-库塔法也可由给定的初值一直递推，而不需其他信息，直到递推出满足精度要求 的结果为止。欧拉法和龙格 -库塔法在递推时只需要前一步的y和f值

13、，故称单步法。在控制系统中，根据实际经验，一般选取采样时间步长为h=1/(5 c)(7-2-12)其中：c为系统开环频率特性的剪切频率。(2)多步法用多步法求解yi时，可能需要y及f(t,y)在tn,tn4,t2,t各时刻的值。1)亚当斯-巴什福思(Adams-Bashforth)显式公式其递推计算公式如下：h2f _ fynyn(h/2)(3fn- fz) hfn)2 h(7-2-13)它是由泰勒级数向前展开式推导得到的。由于yn 1可由yn, fn, fn 4等确定，因此是显式解。但是为了求解一个新的4阶 ：(7-2-11)其中：k1=f(tn,yn)yn 1二yn(h/6)(ki2k22

14、k3k4)y值，需要f的两个值fn及fnj，因此这一递推式不能从t=0自起步。一般常用同阶的龙格一库塔法来启动。2)亚当斯-莫尔顿(Adams-Moulton)隐式公式和梯形积分法Adams-Moult on隐式公式的递推公式如下：yn1= ynhfn 1(7-2-14)它是由向后展开的泰勒级数公式推导得到的。由于yn 1式中包含fn 1，而fn 1计算时又要用到Yn 1，因此，要解出Yn 1就要用迭 代法。其步骤是先估算一个yn 1，计算fn 1，而后用上式求得Yn 1的新估值，重复迭代，直到前后两次yn1值之间的误差在要求范围内为止。由于进行多次迭代运算，因此解的 精度较高，但费时多。梯形

15、积分法的递推公式为yn 1二yn(h/2)( fnfn 1)(7-2-15)这种方法的几何意义比较清楚，它是以梯形面积(h/2)(fn fn 1)来近似原f(t,y)在tn到tn 1之间曲线下的面积。由于积分公式中需有fn -1项，因此与亚当斯隐式公式类似，需用迭代运算方法来求解yn 1值。预估我正法隐式公式一般精度较高，但因其右端包含未知项，所以需要先用另一显式公式估计一个初值，然后再用隐式公式进行迭代运算(或者说，进行校正)，这种方法就称为预估 -校正法。应用此法时需注意：显式和隐式公式的阶数要一致。常用的预估-校正法有：1)采用欧拉法“预估”，梯形法“校正”。2)用2阶亚当斯显式和隐式公

16、式组成预估-校正法。722直接离散化法直接离散化方法有Z变换法、带有零阶保持器的Z变换法、差分反演法及双线性变换法等。数值积分方法把微分方程离散化成不同的迭代算式，其缺点是由于迭代算式 中的系数每一步都要重新计算，因此一般计算量比较大，但适于非线性系统的离散化。而直接用离散化模型代替连续系统数学模型的方法，实质上，就是以常系数差分方程近 似“等效”原来的常系数微分方程。由于差分方程可以直接用迭代方法在计算机上求解，因此非常方便。连续系统离散化的含义是：假设有一连续系统，其输入为u(t)，输出为y(t)。现若用一周期为h的采样开关将输入、输出分别离散化， 要求输出y*( t)在采样时刻的值等于原

17、输出y(t)在同一时刻的值。(1) Z变换法Z变换是现代控制理论中所用到的一种变换，它是把脉冲序列f*(t)的拉氏变换式F*(s)中的ehs换成Z而后得到的F(Z)，我们称之为f*(t)的Z变换。例如，f*(t)是函数f(t)经采样后的脉冲序列，即f*(t)二厂f (kh)、(t - kh)f*(t)的拉氏变换为*-_khsF (s)二、f(kh)e令ehs=z，则匚kF(Z)二 f(kh)Z上k=0,1,2, ,(7-2-16)k=0观察式(7-2-16)可知，ZJ在时域内相当于、；(t - h)，而在采样时刻t = kh则相当于、：(k -1)h)，所以任何一个函数乘以ZJ，则在时域内相当

18、于将该函数提前一个采样周期。例如ZU在时域内相当于U(k -1)或Uk注意这里的k-1为(k_1)h的简写。(2)带有零阶保持器的Z变换在数字采样系统中，为了使系统在非采样周期也能保持稳定，常常在系统的前面增加一个保持器，以便使离散型控制信号转换为连续信号。零阶保持器是一个按常数关系外延的装置，也就是说，零阶保持器使信号在一周期(h)内保持不变，即使ih时刻的输出值保持到(i+1)h时刻。其传函为(1 -esh)/s。设被控对象为D(s)，则(7-2-17)【例7.2.2】 已知D(s)=a/(s+a)，求D(s)带零阶保持器的Z变换。解：D(Z)=Z(1 -eh)/s)D(s)=-sh44e

19、ash11)=Z(1 e ) (二) =s s + a(1 _Z J1-1-) =1匸 -八1_Z1 -ehZ1 -ehZah 1ah 1、D(Z)Z(1 -eh)/s)D(s)Z(!(1 e )Z /(1 e Z )递推公式：yn=e3yn4(1-eh)Xnv(设D(s) =y(s)/x(s)。(3)差分反演法换法求出相应的数字控制器的脉冲传递函数D(Z)。设h=0.4 s。设U(s)/E(s)=D(s)，用一阶差分近似值来代替所有的导数，用d/dt代替s后再换成差分表达式：dE/dtE(k)_E(k_1)/h;dU/dtU(k) _U(k _1)/h其中：E(k), U(k)分别表示其kh

20、时刻的值。因E(k -1)=Z二E(k),故可以求得s=(1 ZJ)/h。于是，差分反演法变换公式为D(Z) =D(s)|卩丄s -h(7-2-18)差分反演变换的特点是：使用比较容易，且不需要将传递函数进行因式分解；D(s)稳定则D(Z)也稳定。【例7.2.3】已知传递函数D(s)=U(s)/e(s)=(1+25 s)/(1+62.5 s), 且h=1 s，试用差分反演法求出数字脉冲传递函数并写出其输出递推公式。解：1+25(1Z)/h26 25Z0.400.39ZD(Z)J11 +62.5(1 -Z)/h63.5-62.5Z10.98Z其输出递推公式为U (k) =0.98U(k -1)0

21、.40e(k) -0.39e(k -1)1通过以上两例可以看出，在求Z变换时分子分母都要化成跟Z有关的多项式，因为ZJ在写递推公式时有实际的物理意义。如果分子分母都通分成含有Z的多项式，则在写递推公式时还需重新变成含有z，的多项式形式，多走了一步，较麻烦。(4)双线性变换法(Tustin变换法)根据Z变换的定义：z=ehs=ehs/2/es/2将此展成Taylor级数，并取前两项：hs/2hs/2e=1+ hs/2 ;e=1 -hs/2于是可得f1 -Z11+Z双线性变换公式为D(Z)二D(s)|2s二h(7-2-19)【例7.2.4】已知某连续控制器的传递函数D(s) =(s 0.5)/(s

22、 - 1)2,试用双线性变7.3 控制系统建模与数字仿真过程（1）控制系统仿真方法选择策略在上节介绍的两类控制系统仿真方法中，数值积分类仿真方法相对于直接离散化方法具有适用于非线性系统的优点，而计算量大则是它的缺点；直接离散化方法具有计算量小的优点，只适用于线性系统是它的缺点。这两类仿真方法的优缺点具有互补性。在一个标准的反馈控制系统中，可以利用这两类方法的优缺点，针对被控对象和控制器本身的特性，分别选择不同的仿真方法进行离散化。一般情况下，在控制系统中，要求被控对象的数学模型具有较高的精度，所以采用非线性数学模型，其离散化方法选择数值积分类方法，因为龙格库塔4阶或5阶方法具有较高的精度，所以作为首选；控制系统中的控制器主要是自己设计的，要求其结构简单，易于实现，一般是线性的，选择直接离散化方法是择其计算量小的优点，因为双线性变换方法具有较高的精度，所以作为首选。这样，对于控制系统的被控对象和控制器选择了两类不同的仿真方 法进行离散化，根据各自的特性充分利用了两类仿真方法的优点，使两类仿真方法共存于一个控制系统仿真中。（2