系统仿真技术_第2章 经典的连续系统仿真建模方法学

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