1、系统仿真技术第第2章章 经典的连续系统仿真建模方法学经典的连续系统仿真建模方法学 xx合肥工业大学机械与汽车工程学院1感谢你的观看2019年8月32.1 离散化原理及要求离散化原理及要求 v问题：数字计算机在数值及时间上的问题：数字计算机在数值及时间上的离散性离散性-被仿真系统数值及时间上的被仿真系统数值及时间上的连续性连续性?v连续系统的仿真，从本质上：对原连续系统连续系统的仿真，从本质上：对原连续系统从时间、数值两个方面从时间、数值两个方面进行离散化进行离散化并选择合并选择合适的数值计算方法来适的数值计算方法来近似近似积分运算积分运算 v 离散模型离散模型原连续模型？原连续模型？ 2感谢你
2、的观看2019年8月3相似原理相似原理 设系统模型为：设系统模型为： ，其中，其中u(t)为输为输入变量，入变量，y(t)为系统变量；令仿真时间间隔为为系统变量；令仿真时间间隔为h，离散化后的输入变量为，离散化后的输入变量为 ，系统变量，系统变量为为 ，其中，其中 表示表示t=nh。 如果如果 , 且即即 , （对所有（对所有n=0,1,2,） 则可认为两模型等价。则可认为两模型等价。 ),(tuyfy )( ntu)( ntynt)()( nntutu)()( nntyty0)()( )(nnnututute0)()( )(nnnytytyte3感谢你的观看2019年8月3u(t)h y(t
3、) -+图图2.1 2.1 相似原理相似原理原连续模型 ),(tuyfy 仿真模型 ), , (ntuyfy )( nty0)(nyte)( ntu4感谢你的观看2019年8月3对仿真建模方法三个基本要求对仿真建模方法三个基本要求 （1）稳定性稳定性：若原连续系统是稳定的，则离散化：若原连续系统是稳定的，则离散化后后得到的仿真模型也应是稳定的。得到的仿真模型也应是稳定的。 （2）准确性准确性：有不同的准确性评价准则，最基本：有不同的准确性评价准则，最基本的的准则是：准则是： v绝对误差准则：绝对误差准则： v相对误差准则：相对误差准则： 其中其中 规定精度的误差量。规定精度的误差量。 )()(
4、 )(nnnytytyte)( )()( )(nnnnytytytyte5感谢你的观看2019年8月3对仿真建模方法三个基本要求对仿真建模方法三个基本要求(续续) （3）快速性快速性：若第：若第n步计算对应的系统时间间隔步计算对应的系统时间间隔为为 计算机由计算需要的时间为计算机由计算需要的时间为Tn，若，若 Tn=hn 称为称为实时仿真实时仿真，Tn hn称为称为超实时仿真超实时仿真 Tn hn 称为称为亚实时仿真亚实时仿真，对应于，对应于离线仿真离线仿真 ,1nnntth6感谢你的观看2019年8月3数值积分算法数值积分算法 对对 ，已知系统变量，已知系统变量y的初始条件的初始条件要求要求
5、y随时间变化的过程随时间变化的过程初值问题初值问题 计算过程：由初始点计算过程：由初始点 的的 欧拉法欧拉法对任意时刻对任意时刻tn+1 截断误差正比于截断误差正比于 ),(tuyfy 00)(yty00)(yty)(00ytf，ttdtytfyty0),()(0yy tytf ty11000( )()，)()()(111nnnnnnnytfttytyy，h2f(t,y)f(t0,yo)ttt0t17感谢你的观看2019年8月3数值积分算法（续）数值积分算法（续） 梯形法梯形法: 是隐函数形式。预报是隐函数形式。预报欧拉法估计初值，校欧拉法估计初值，校正正用梯形法校正：用梯形法校正：v 校正公
6、式校正公式 v 预报公式预报公式 反复迭代，直到满足反复迭代，直到满足经典的数值积分法分为两类：单步法与多步法经典的数值积分法分为两类：单步法与多步法),(),(21)(1111nnnnnnnytfytfhytyy),(),(21)(11)1(1innnnninytfytfhyy),()(1nnninytfhyyininyy1118感谢你的观看2019年8月32.2 龙格库塔法龙格库塔法 2.2.1龙格龙格-库塔法基本原理库塔法基本原理 对对 若令：若令： 则有则有 的数值求解：称作的数值求解：称作“右端函数右端函数”计算问计算问题。题。在在 附近展开泰勒级数，只保留附近展开泰勒级数，只保留
7、项，则有：项，则有：1),()()(1nnttnndtytftyty)(nntyy 1),(QnnttndtytfnnnnyytyQ)(11nQt0h2200010)(21),(htfdtdyyfhytfyyt9感谢你的观看2019年8月3龙格龙格-库塔法基本原理（续）库塔法基本原理（续） 假设这个解可以写成如下形式：假设这个解可以写成如下形式： 其中其中 对对 式右端的函数展成泰勒级数，保留式右端的函数展成泰勒级数，保留h项，项，可得：可得：代入，则有：代入，则有： yya ka kh101122(),(001ytfk kf tb hyb k h201021()，k2hyfkbtfbytfk
8、t0)(),(121002)(),(),(012100200101hyfkbtfbytfhaythfayyt10感谢你的观看2019年8月3龙格龙格-库塔法基本原理（续）库塔法基本原理（续） 将将(2)式与式与(1)式进行比较，可得：式进行比较，可得：四个未知数四个未知数 但只有三个方程，但只有三个方程，因此有无穷多个解。因此有无穷多个解。 若限定若限定 ，则，则计算公式：计算公式：其中其中 aaa ba b122 12 211 21 2+=,/ ,/ ，2121bbaaaa12aabb1212121，)(22101kkhyy)(),(1002001hkyhtfkytfk，11感谢你的观看20
9、19年8月3龙格龙格-库塔法基本原理（续）库塔法基本原理（续） 若写成一般递推形式，即为：若写成一般递推形式，即为：其中其中 （1）截断误差正比于）截断误差正比于h3，称为二阶龙格，称为二阶龙格-库塔法库塔法（简称（简称RK-2）。）。 （2）截断误差正比于）截断误差正比于h5的四阶龙格的四阶龙格-库塔法（简称库塔法（简称RK-4）公式：）公式：其中：其中： )(2)(2111kkhyytynnn),(),(121hkyhtfkytfknnnn，)22(6)(432111kkkkhyytynnn),(1nnytfk )22(12khyhtfknn，)22(23khyhtfknn，)(34hky
10、htfknn，12感谢你的观看2019年8月32.2.2龙格龙格-库塔法的特点库塔法的特点 1.形式多样性形式多样性 例：例： 非唯一解，可以得到许多非唯一解，可以得到许多种龙格种龙格-库塔公式：库塔公式： (中点公式中点公式) 其中其中 各种龙格各种龙格-库塔法可以写成如下一般形式：库塔法可以写成如下一般形式：其中其中 2121bbaa，hkyynn21),(1nnytfk )22(12khyhtfknn，siiinnkChyy11)(11ijjijninikbhyhatfk，si， 2113感谢你的观看2019年8月3龙格龙格-库塔法的特点库塔法的特点(续）续） 式中各系数满足以下关系式中
11、各系数满足以下关系 s称为级数，表示每步计算右端函数称为级数，表示每步计算右端函数f的最少次数。的最少次数。可以证明，可以证明，1阶公式至少要计算一次，阶公式至少要计算一次，2阶公式阶公式 ;.；4阶公式阶公式 ；依此类推。有时为了某种特殊；依此类推。有时为了某种特殊需要，可以选择需要，可以选择 的计算公式。的计算公式。 aabisCiijjiiis11110231， ，smin 2smin 4minss 14感谢你的观看2019年8月3龙格龙格-库塔法的特点库塔法的特点(续）续）2.单步法单步法 在计算在计算 时只用到时只用到 ，而不直接用，而不直接用 等等项。优点：存储量减小，可以自启动项
12、。优点：存储量减小，可以自启动.3.可变步长可变步长 步长步长h在整个计算中并不要求固定，可以根据精度在整个计算中并不要求固定，可以根据精度要求改变要求改变 ，但是在一步中，为计算若干个系数但是在一步中，为计算若干个系数 ，则则必须用同一个步长必须用同一个步长h。 1nyny21nnyy，ik15感谢你的观看2019年8月3龙格龙格-库塔法的特点库塔法的特点(续）续）4.速度与精度速度与精度 四阶方法的四阶方法的h可以比二阶方法的可以比二阶方法的h大大10倍，每步计倍，每步计算量仅比二阶方法大一倍，高于四阶的方法由于每步算量仅比二阶方法大一倍，高于四阶的方法由于每步计算量将增加较多，而精度提高
13、不快。计算量将增加较多，而精度提高不快。 16感谢你的观看2019年8月32.2.3 实时龙格库塔法实时龙格库塔法 实时仿真实时仿真：要求仿真模型的运行速度往往与实际：要求仿真模型的运行速度往往与实际系统运行的速度保持一致。系统运行的速度保持一致。 一般的数值积分法一般的数值积分法难以满足难以满足实时仿真的要求，这实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速度较不仅仅是因为由这些方法所得到的模型的执行速度较慢，而且这些方法的机理不符合实时仿真的特点。慢，而且这些方法的机理不符合实时仿真的特点。 考虑系统考虑系统 )(tuyfdtdy，17感谢你的观看2019年8月3实时龙格库塔
14、法实时龙格库塔法(续）续） RK-2公式如下：公式如下： 一个计算步内分两子步：一个计算步内分两子步： vtn时刻：利用当前的时刻：利用当前的un，yn计算计算k1-计算一次右端函数计算一次右端函数f需需 。vtn+h/2时刻：应计算时刻：应计算k2，尽管此时，尽管此时yn +1/2已经得到，但已经得到，但un +1则无法得到。则无法得到。(若对若对un +1也进行预报也进行预报加大仿真误差加大仿真误差)。仿。仿真执行延迟真执行延迟h/2输出要迟后半个计算步距。输出要迟后半个计算步距。),(),()(211121211hkyutfkyutfkkkhyynnnnnnnnh / 218感谢你的观看
15、2019年8月3实时龙格库塔法实时龙格库塔法(续）续）RK-2 2的计算流程的计算流程并输出计算1ny1nu采入1k计算nt2htn1nntht21htn21nntht2k计算1k计算下一个19感谢你的观看2019年8月3实时龙格库塔法实时龙格库塔法(续）续） 实时实时2阶龙格库塔法：阶龙格库塔法： tn时刻，应计算时刻，应计算k1，利用当前的，利用当前的un，yn，需要，需要 ；tn+h/2时刻，应计算时刻，应计算k2，此时，此时yn +1/2已经得到，已经得到，un +1/2也也可得到，可得到，k2的计算就不会引入新的误差。计算一次右的计算就不会引入新的误差。计算一次右端函数端函数 需要需
16、要 ，可实时输出，可实时输出yn +1。 )2,(),(12/12/12121khyutfkyutfkhkyynnnnnnnnh / 2fh / 220感谢你的观看2019年8月3实时龙格库塔法实时龙格库塔法(续）续）实时实时RK-2RK-2公式计算流程公式计算流程)2/(htun采入并输出计算1nynt2htn1nntht21htn21nntht1k计算2k计算1k计算下一个21感谢你的观看2019年8月32.3 线性多步法线性多步法 2.3.1 线性多步法基本原理线性多步法基本原理 基本原理基本原理：利用一个多项式去匹配变量若干已知：利用一个多项式去匹配变量若干已知值和它们的导数值。值和它
17、们的导数值。 设：设： 时刻的时刻的 和和 已知已知; 预报：由预报：由 和和 来计算来计算 校正：若校正：若 也已知，由它们来计算也已知，由它们来计算 tttnnn k, 11yyynnn k, 11 , , yyynnn k 11yyynnn k, 11 , , yyynnn k 11yyn kn k, kny yn k22感谢你的观看2019年8月3线性多步法线性多步法(续续) 采用的多项式具有以下形式采用的多项式具有以下形式(m阶阶) 其中：其中： 是待定系数，是待定系数，在在 时刻，时刻， ，可得到：，可得到： （2-1） ytdtthdytidtthidmin kiiiimimmh
18、in kimihiiim( ) ( ) 00111111dihttkntn k0ydn k0 ydn kh 1123感谢你的观看2019年8月3线性多步法线性多步法(续续) 由由 和和 确定确定 ，需要需要m+1个独立方程。该个独立方程。该m+1个方程可由以下个方程可由以下等式导出：等式导出： （2-2） yyynnn k, 11 , , yyynnn k 11),2 , 1 , 0(midiytydtthd jytyidtthid jmn kjin kn kjiiiimimmn kjhin kn kjimihiiim( ) ( ) 0011111124感谢你的观看2019年8月3线性多步法线
19、性多步法(续续)1、预报公式、预报公式 令令m=2k-1，从（，从（2-2）式得到如下方程组：）式得到如下方程组： ddddypppmpn k0121 ddddypppmmpn k0122242 ddddypppmmpn k01223333 ddmdhyppmpn k1212 ddmdhyppmmpn k1212222 ddmdhyppmmpn k1213233 25感谢你的观看2019年8月3 （2-3） 将其写成矩阵形式将其写成矩阵形式： （2-4） 其中上标其中上标p表示预报。表示预报。 其解为其解为： （2-5） V dZpppnknknknnknknknmkkkkmmmmmmyhyh
20、yhyhyyyyddddddddkmkkmmmkkkk32132121121012121232323232103233210223221032101333312222111111dVZpp1p26感谢你的观看2019年8月3 由于由于 为常数阵，其逆存在，为常数阵，其逆存在，Z向量中的各元向量中的各元素为已知值，因而素为已知值，因而d向量的各元素值可计算得到，向量的各元素值可计算得到，从而由从而由 ,得到下一时刻的预报值。得到下一时刻的预报值。 Vpydydn kn kh 011, 缺点：只有缺点：只有 是所需要的，其它元素的计是所需要的，其它元素的计算成为多余算成为多余 ,得不到得不到 与与
21、 和和 显式表达式显式表达式 。dd01,yyn kn k, yyynnn k, 11 , , yyynnn k 1127感谢你的观看2019年8月3线性多步法线性多步法(续续)定义：定义： （m+1） 1的列向量的列向量 (2-6)定义辅助变量定义辅助变量 (2-7)此式可改写为此式可改写为 (2-8)向量向量 的元素可划分为两个组的元素可划分为两个组 (2-9) T)0 , 0 , 1 (0eyn ke deVZ0Tp0Tp1p pTeV0Tp1VepT0 p p pppkpppkpTaaabbb1212,28感谢你的观看2019年8月3例：例：k=3，则，则(2-8)式为：式为：可计算得
22、到：可计算得到： 1110001231111232222312333233123442431235525310000022332244335544123123 aaabbbpppppp pppppppTaaabbb 12312318 910918 3, ,29感谢你的观看2019年8月3 只依赖于只依赖于k，即先前，即先前 和和 的个数，而的个数，而与它们的数值无关。这样，可以预先求解（与它们的数值无关。这样，可以预先求解（2-8）式）式得到得到从而得到从而得到 的显式表达式：的显式表达式： pyn kj yn kj yn kpTe dZ0Tpp yn kya yhb yn kjpn kjjp
23、n kjjkjk 1130感谢你的观看2019年8月3例：试推导用例：试推导用 预报预报 公式公式 条件：已知条件：已知 由由 yyyn kn kn k 112, , yn kVepT0 pVp111012014yyyn kn kn k 112, , pTpppabb112,VpT ppppabb100111124100112 abbppp11213 21 2 ,/ ,/yyhyyn kn kn kn k 11212331感谢你的观看2019年8月32、校正公式、校正公式 预报公式预报公式显式公式，未包括显式公式，未包括 。 校正校正：对该预报值应进行校正，即先预报：对该预报值应进行校正，即先
24、预报得到得到 ，然后再用此值推出，然后再用此值推出 。 由由 和和 以及以及 来预来预报报 ，可令，可令m=2k-1，从（，从（2-2）式得到如下）式得到如下方程组：方程组： yn kyn k yn kyyynnn k, 11,kny , , yyynnn k 11 yn kkncyd032感谢你的观看2019年8月3将其写成矩阵形式：将其写成矩阵形式： （2-10） ddddycccmcn k0121 ddddycccmmcn k0122242 ddddycccmmcn k01223333 ddmdhypcmcn k1212 ddmdhyccmmcn k1212222 ddmdhycpmmc
25、n k1213233 V dZccc33感谢你的观看2019年8月3校正公式校正公式(续续)其中上标其中上标c表示校正，可得表示校正，可得 (2-11) dVZcc1c1000011111122221333310123012232201233332323232121012312mmmmmkkkmkkkkmmmdddddddd yyyyyhyhyhyn kn kn kn knn kn kn k12312334感谢你的观看2019年8月3校正公式校正公式(续续) 定义：定义： 为（为（m+1） 1的列向量，的列向量，上标上标T表示转置。将表示转置。将 左乘（左乘（2-10）式可得：）式可得： （2
26、-12） 定义定义 （2-13） 可改写为可改写为 （2-14） （2-15）e1T ( , , , )0100e1Thyn ke deVZ0Tc1Tc1c cTeV1Tc1VecT1 c ccckccckcTaaabbb0112,35感谢你的观看2019年8月3例：例：k=3 同样，同样， 只依赖于只依赖于k，即先前，即先前 和和 的个数，的个数，而与它们的数值无关。这样而与它们的数值无关。这样 （2-16）从而从而 （2-17） 000010 352552103424421033233210322222101112100001113213214453342232ppppppbbbaaa c
27、yn kj yn kj hyn kcTe dZ1Tcc ()yha yhb yn kjcn kjjcn kjjkjk 11036感谢你的观看2019年8月3例：已知例：已知 ，预估，预估 ，然后用，然后用 校正校正 。预估预估 预报公式为预报公式为校正校正校正公式为校正公式为 yyyn kn kx k 123,kny21,knknknyyy yn k111124139012123 dddyyypppn kn kn k pTpppaaa1233 31,32133knknknknyyyy21210 421111001knknkncccyyyddd 2123210, 2 ,cccTcaaa)2(22
28、11231knknknhknyyyy 37感谢你的观看2019年8月32.3.2 线性多步法误差分析线性多步法误差分析 为了便于分析，对预报公式和校正公式，定义统一的表达为了便于分析，对预报公式和校正公式，定义统一的表达式：式： （2-18） -显式预报显式预报 显式预报显式预报 时称为后向差分公式（时称为后向差分公式（BDF） 同时均不等于同时均不等于0时为隐式校正公式时为隐式校正公式 k称为公式的阶次称为公式的阶次。假设变量各时间的。假设变量各时间的精确值精确值已经得到，已经得到，将其代入（将其代入（21）式，可得：）式，可得： （2-19） kikiikniikniyhy00000 yn
29、 k00yn k021k00,kikiknikniihtyhihty00)()(38感谢你的观看2019年8月3线性多步法误差分析（续）线性多步法误差分析（续） 在在 附近，将每个函数展开成泰勒级数：附近，将每个函数展开成泰勒级数： （2-20） 对所有对所有i(i=0,1,2,k)，将（将（2-20）式代入（）式代入（2-19）式，合并）式，合并同类项，可得同类项，可得 （2-21）tn k )()()()(!2)(! 12knihknihknkntytytyihty )()()()(!2)(! 12knihknihknkntytytyihty )()()()()(2210knmmmknkn
30、kntyChtyChtyhCtyC39感谢你的观看2019年8月3线性多步法误差分析（续）线性多步法误差分析（续）其中其中 （2-22） 如果如果 均为均为0，则称为，则称为p阶公式阶公式 （2-23） kC2100)()2() 1(10!0121! 1111kkkC)2()2() 1(21! 112221!2122kkkkC)2()2() 1(2221!213231!3133kkkkC)2()2() 1(1211)!1(121!1kmmmkmmmmmkkCC C CCp012,)()()()1(1100knpppkikikniknityChihtyhihty40感谢你的观看2019年8月3线
31、性多步法误差分析（续）线性多步法误差分析（续） 下面，如果我们能证明上一节推导出的公式若能下面，如果我们能证明上一节推导出的公式若能满满足求足求 均为均为0的条件，则就得出了这些公式的条件，则就得出了这些公式的截断误差满足（的截断误差满足（2-23）式。）式。 以三阶公式为例，将（以三阶公式为例，将（2-22）式与）式与 相关表达式表示成右端为零，可得：相关表达式表示成右端为零，可得： CC CCp012,32103210,00000000 37271703210362616032103525150321034241403210332313032103222120321011113210000
32、01111321032106665554445533344222332241感谢你的观看2019年8月3 先讨论预报公式，由于先讨论预报公式，由于 和和 ，这意味，这意味着要将上述矩阵的第着要将上述矩阵的第1列移到等式的右边，并去掉列移到等式的右边，并去掉第第5列。为了使矩阵成为方阵，将其最后两行也去列。为了使矩阵成为方阵，将其最后两行也去掉。掉。 100000000000 372717032103626160321035251503210342414032103323130321032221203210111132100000111132103210666555444553334422233
33、2242感谢你的观看2019年8月300000000 37271703210362616032103525150321034241403210332313032103222120321011113210000011113210321066655544455333442223322 其结果与用于推导预报公式的矩阵方程其结果与用于推导预报公式的矩阵方程完全一样。完全一样。 43感谢你的观看2019年8月3 对校正公式，可采用类似的办法，只是对校正公式，可采用类似的办法，只是 ，这样要将，这样要将第第第5列移到等式的右边；若假定列移到等式的右边；若假定 ，则可去掉第，则可去掉第4列。同样列。同样为得
34、到方阵，去掉最后两行，结果就是为得到方阵，去掉最后两行，结果就是3阶校正公式。阶校正公式。 103000000000 37271703210362616032103525150321034241403210332313032103222120321011113210000011113210321066655544455333442223322 这就表明，上一节导出的预报与校正公式的截断误差系这就表明，上一节导出的预报与校正公式的截断误差系数可以用下式来计算数可以用下式来计算 )2()2() 1(21!11211)!1(111kpppkpppppkkC（2-24）44感谢你的观看2019年8月3
35、2.4 稳定性分析稳定性分析 仿真方法选择的基本要求：仿真方法选择的基本要求：仿真计算不改仿真计算不改变原系统的绝对稳定性。变原系统的绝对稳定性。原系统是稳定的。观察欧拉法仿真递推公式原系统是稳定的。观察欧拉法仿真递推公式故有故有 （2-25） yn (n=0,1,2,)为它的一个仿真解。为它的一个仿真解。 iydtdy,/Re 0),(1nnnnythfyynnnyhyy145感谢你的观看2019年8月3稳定性分析（续）稳定性分析（续） 设设 为其准确解，即为其准确解，即 （2-26）用用(2-26)式减去式减去(2-25)式，可得：式，可得：即即 （2-27）特征方程为特征方程为 （2-2
36、8） 显然，为了使扰动序列显然，为了使扰动序列 n不随不随n增加而增长，必须增加而增长，必须要求：要求： 我们称它所对应的域就是该算法的稳定域：我们称它所对应的域就是该算法的稳定域： h 2 1/，即，即h小于等于系统时间常数的两倍。小于等于系统时间常数的两倍。 nny)()()(11nnnnnnyhyynnnh10)1 (1nnhzh()1011 h46感谢你的观看2019年8月3确定数值积分法稳定域的一般方法确定数值积分法稳定域的一般方法 测试方程：测试方程： 数值积分公式数值积分公式 (2-29)其中其中 是一个关于是一个关于 高阶多项式函数，高阶多项式函数，则只有当则只有当 时，算法才稳定。时，算法才稳定。 0 -1 -2ImhReh图图 2.8 2.8 稳定域稳定域iydtdy,/nnyhpy)(1p()hh1)(hp47感谢你的观看2019年8月3
