控制系统CAD第4章控制系统数字仿真

1第4章控制系统数字仿真数字仿真就是采用数学模型，在数字计算机上借助数值解法所进行的仿真实验。所谓数值解法，就是寻求y(t)在[a,b]区间内的一系列离散节点t1t2…tmtm+1…上的近似值y1,y2,…，ym,ym+1,…，即求取ym+1y(tm+1)。相邻节点的间距h=tm+1-tm称为步长，这里假定h为定值，即tm=t0+mh，m=0,1,2,…。本章主要讲述数字仿真的基本理论和方法。4.1数值积分法系统的动态特性通常用一阶微分方程组来描述，也即状态空间表达式。一般来说，只有极少数的微分方程能用到初等方法得其解析解（或用解析的方法得到精确解），多数只能用近似数值求解。利用计算机求微分方程主要使用数值积分法，它是系统仿真的一个重要方法。在这里，我们主要研究一阶微分方程的形式，如：200)(),(yttyytfdtdy（1），求y(t)解：ttottttdtytftytydtytfdydtytfdyytfdtdy000),()()(),(),(),(当t=tm+1，t0=tm时31),()()(1mmttmmdtytftyty（2）数值积分法时在已知初值的情况下，对f(t,y)进行近似积分，从而对y(t)进行数值求解的方法。下面介绍几种在数字仿真常用的数字积分法。1.欧拉法欧拉法又称为折线法，是最简单，也是最早的一种数值计算方法。对于式（2），如果积分间隔h=tm+1-tm取得足够小，使得在tm与tm+1之间的f(t,y)可近似看做常数f(tm,ym)。这样式（2）可化为：),()()(1mmmmythftyty即),(1mmmmythfyy（3）4（3）式即为欧拉公式。欧拉公式的几何解释：对于微分方程（1）的解y(t)看作是一条曲线，在任一步长内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长起点的斜率。基于上述的几何解释，我们从初始点(t0,y0)出发向前推进(t1,y1)点，(t2,y2)点，…),(),(),(111120001mmmmythfyyythfyyythfyy图中阴影部分即为误差。欧拉法的特点是：计算简单，但精度较低。…t0t1t2tm+1t0y0y1y2ym+1y…5例：用欧拉法求解初值问题1)0(2yytyy（10t），h=0.1。解：因为)2(),(1mmmmmmmmytyhyythfyy则，67848.1)2(1918.1)1.11.021.1(1.01.1)2(1.111.01)2(9999101111200001ytyhyyytyhyyytyhyy该题解为：ty21，将准确解y(tm)与近似解ym一起放入下表，可得：tmymy(tm)tmymy(tm)0.11.10001.09540.21.19181.18320.31.27741.26490.41.35821.34160.51.43511.41420.61.50901.48320.71.58031.59420.81.64981.61250.91.71781.67331.01.78481.73217由此表可以看出欧拉公式的精度很差。2.后退的欧拉法若用f(tm+1,ym+1)来代替f(tm,ym)，则（3）式可变为：),(111mmmmythfyy（4）则（4）式称为后退的欧拉公式。后退的欧拉公式是隐式的（因为（4）式右边的ym+1是未知的），此时通常需要用迭代法求解，即：),(),()0(111)0(1mmmmmmmmythfyyythfyy（5）后退的欧拉公式的几何解释：…t0t1t2tm+1t0y0y1y2ym+1y…8在任一步长内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长终点的斜率。例：用后退的欧拉法求解初值问题1)0(2yytyy（10t），h=0.1。解：)2(),()2(),()0(11)0(1)0(111)0(1mmmmmmmmmmmmmmmmytyhyythfyyytyhyythfyy91918.1)1.11.021.1(1.011.1)101(1.011)0(1yy2903.1)2942.12.022942.1(1.01918.12942.1)1918.11.021918.1(1.01918.12)0(2yy后退的欧拉公式和欧拉公式的精度相同，都是一阶精度。103.梯形法比较欧拉公式和后退的欧拉公式可知，如果对这两种方法进行算术平均，即可大大消除主要误差，从而获得更大的精度，这种方法通常称为梯形法，其计算公式为：)],(),([2111mmmmmmytfytfhyy（6）同后退的欧拉公式一样，梯形公式也是隐式的。此时通常采用欧拉公式先预报一个)0(1my，再将预报的)0(1my代入（6）式进行校正，求出ym+1。梯形法的迭代公式为：)],(),([2),()0(111)0(1mmmmmmmmmmytfytfhyyytfhyy（7）11（7）式又被称为预估-校正公式。显然，梯形法要比欧拉法和后退的欧拉法精度更高，但计算量比欧拉法大。例：用梯形法求解初值问题1)0(2yytyy（10t），h=0.1。解：)]2()2[(2)],(),([2)2(),()0(11)0(1)0(111)0(1mmmmmmmmmmmmmmmmmmmmmytyytyhyytfytfhyyytyhyythfyy120959.1)]1.11.021.1()1021[(21.01)]2()2[(21.1)101(1.01)2()0(11)0(1000010000)0(1ytyytyhyyytyhyytmymy(tm)tmymy(tm)0.11.09591.09540.21.18411.18320.31.26621.26490.41.34341.34160.51.41641.41420.61.48601.48320.71.55251.59420.81.61531.61250.91.67821.67331.01.73791.7321134.龙格库塔法10),()()(1htyhtytymmm则：)](,[)()()()(1htyhtfhtyhtyhtytymmmmmm我们把)](,[htyhtfmm称作[tm,tm+1]上的平均斜率，由此可见，只要给平均斜率提供一种算法，便可导出一种计算公式。当θ=0时，平均斜率为f(tm,ym)，为欧拉公式；当θ=1时，平均斜率为f(tm+1,ym+1)，为后退的欧拉公式；14若采用tm和tm+1两点斜率值的算术平均作为平均斜率，则平均斜率为)],(),([2111mmmmytfytf，此时为梯形公式（精度更高）。这给我们启示，如果设法在[tm,tm+1]区间内多预测几个点的斜率值，然后将它们进行加权平均作为平均斜率，则可构造出精度更高的计算公式，这就是龙格库塔的基本思想。（1）二阶龙格库塔法随意考查区间[tm,tm+1]区间内一点：phttmpm，10p用tm和tm+p两个点的斜率值K1，K2线性组合得到平均斜率，则15),(),()(12122111KphytfKytfKKKhyympmmmmm（8）因为：mmyyttmmmpmyfphKtfphytfKphytfK)(),(),(11200)(),(),(0000yyxxyfyxfxytfyyxxf注：所以：16mmmmyyttmmmyyttmmmmmmyfftfphytfhyyfphKtfphytfytfhyy)(),()(]})(),([),({22211211（8）又因为：mmyyttmmmmmmmyfftfhytfhyyhyhyy)(2),(!2221（9）（根据泰勒公式）17比较（8）式跟（9）式，可得：211221p（10）满足（10）式的一族公式（8）统称为二阶龙格库塔公式。当p=1时，此时21,2121，此时：18),(),()(21121211mmmmmmytfKytfKKKhyy即梯形公式是二阶龙格库塔的一个特例。另外给出一个常用的二阶龙格库塔公式：（1,0,2121p）19)2,2(),(12121KhyhtfKytfKhKyymmmmmm（2）三阶龙格库塔法为了进一步提高精度，取[tm,tm+1]区间内三个点，即tm，tm+p=tm+ph，tm+q=tm+qh，其中，0p≤q≤1，用这三点的斜率值K1，K2，K3的线性组合得到平均斜率，此时：20)](,[),(),()(2131213322111sKrKqhytfKKphytfKytfKKKKhyymqmmpmmmmm（11）经过推导可得：21613121113232232321pqsqpqpsr（12）满足（12）式的一族公式（11）统称为三阶龙格库塔公式。另外给出一个常用的三阶龙格库塔公式：22)2,()2,2(),()4(62131213211hKhKyhtfKKhyhtfKytfKKKKhyymmmmmmmm（3）四阶龙格库塔法同理，取[tm,tm+1]区间内四个点的斜率值进行线性组合来得到平均斜率，即为四阶龙格库塔法。因其计算很复杂，直接给出最常用的四阶龙格库塔公式：23),()2,2()2,2(),()22(6342312143211hKyhtfKKhyhtfKKhyhtfKytfKKKKKhyymmmmmmmmmm（13）例：用四阶龙格库塔法求解初值问题1)0(2yytyy（10t），h=0.2。解：24334223112143211)(22)2(222)2(222)22(6hKyhthKyKKhyhtKhyKKhyhtKhyKytyKKKKKhyymmmmmmmmmmmmmm251832.1)22(68432.09086.02.01)2.00(29086.022.01)(29086.09182.022.01)1.00(29182.022.012)2(229182.0122.01)1.00(2122.012)2(221102124321013003042002031001020001KKKKhyyhKyhthKyKKhyhtKhyKKhyhtKhyKytyK26tmym(1)ym(2)y(tm)0.10.20.30.40.50.60.70.80.91.01.09591.18411.26621.34341.41641.48601.55251.61531.67821.73791.18321.34171.48331.61251.73211.18321.34161.48321.61251