系统仿真技术Chapter2经典的连续系统仿真建模方法学

1、控制系统仿真技术控制系统仿真技术盛盛 立立中国石油大学自动化系chapter 2 第2章经典的连续系统仿真建模方法学 对下面的控制系统描述，需要放在计算机上求解对下面的控制系统描述，需要放在计算机上求解p常微分方程常微分方程 p传递函数传递函数p状态空间描述状态空间描述 11201112112nnnnnnnnnnnd ydydyduduaaaa yccc udtdtdtdtdt121211011( )( )( )nnnnnnnnc sc scscy sg su sa sa sasabuaxxcxy方法一：方法一：ode23,ode45可解一阶微分方程组，可解一阶微分方程组， 状态空间描述是一阶

2、微分方程组状态空间描述是一阶微分方程组常微分方程，传递函数常微分方程，传递函数 状态空间表达式状态空间表达式求解？求解？ ode23,ode45可解一阶微分方程组，可解一阶微分方程组，原理是什么？原理是什么？对下面的控制系统描述，需要放在计算机上求解对下面的控制系统描述，需要放在计算机上求解p常微分方程常微分方程 p传递函数传递函数p状态空间描述状态空间描述 buaxxcxy方法一：方法一：ode23,ode45可解一阶微分方程组，可解一阶微分方程组， 状态空间描述是一阶微分方程组状态空间描述是一阶微分方程组常微分方程，传递函数常微分方程，传递函数 状态空间表达式状态空间表达式原理原理: 一阶

3、微分方程（线性，非线性）一阶微分方程（线性，非线性）数值求解方法数值求解方法11201112112nnnnnnnnnnnd ydydyduduaaaa yccc udtdtdtdtdt121211011( )( )( )nnnnnnnnc sc scscy sg su sa sa sasa 数数值值求求解解方方法法 欧拉法欧拉法 梯形法梯形法 龙格库塔法龙格库塔法rk2rk42.1 离散化原理及要求离散化原理及要求 p问题：数字计算机在数值及时间上的问题：数字计算机在数值及时间上的离散性离散性-被仿真系统数值及时间上的被仿真系统数值及时间上的连续性连续性?p连续系统的仿真，从本质上：对原连续系

4、统从时连续系统的仿真，从本质上：对原连续系统从时间、数值两个方面间、数值两个方面对原系统进行离散化对原系统进行离散化并选择合并选择合适的数值计算方法来适的数值计算方法来近似近似积分运算积分运算 p 离散模型离散模型原连续模型？原连续模型？ 相似原理相似原理 p设系统模型为：设系统模型为： ，其中，其中u(t)为输为输入变量，入变量，y(t)为系统变量；令仿真时间间隔为为系统变量；令仿真时间间隔为h，离散化后的输入变量为，离散化后的输入变量为 ，系统变量，系统变量为为 ，其中，其中 表示表示t=nh。p如果如果 , 且即即 , （对所有（对所有n=0,1,2,） ,则可认为两模型等价。则可认为两

