控制系统CAD控制系统数字仿真

第4章控制系统数字仿真数字仿真就是采用数学模型，在数字计算机上借助数值解法所进行的仿真实验。

所谓数值解法，就是寻求y(t)在[a,b]区间内的一系列离散节点t1<t2<…<t m<t m+1<…上的近似值y1, y2, …，y m, y m+1, …，即求取y m+1 y(t m+1)。

相邻节点的间距h=t m+1-t m称为步长，这里假定h为定值，即t m=t0+mh，m=0, 1, 2,…。

本章主要讲述数字仿真的基本理论和方法。

4.1 数值积分法系统的动态特性通常用一阶微分方程组来描述，也即状态空间表达式。

一般来说，只有极少数的微分方程能用到初等方法得其解析解（或用解析的方法得到精确解），多数只能用近似数值求解。

利用计算机求微分方程主要使用数值积分法，它是系统仿真的一个重要方法。

在这里，我们主要研究一阶微分方程的形式，如：⎪⎩⎪⎨⎧===00)(),(y t t y y t f dtdy（1）， 求y(t)解：⎰⎰⎰+====tt o tt tt dty t f t y t y dty t f dy dty t f dy y t f dtdy),()()(),(),(),(当t=t m+1，t 0=t m 时⎰++=+1),()()(1m mt t m m dty t f t y t y（2）数值积分法时在已知初值的情况下，对f(t, y)进行近似积分，从而对y(t)进行数值求解的方法。

下面介绍几种在数字仿真常用的数字积分法。

1.欧拉法欧拉法又称为折线法，是最简单，也是最早的一种数值计算方法。

对于式（2），如果积分间隔h=t m+1 - t m 取得足够小，使得在t m 与t m+1之间的f(t, y)可近似看做常数f(t m , y m )。

这样式（2）可化为：),()()(1m m m m y t hf t y t y +≈+即),(1m m m m y t hf y y +≈+ （3）（3）式即为欧拉公式。

欧拉公式的几何解释：对于微分方程（1）的解y(t)看作是一条曲线，在任一步长内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长起点的斜率。

基于上述的几何解释，我们从初始点 (t 0, y 0)出发向前推进(t 1, y 1)点，(t 2, y 2)点，…),(),(),(111120001m m m m y t hf y y y t hf y y y t hf y y +=+=+=+M图中阴影部分即为误差。

欧拉法的特点是：计算简单，但精度较低。

例：用欧拉法求解初值问题⎪⎩⎪⎨⎧=-='1)0(2y y t y y （10≤≤t ），h=0.1。

解：因为)2(),(1m mm m m m m m y t y h y y t hf y y -+=+=+则，7848.1)2(1918.1)1.11.021.1(1.01.1)2(1.111.01)2(999910111120001=-+==⨯-+=-+==⨯+=-+=y t y h y y y t y h y y y t y h y y M 该题解为：ty 21+=，将准确解y(t m )与近似解y m 一起放入下表，可得：由此表可以看出欧拉公式的精度很差。

2.后退的欧拉法若用f(t m+1,y m+1)来代替f(t m ,y m )，则（3）式可变为：),(111++++=m m m m y t hf y y （4）则（4）式称为后退的欧拉公式。

后退的欧拉公式是隐式的（因为（4）式右边的y m+1是未知的），此时通常需要用迭代法求解，即：⎪⎩⎪⎨⎧+=+=++++),(),()0(111)0(1m m m m m m m m y t hf y y y t hf y y （5）后退的欧拉公式的几何解释：在任一步长内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长终点的斜率。

例：用后退的欧拉法求解初值问题⎪⎩⎪⎨⎧=-='1)0(2y y t y y （10≤≤t ），h=0.1。

解：⎪⎪⎩⎪⎪⎨⎧-+=+=-+=+=+++++++)2(),()2(),()0(11)0(1)0(111)0(1m m m m m m m m m m m m m m m m y t y h y y t hf y y y t y h y y t hf y y⎪⎪⎩⎪⎪⎨⎧=⨯-⨯+==-⨯+=1918.1)1.11.021.1(1.011.1)101(1.011)0(1y yM⎪⎪⎩⎪⎪⎨⎧=⨯-⨯+==⨯-⨯+=2903.1)2942.12.022942.1(1.01918.12942.1)1918.11.021918.1(1.01918.12)0(2y y后退的欧拉公式和欧拉公式的精度相同，都是一阶精度。

3.梯形法比较欧拉公式和后退的欧拉公式可知，如果对这两种方法进行算术平均，即可大大消除主要误差，从而获得更大的精度，这种方法通常称为梯形法，其计算公式为：)],(),([2111+++++=m m m m m m y t f y t f hy y （6） 同后退的欧拉公式一样，梯形公式也是隐式的。

此时通常采用欧拉公式先预报一个)0(1+m y，再将预报的)0(1+m y代入（6）式进行校正，求出y m+1。

梯形法的迭代公式为：⎪⎩⎪⎨⎧++=⋅+=++++)],(),([2),()0(111)0(1m m m m m m m m m m y t f y t f h y y y t f h y y （7）（7）式又被称为预估-校正公式。

显然，梯形法要比欧拉法和后退的欧拉法精度更高，但计算量比欧拉法大。

例：用梯形法求解初值问题⎪⎩⎪⎨⎧=-='1)0(2y y t y y （10≤≤t ），h=0.1。

解：⎪⎪⎪⎩⎪⎪⎪⎨⎧-+-+=++=-+=+=+++++++)]2()2[(2)],(),([2)2(),()0(11)0(1)0(111)0(1m m m m m m m m m m m m m mm m m m m m m y t y y t y h y y t f y t f h y y y t y h y y t hf y yM⎪⎪⎪⎩⎪⎪⎪⎨⎧=⨯-+⨯-+=-+-+==-⨯+=-+=0959.1)]1.11.021.1()1021[(21.01)]2()2[(21.1)101(1.01)2()0(11)0(1000010000)0(1y t y y t y h y y y t y h y y4.龙格库塔法10),()()(1<<+'=-+θθh t y ht y t y m m m 则：)](,[)()()()(1h t y h t f h t y h t y h t y t y m m m m m m θθθ++⋅+=+'⋅+=+我们把)](,[h t y h t f m m θθ++称作[t m ,t m+1]上的平均斜率，由此可见，只要给平均斜率提供一种算法，便可导出一种计算公式。

当θ=0时，平均斜率为f(t m ,y m )，为欧拉公式；当θ=1时，平均斜率为f(t m+1,y m+1)，为后退的欧拉公式；若采用t m 和t m+1两点斜率值的算术平均作为平均斜率，则平均斜率为)],(),([2111+++m m m m y t f y t f ，此时为梯形公式（精度更高）。

这给我们启示，如果设法在[t m ,t m+1]区间内多预测几个点的斜率值，然后将它们进行加权平均作为平均斜率，则可构造出精度更高的计算公式，这就是龙格库塔的基本思想。

（1）二阶龙格库塔法随意考查区间[t m ,t m+1]区间内一点：ph t t m p m +=+，10≤<p用t m 和t m+p 两个点的斜率值K 1，K 2线性组合得到平均斜率，则⎪⎩⎪⎨⎧⋅+==++=++),(),()(12122111K ph y t f K y t f K K K h y y m p m m m m m λλ（8）因为：mmy y t t m m m p m y f phK t f ph y t f K ph y t f K ==+∂∂+∂∂+=⋅+=)(),(),(112Λ+∂∂∆+∂∂∆+=∆+∆+==0)(),(),(0000y y x x yfy x f x y t f y y x x f 注：所以：mmmmy y t t m m m y y t t m m m m m m y ff t f ph y t f h y y fphK t f ph y t f y t f h y y ====+∂∂+∂∂+++=∂∂+∂∂+++=)(),()(]})(),([),({22211211λλλλλ（8）又因为：+∂∂+∂∂+⋅+=+''⋅+'⋅+===+mmy y t t m m m m mm m y ff t f h y t f h y y hy h y y )(2),(!2221Λ（9）（根据泰勒公式）比较（8）式跟（9）式，可得：⎪⎩⎪⎨⎧==+211221p λλλ（10）满足（10）式的一族公式（8）统称为二阶龙格库塔公式。

当p=1时，此时21,2121==λλ，此时：⎪⎪⎩⎪⎪⎨⎧==++=+++),(),()(21121211m m m m m m y t f K y t f K K K h y y即梯形公式是二阶龙格库塔的一个特例。

另外给出一个常用的二阶龙格库塔公式：（1,0,2121===λλp ）⎪⎪⎩⎪⎪⎨⎧++==+=+)2,2(),(12121K h y h t f K y t f K hK y y m m m m m m（2）三阶龙格库塔法为了进一步提高精度，取[t m ,t m+1]区间内三个点，即t m ，t m+p = t m +ph ，t m+q = t m +qh ，其中，0<p ≤q ≤1，用这三点的斜率值K 1，K 2，K 3的线性组合得到平均斜率，此时：⎪⎪⎩⎪⎪⎨⎧++=⋅+==+++=+++)](,[),(),()(2131213322111sK rK qh y t f K K ph y t f K y t f K K K K h y y m q m m p m m m m m λλλ（11）经过推导可得：⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧==+=+=+=++613121113232232321pqs q p q p s r λλλλλλλλ（12）满足（12）式的一族公式（11）统称为三阶龙格库塔公式。

另外给出一个常用的三阶龙格库塔公式：⎪⎪⎪⎩⎪⎪⎪⎨⎧+-+=++==+++=+)2,()2,2(),()4(62131213211hK hK y h t f K K h y h t f K y t f K K K K h y y m m m m m m m m（3）四阶龙格库塔法同理，取 [t m ,t m+1]区间内四个点的斜率值进行线性组合来得到平均斜率，即为四阶龙格库塔法。