20-经典的连续系统仿真建模方法学

第二章经典的连续系统仿真建模方法学包括：连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。•问题：数字计算机在数值及时间上的离散性----被仿真系统数值及时间上的连续性连续系统仿真，从本质上：–对原连续系统从时间、数值两个方面对原系统进行离散化并选择合适的数值计算方法来近似积分运算离散模型≈原连续模型？•相似原理:设系统模型为：，其中u(t)为输入变量，y(t)为系统变量；令仿真时间间隔为h，离散化后的输入变量为，系统变量为，其中表示tn=nh。如果，•即，（对所有n=0,1,2,…）则可认为两模型等价。2.1离散化原理及要求),,(tuyfy)(ˆntu)(ˆnty)()(ˆnntutu)()(ˆnntyty0)()(ˆ)(nnnututute0)()(ˆ)(nnnytytyte•对仿真建模方法三个基本要求：（1）稳定性：不改变原系统的稳定性若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。若原连续系统是不稳定的，则离散化后得到的仿真模型也应是不稳定的。（2）准确性：有不同的准确性评价准则，最基本的准则是：--绝对误差准则：–相对误差准则：其中为规定精度的误差量。原连续模型仿真模型u(t)hy(t)-+图2.1相似原理),,(tuyfy),ˆ,ˆ(ˆntuyfy)(ˆnty0)(nyte)()(ˆ)(nnnytytyte)(ˆ)()(ˆ)(nnnnytytytyte•（3）快速性：若第n步计算对应的系统时间间隔为计算机由计算需要的时间为Tn，若Tn=hn称为实时仿真；Tnhn称为超实时仿真；Tnhn称为亚实时仿真，对应于离线仿真•数值积分算法：对已知系统变量y的初始条件求y随时间变化的过程――初值问题计算过程：初始点的,1nnntth)(nty)(1nty),,(tuyfy00)(yty)(ty00)(yty)(00ytf，ttdtytfyty0),()(0f(t,y)f(t0,yo)ttt0t1•欧拉法t1时刻•对任意tn+1时刻•截断误差正比于•梯形法:是隐函数形式。预报－—欧拉法估计初值，校正－—用梯形法校正：预报公式校正公式•反复迭代，直到满足。经典的数值积分法可分为两类：单步法与多步法yytytfty11000()()，)()()(111nnnnnnnytfttytyy，h2)],(),([21)(1111nnnnnnnytfytfhytyy),()(1nnninytfhyy)],(),([21)(11)1(1innnnninytfytfhyyininyy1112.2龙格库塔法•2.2.1龙格-库塔法基本原理•对•若令：•则有•的数值求解：称作“右端函数”计算问题。•经典数值积分法原理：在附近展开台劳级数，只保留项，则有：•(1)•设这个解可以写成如下形式：•其中1),()()(1nnttnndtytftyty)(nntyy1),(QnnttndtytfnnnnyytyQ)(11nQt0h2200010)(21),(htfdtdyyfhytfyytyyakakh101122()),(001ytfkkftbhybkh201021()，•对式右端的函数展成台劳级数，保留h项，可得：•代入，则有：•(2)(2)式与(1)式进行比较，可得：•四个未知数但只有三个方程，因此有无穷多个解。若限定，则•计算公式：•其中•若写成一般递推形式，即为：k2hyfkbtfbytfkt0)(),(121002])(),([),(012100200101hyfkbtfbytfhaythfayytaaabab12212211212+===,/,/，，，，2121bbaaaa12aabb1212121，)(22101kkhyy)(),(1002001hkyhtfkytfk，，)(2)(2111kkhyytynnn•其中•截断误差正比于h3，称为二阶龙格-库塔法（简称RK-2）。•截断误差正比于h5的四阶龙格--库塔法（简称RK-4）公式：•其中：••龙格--库塔法的特点：•1.形式多样性–例：对二阶，非唯一解，可以得到许多种龙格--库塔公式：),(),(121hkyhtfkytfknnnn，)22(61)(432111kkkkyytynnn),(1nnytfk)22(12khyhtfknn，)22(23khyhtfknn，)(34hkyhtfknn，2121bbaa，，，hkyynn21•其中各种龙格---库塔法可以写成如下一般形式：•其中••式中各系数满足以下关系•s称为级数，表示每步计算右端函数的最少次数。可以证明，1阶公式至少要计算一次，2阶公式；….；4阶公式；依此类推。有时为了某种特殊需要，可以选择的计算公式。),(1nnytfk)22(12khyhtfknn，siiinnkChyy11)(11ijjijninikbhyhatfk，si，，，21aabisCiijjiiis11110231，，，fsmin2smin4minss•2.单步法，可以自启动，存储量小–在计算时只用到，而不直接用等项。•3.可变步长–步长h在整个计算中并不要求固定，可以根据精度要求改变–但是在一步中，为计算若干个系数，则必须用同一个步长h。•4.速度与精度–四阶方法的h可以比二阶方法h的大10倍，每步计算量仅比二阶方法大一倍–高于四阶的方法由于每步计算量将增加较多，而精度提高不快。1nyny21nnyy，ik2.2.2实时龙格－库塔法–实时仿真：要求仿真模型的执行速度与实际系统运行的速度保持一致–一般的数值积分法难以满足实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速–速度较慢，而且这些方法的机理不符合实时仿真的特点。–考虑系统–分析RK-2公式：–一个计算步内分两子步：–tn时刻：利用当前的un，yn计算k1----计算一次右端函数f需h/2–tn+h/2时刻：应计算k2，尽管此时yn+1/2已经得到，但un+1则无法得到。(若对un+1也进行预报――加大仿真误差)。仿真执行延迟h/2――输出要迟后半个计算步距))((tuyfdtdy，),,(),,()(211121211hkyutfkyutfkkkhyynnnnnnnn图2.4RK-2的计算流程1nu采入并输出计算1ny1k计算2k计算1k计算下一个nt2htn1nntht21htn21nntht•实时2阶龙格－库塔法：•tn时刻：应计算k1，利用当前的un，yn，需要；tn+h/2时刻，应计算k2，此时yn+1/2已经得到，un+1/2也可得到，k2的计算就不会引入新的误差。计算一次右端函数f需要h/2，可实时输出yn+1。)2,,(),,(12/12/12121khyutfkyutfkhkyynnnnnnnn2.3线性多步法2.3.1线性多步法基本原理基本原理：利用一个多项式去匹配变量若干已知值和它们的导数值。设：时刻的和已知tttnnnk,,,11yyynnnk,,,11,,,yyynnnk11yyynnnk,,,11,,,yyynnnk11yynknk,knyynkytdtthdytidtthidminkiiiimimmhinkimihiiim()()00111111预报：由和来计算校正：若也已知，由它们来计算采用t为自变量的多项式，具有以下形式(m阶)：ditthnktnkydnk0ydnkh11yyynnnk,,,11,,,yyynnnk11diytydtthdjytyidtthidjmnkjinknkjiiiimimmnkjhinknkjimihiiim()()00111111其中：是待定系数，＝在时刻，＝0，可得到：由和确定需要m+1个独立方程。该m+1个方程可由以下等式导出：(1)(i=0,1,2,…m)，（j=1,2,...,k）（2）ddddypppmpnk0121ddddypppmmpnk0122242ddddypppmmpnk01223333ddmdhyppmpnk1212ddmdhyppmmpnk1212222ddmdhyppmmpnk12132331、预报公式令m=2k-1，从（2）式得到如下方程组：nknknknnknknknmkkkkmmmmmmyhyhyhyhyyyyddddddddkmkkmmmkkkk32132121121012121232323232103233210223221032101333312222111111VdZppp将其写成矩阵形式：其中上标p表示预报。(3)（4）dVZpp1pVpydydnknkh011,dd01,yynknk,yyynnnk,,,11,,,yyynnnk11T)0,,0,1(0eynkedeVZ0Tp0Tp1ppTeV0Tp1VepT0pppppkpppkpTaaabbb1212,,,,,,其解为：由于因而d向量的各元素值可计算得到，从而由,得到下一时刻的预报值。是所需要的，其它元素的计算成为多余与和显式表达式（m+1）1的列向量定义辅助变量此式可改写为向量的元素可划分为两个组(9)缺点：（1）只有（2）得不到定义：（6）(7)(8)(5)为常数阵，其逆存在，Z向量中的各元素为已知值，1110001231111232222312333233123442431235525310000022332244335544123123aaabbbpppppp123123,,,,,18,9,10,9,18,3TTpppppppaaabbbpynkjynkjynkpTedZ0Tppynkyayhbynkjpnkjjpnkjjkjk11例：k=3，则(8)式为：可计算得到：只依赖于k，即先前和的个数，而与它们的数值无关。从而得到的显式表达式：(11)这样，可以预先求解（8）式得到。(10)yyynknknk112,,ynkyyynknknk112,,pTpppabb112,,VepT0pVp111012014VpTppppabb100111124100112abbppp11213212,/,/yyhyynknknknk112123ynk~ynk