5、模型等价。 ),(tuyfy )( ntu)( ntynt)()( nntutu)()( nntyty0)()( )(nnnututute0)()( )(nnnytytyteu(t)h y(t) -+图图2.1 2.1 相似原理相似原理原连续模型 ),(tuyfy 仿真模型 ), , (ntuyfy )( nty0)(nyte)( ntu相似原理相似原理 对仿真建模方法三个基本要求对仿真建模方法三个基本要求 p（1）稳定性稳定性：若原连续系统是稳定的，则离散化后得到：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。的仿真模型也应是稳定的。 p（2）准确性准确性：有不同的准确性评价准

6、则，最基本的准则：有不同的准确性评价准则，最基本的准则是：是： 绝对误差准则：绝对误差准则： 相对误差准则：相对误差准则： 其中其中 规定精度的误差量。规定精度的误差量。 )()( )(nnnytytyte)( )()( )(nnnnytytytytep（3）快速性快速性：若第：若第k步计算对应的系统时间间隔步计算对应的系统时间间隔为为 计算机由计算需要的时间为计算机由计算需要的时间为 ，若，若 tn=hn 称为称为实时仿真实时仿真，tn hn称为称为超实时仿超实时仿真真, tn hn 称为称为亚实时仿真。亚实时仿真。1,nnnhttnt,nnth,nnth,nnth 系统仿真中最常用、最基本

7、的求解常微分方程数值解系统仿真中最常用、最基本的求解常微分方程数值解的方法主要是数值积分法。的方法主要是数值积分法。 设系统常微分方程为：设系统常微分方程为： （2-1） 为包含有时间为包含有时间t和函数和函数y的表达式，的表达式，y0为函数为函数y在初在初始时刻始时刻t0时的对应初值。我们将求解方程（时的对应初值。我们将求解方程（2-1）中函数）中函数 的问题称为常微分方程数值求解问题。的问题称为常微分方程数值求解问题。00)(),(ytyytfdtdy),( ytf)(ty2.2 数值积分法数值积分法 1欧拉公式的推导 将（2-1）式在小区间上进行积分可得： 1),(1kkttkkdtyt

8、fyy),(),(1kkttythfdtytfkk其几何意义是把 ),( ytf在,1kktt 区间内的曲边面积用矩形面积近似代替，如图2-1所示。 2.2.1 欧拉法欧拉法euler00)(),(ytyytfdtdy1),(1kkttkkdtytfyytf(t,y)0fktktk+1h图2-1 欧拉法数值积分 欧拉法欧拉法euler当h很小时，可以认为造成的误差是允许的。所以有：),(1kkkkythfyy称之为欧拉公式。 p截断误差正比于截断误差正比于 h2欧拉欧拉euler公式：公式：),(1kkkkythfyy00)(ytyp截断误差正比于截断误差正比于 h22 . 欧拉法具备以下特点

9、：欧拉法具备以下特点：（1）欧拉法实际上是采用折线代替了实际曲线，也称之为）欧拉法实际上是采用折线代替了实际曲线，也称之为折线法。折线法。（2）欧拉法计算简单，容易实现。由前一点值仅一步递推）欧拉法计算简单，容易实现。由前一点值仅一步递推就可以求出后一点值，所以称为单步法。就可以求出后一点值，所以称为单步法。（3）欧拉法计算只要给定初始值，即可开始进行递推运算，）欧拉法计算只要给定初始值，即可开始进行递推运算，不需要其它信息，因此它属于自启动模式。不需要其它信息，因此它属于自启动模式。（4）欧拉法是一种近似的处理，存在计算误差，所以系统）欧拉法是一种近似的处理，存在计算误差，所以系统的计算精度

10、较低。的计算精度较低。 欧拉法欧拉法euler 数数值值求求解解方方法法 欧拉法欧拉法 梯形法梯形法 龙格库塔法龙格库塔法rk2rk4p截断误差正比于截断误差正比于 h21梯形公式 为了弥补欧拉法计算精度较低的不足，可以采用梯形面积公式来代替曲线下的定积分计算，如图2-2所示。 依然对式（2-1）进行求解，采用梯形法作相应近似处理之后，其输出为： ),(),(2111kkkkkkytfytfhyy 称为梯形积分公式 。2.2.2 梯形法梯形法tf(t,y)0fktktk+1hfk+1图2-2 梯形法数值积分 梯形法梯形法 从中可以看到，在计算 时，其右端函数中也含有 ，这种公式称为隐式公式，不

11、能靠自身解决，需要采用迭代方法来启动，称之为多步法。可以先采用欧拉公式进行预报，再利用梯形公式进行校正。即梯形法的预报校正公式 ： 1ky1ky),(),(21),(1)0(111)0(kkkkkkkkkkytfytfhyyythfyy梯形法梯形法2. 梯形法具备以下特点：梯形法具备以下特点：（1）采用梯形代替欧拉法的矩形来计算积分面积，其计）采用梯形代替欧拉法的矩形来计算积分面积，其计算精度要高于欧拉法。算精度要高于欧拉法。（2）采用预报）采用预报校正公式，每求一个校正公式，每求一个 ，计算量要比，计算量要比欧拉法多一倍。因此计算速度较慢。欧拉法多一倍。因此计算速度较慢。（3）梯形公式中的右

12、端函数含有未知数，不能直接计算）梯形公式中的右端函数含有未知数，不能直接计算左端的变量值，这是一种隐式处理，要利用迭代法求解。左端的变量值，这是一种隐式处理，要利用迭代法求解。即梯形法不能自启动，要靠多步法来实现计算。即梯形法不能自启动，要靠多步法来实现计算。ky梯形法梯形法 数数值值求求解解方方法法 欧拉法欧拉法 梯形法梯形法 龙格库塔法龙格库塔法rk2rk4截断误差正比于截断误差正比于 ，记为，记为 h2)(2ho截断误差正比于截断误差正比于 ，记为，记为 3h)(3ho欧拉欧拉euler公式：公式：),(1kkkkythfyy00)(ytyp截断误差正比于截断误差正比于 h2),(),(

13、21),(1)0(111)0(kkkkkkkkkkytfytfhyyythfyy梯形法公式梯形法公式p截断误差正比于截断误差正比于 3h2.2.3 龙格库塔法龙格库塔法 1 龙格龙格-库塔法基本原理库塔法基本原理 对对 的数值求解：称作的数值求解：称作“右端函数右端函数”计算问题。计算问题。 将将 在在 附近展开附近展开taylor级数，只保留级数，只保留 项，则有：项，则有：1),()()(1nnttnndtytftyty)(nntyy 1),(qnnttndtytfnnnnyytyq)(11nqt0h2200010)(21),(htfdtdyyfhytfyyt若令：若令：则有则有 00)(

14、),(ytyytfdtdy1y 假设这个解可以写成如下形式：假设这个解可以写成如下形式： 其中其中 对对 式右端的函数展成式右端的函数展成taylor级数，保留级数，保留h项，项，可得：可得： 代入，则有：代入，则有： yya ka kh101122(),(001ytfk kf tb hyb k h201021()，k2hyfkbtfbytfkt0)(),(121002)(),(),(012100200101hyfkbtfbytfhaythfayyt200010)(21),(htfdtdyyfhytfyyt200010)(21),(htfdtdyyfhytfyyt龙格龙格-库塔法基本原理（续）

15、库塔法基本原理（续） 进行比较，可得：进行比较，可得： 四个未知数四个未知数 但只有三个方程，但只有三个方程，因此有无穷多个解。因此有无穷多个解。 若限定若限定 ，则，则 计算公式：计算公式： 其中其中 aaa ba b122 12 211 21 2+=,/ ,/ ，2121bbaaaa12aabb1212121，)(22101kkhyy)(),(1002001hkyhtfkytfk，龙格龙格-库塔法基本原理（续）库塔法基本原理（续）若写成一般递推形式，即为：若写成一般递推形式，即为： 其中其中p截断误差正比于截断误差正比于h3，称为二阶龙格，称为二阶龙格-库塔法（简称库塔法（简称rk-2）。

16、）。 )(2)(2111kkhyytynnn),(),(121hkyhtfkytfknnnn，二阶龙格二阶龙格-库塔公式库塔公式)(2)(2111kkhyytynnn),(),(121hkyhtfkytfknnnn四阶龙格库塔公式 ：),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkkkkkk四阶龙格四阶龙格库塔（库塔（rungekutta）法）法p截断误差正比于截断误差正比于 5h（1）为单步法，并且可自启动。）为单步法，并且可自启动。（2）改变仿真步长比较方便，可根据精度要求而定。）改变仿真步长比较方

17、便，可根据精度要求而定。（3）仿真计算量与仿真步长）仿真计算量与仿真步长h的大小密切相关，的大小密切相关，h值越值越小计算精度越高，但所需仿真时间也就越长。小计算精度越高，但所需仿真时间也就越长。（4）用泰勒级数展开龙格库塔法计算公式时，只取）用泰勒级数展开龙格库塔法计算公式时，只取h的一次项，即为欧拉法计算公式；若取到的一次项，即为欧拉法计算公式；若取到h2项，则为项，则为二阶龙格库塔法计算公式；若取到二阶龙格库塔法计算公式；若取到h4项，则为四阶龙项，则为四阶龙格库塔法计算公式。格库塔法计算公式。 龙格库塔法特点龙格库塔法特点 【例2.1】 已知一阶系统的微分方程为： ，初始条件 ，取仿真

18、步长h=0.1，分别用欧拉法、梯形法和龙格库塔法计算该系统仿真第一步的值。102ydtdy1)(00 yty解：原方程可变为: ),(210kkytfydtdy即 1210),(0yyytfkkk2.2.4 数值积分公式应用数值积分公式应用 （1）用欧拉法计算 根据欧拉公式，将函数表达式及其初始值代入后，可得该系统仿真第一步的值：8 .1)1210(1 .01),(0001ythfyy数值积分公式应用数值积分公式应用 ),(1kkkkythfyy（2）用梯形法计算： 根据预报校正公式，将函数表达式及其初始值代入后，可得仿真第一步的值。 用预报公式求起始值：8 . 1) 1210(1 . 01)

19、,(0001)0(ythfyy数值积分公式应用数值积分公式应用 ),(),(21),(1)0(111)0(kkkkkkkkkkytfytfhyyythfyy再用校正公式得到系统仿真第一步的值：72. 1)8 . 1210() 1210(1 . 0211),(),(211)0(10001ytfytfhyy数值积分公式应用数值积分公式应用 二阶龙格二阶龙格-库塔公式库塔公式)(2)(2111kkhyytynnn),(),(121hkyhtfkytfknnnn（3）用二阶龙格库塔法计算 根据公式先计算出两个系数，再计算仿真第一步的值： 812100),(001yytfk4 . 6)81 . 01 (

20、210)(210),(101002hkyhkyhtfk数值积分公式应用数值积分公式应用 则系统仿真第一步的值为：72. 1)4 . 68(1 . 0211)(22101kkhyy数值积分公式应用数值积分公式应用 （4）用四阶龙格库塔公式计算根据公式先计算出4个系数，再计算仿真第一步的值： 81210210),(0001yytfk2 . 7)81 . 0211 (210)2(210)2,2(101002khykhyhtfk数值积分公式应用数值积分公式应用 ),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkk

21、kkkk四阶龙格四阶龙格库塔（库塔（rungekutta）法）法5h28. 7)2 . 71 . 0211 (210)2(210)2,2(202003khykhyhtfk544. 6)28. 71 . 01 (210)(210),(303004hkyhkyhtfk数值积分公式应用数值积分公式应用 则系统仿真第一步的值为：725067. 1)544. 628. 722 . 728(1 . 0611)22(6432101kkkkhyy数值积分公式应用数值积分公式应用 从上述结果可以看出从上述结果可以看出: 对于同一个系统进行仿真计算时，其值的精度是随对于同一个系统进行仿真计算时，其值的精度是随着数

22、值积分公式的变化而改变的，其中欧拉法计算精着数值积分公式的变化而改变的，其中欧拉法计算精度最低，其次为梯形法和二阶龙格度最低，其次为梯形法和二阶龙格库塔法，四阶龙库塔法，四阶龙格格库塔法计算精度最高。库塔法计算精度最高。数值积分公式应用数值积分公式应用 数值积分公式在状态方程中应用数值积分公式在状态方程中应用 )22(6),()2,2()2,2()(),(432113423121kkkkhxxhkxhtfkkhxhtfkkhxhtfktbuaxxxtfkkkkkkkkkkkkk11kkcxycxybaxtxfxu),(2.3.1 仿真精度与系统稳定性仿真精度与系统稳定性1. 仿真过程的误差（1

23、）初始误差:现场采集数据不一定很准，会造成仿真过程中产生误差，称为初始误差。应对现场数据进行准确的检测，也可多次采集，以其平均值作为参考初始数据。（2）舍入误差:由于不同档次的计算机其计算结果的有效值不一致，导致仿真过程出现舍入误差。 应选择挡次高的计算机，其字长越长，仿真数值结果尾数的舍入误差就越小。（3）截断误差:仿真步距确定后，数值积分公式的阶次将导致系统仿真时产生截断误差，阶次越高，截断误差越小。仿真时多采用四阶龙格库塔法，其截断误差较小。2.3 数值积分法性能分析数值积分法性能分析 2. 仿真过程的稳定性仿真过程的稳定性 计算结果对系统仿真的计算误差反应不敏感，称之为算计算结果对系统

24、仿真的计算误差反应不敏感，称之为算法稳定，否则称算法不稳定。对于不稳定的算法，误差会不法稳定，否则称算法不稳定。对于不稳定的算法，误差会不断积累，最终可能导致仿真计算达不到系统要求而失败。断积累，最终可能导致仿真计算达不到系统要求而失败。（1）系统的稳定性与仿真步长的关系）系统的稳定性与仿真步长的关系 一个数值解是否稳定，取决于该系统微分方程的特征根一个数值解是否稳定，取决于该系统微分方程的特征根是否满足稳定性要求，而不同的数值积分公式具有不同的稳是否满足稳定性要求，而不同的数值积分公式具有不同的稳定区域，在仿真时要保证稳定就要合理选择仿真步长，使微定区域，在仿真时要保证稳定就要合理选择仿真步

25、长，使微分方程的解处于稳定区域之中。分方程的解处于稳定区域之中。数值积分法性能分析数值积分法性能分析 （2）积分步长的选择）积分步长的选择 由于积分步长直接与系统的仿真精度和稳定性密由于积分步长直接与系统的仿真精度和稳定性密切相关，所以应合理地选择积分步长切相关，所以应合理地选择积分步长h的值。的值。 通常遵循两个原则：通常遵循两个原则：使仿真系统的算法稳定。使仿真系统的算法稳定。使仿真系统具备一定的计算精度。使仿真系统具备一定的计算精度。 一般掌握的原则是：在保证计算稳定性及计算精一般掌握的原则是：在保证计算稳定性及计算精度的要求下，尽可能选较大的仿真步长。度的要求下，尽可能选较大的仿真步长

26、。数值积分法性能分析数值积分法性能分析 由于工程系统的仿真处理采用四阶龙格由于工程系统的仿真处理采用四阶龙格库塔法居多，库塔法居多，所以选择仿真积分步长可参考以下公式：所以选择仿真积分步长可参考以下公式： 时域内：时域内： ；其中；其中ts 为系统过渡过程调节时间为系统过渡过程调节时间 频域内：频域内： ；其中；其中 为系统的开环截止频率为系统的开环截止频率40sth ch51c数值积分法性能分析数值积分法性能分析 3.速度与精度速度与精度 p四阶方法的四阶方法的h可以比二阶方法的可以比二阶方法的h大大10倍，每步计倍，每步计算量仅比二阶方法大一倍算量仅比二阶方法大一倍p 高于四阶的方法由于每

27、步计算量将增加较多，而高于四阶的方法由于每步计算量将增加较多，而精度提高不快。精度提高不快。 数值积分法性能分析数值积分法性能分析 仿真步长与稳定性关系仿真步长与稳定性关系2.4 稳定性分析稳定性分析 p仿真方法选择的基本要求：仿真方法选择的基本要求：仿真计算不改变原系统仿真计算不改变原系统的绝对稳定性。的绝对稳定性。p原系统是稳定的。观察欧拉法仿真递推公式原系统是稳定的。观察欧拉法仿真递推公式p故有故有 (i) yn (n=0,1,2,)为它的一个仿真解，为它的一个仿真解， iydtdy,/re 0),(1nnnnythfyynnnyhyy1稳定性分析（续）稳定性分析（续）设设 为其准确解，即为其准确解，即 （ii）用用(ii)式减去式减去(i)式，可得：式，可得：即即特征方程为特征方程为nny)(