控制系统仿真及MATLAB语言-连续系统的离散化方法

1、第四章第四章 连续系统的离散化方法连续系统的离散化方法4.1 常微分方程的数值解法00( , )( )dxf t xdtx tx一一. 数值求解的基本概念数值求解的基本概念设微分方程为则求解方程中函数x(t)问题的常微分方程初值问题所谓数值求解就是要在时间区间a, b中取若干离散点 求出微分方程在这些时刻的近似值012Nx x xx01Natttb2(2)( )00000()( )( )( )( )2!kkhhx thx th x txtxtknnntth1取前两项近似：1( ,)kkkkxhf txx这种方法的几何意义就是把f(t,x)在区间tk,tk+1内的曲边面积用矩形面积近似代替。计算

2、简单，计算量小，而且可以自启动。当h很小时，造成的误差是允许的。该算法具有一阶精度。取k=0,1,2,N,从t0开始，逐点递推求解t1时的y1， t2时的y2,直至tn时的yn，称之为欧拉递推公式。矩形面积矩形面积1. 欧拉法欧拉法 欧拉法的特点：导出简单，几何意导出简单，几何意义明显，便于理解，能说明构造数义明显，便于理解，能说明构造数值解法一般计算公式的基本思想值解法一般计算公式的基本思想。通常用它来说明有关的基本概念。 例例 设系统方程为 用Euler法求其数值解（取步长 , )1 . 0h 2,01xxx 10 t递推公式为1,10.1nnnnnnxxhftxxx则4628. 01 .

3、 01, 0 . 17519. 0819. 01 . 01819. 01 . 01, 3 . 0819. 091. 09 . 01 . 01, 2 . 09 . 01 . 01, 1 . 01, 099101022331122001100 yyytyyytyyytyyytyt已知方程的解析解为 精确解和解析解作比较： 误差在 数量级，精度较差。 210ty11t00.10.20.30.40.51.0精确解10.909 0.833 0.769 0.666 0.625 0.5数值解10.90.819 0.752 0.659 0.6470.4632. 龙格库塔法龙格库塔法 基本思想：基本思想：取Ta

4、ylor级数展开式前三项近似求解，并利用线性组合代替导数的求解。 既可避免计算高阶导数，又可提高数值积分的精度，这就是Runge-Kutta法的基本思想。2. 龙格库塔法龙格库塔法kktxffr为精度阶次，为精度阶次，ai为待定系数，由精度确定；为待定系数，由精度确定；ki用下用下式表示式表示线性组合线性组合1121(,) ,2,3,iikkjjkf tb h xhbkir21()2!kkkkktxkhxxh ffff11rkkiiixxha k2(2)0000()( )( )( )2hx thx th x txt等各阶导数不易计算，用下式中等各阶导数不易计算，用下式中kiki的线性组合代替的

5、线性组合代替1)当r=1时：11 1,kkxxha k1( ,),kkkf tx与Taloy展开式相比较，可得a1=1,则上式成为11( ,),kkkkkxxhkxhf tx欧拉递推公式2)当r=2时：12121kkkkkft ,xkftb h,xhb k2121kkktxkfbhfb hk f12 1kkf tbh,xhbkkkt , x将 在点 展成Taylor级数与台劳公式的二阶展开近似公式相比，可得以下关系：与台劳公式的二阶展开近似公式相比，可得以下关系：12222 111 21 2aab aa b三个方程，四个未知数，解不唯一三个方程，四个未知数，解不唯一各个系数的几种取法各个系数的

6、几种取法见书上。见书上。3) r=4时，四阶龙格库塔公式时，四阶龙格库塔公式-最常用：最常用：仿真中遇到的大多数工程实际问题，四阶龙格库塔法以能满足精度要求，其截断误差o(h5) 与h5同数量级。该法可以自启动。1123412132432262222kkkkkkkkkkhxx(KKKK )Kf t ,xhhKft,xKhhKft,xKKf th,xhK 4)、状态空间四阶龙格、状态空间四阶龙格-库塔递推式库塔递推式 若单输入单输出系统的状态空间表达式为：若单输入单输出系统的状态空间表达式为：XAXBUYCXDU 在仿真中，对于在仿真中，对于n阶系统，状态方程可以写成一阶微分方程阶系统，状态方程

7、可以写成一阶微分方程1 122123( ,)(1,2, )iiiinniinxa xa xa xbuf t x x xxin 根据四阶龙格根据四阶龙格-库塔公式库塔公式,有有1123411 122211113112241133(22)(1,2, )6( )()()()222()()()222()(kkiiiiiikkkiiiinnikkkiiiinniikkkiiiinniikkkiiiinnihxxkkkkinka xa xa xbu thhhkaxkaxkbu thhhkaxkaxkbu tkaxhkaxhk)()ikbu thT=tT=tk k时刻的时刻的x xi i值值T=tT=tk

8、k+h+h时刻的时刻的x xi i值值 另另111121211211212222313233414244 , , ,kkkTkkk TnnTTnnTTnnxxxxxxkkkkkkkkkkkkk+1k1xxKKKK 状态方程的四阶龙格状态方程的四阶龙格-库塔公式如下：库塔公式如下：23423243(22)6( )()()22()()22()()kkkkhu thhu thhu thu thk+1k11kk1kkk+1k+1xxKKKKKAxBKA xKBKA xKBKA xKByCxRK法的特点：法的特点： 1 需要存储的数据少，占用的存储空间少；2 只需知道初值，即可启动递推公式进行计 算，可

9、自启动；3 容易实现变步长运算。4 每积分一步需要计算多次右函数，计算量 大。基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：t, x=ode23(xfun, t0, tf ,x0)t, x=ode45(xfun, t0, tf ,x0)常微分方常微分方程函数名程函数名起始起始时间时间终止终止时间时间初始状初始状态向量态向量输入输入输出输出4/54/5阶龙格阶龙格- -库塔算法库塔算法2/32/3阶龙格阶龙格- -库塔算法库塔算法3.3.常微分方程常微分方程MatlabMatlab求解求解解解: 令令 y1=x，y2=x1、建立、建立M-文件文件vdp.m如下：如下

10、： function dy=vdp(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=2*(1-y(1)2)*y(2)-y(1); 2、取、取t0=0，tf=20，输入命令：，输入命令： T,Y=ode45(vdp,0 10,1;1); plot(T,Y(:,1),-, T,YT,Y(:,2)3、结果、结果1)0(, 1)0(, 0)0(51. 0321213312321yyyyyyyyyyyy解解 1、建立、建立m-文件文件rigid.m如下：如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(

11、2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);2、取取t0=0，tf=12，输入命令：，输入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)3、结果如图、结果如图图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.4.2 数值算法的稳定性及求解原则1 1. .数值算法的稳定性数值算法的稳定性 特征根在特征根在s s平面的左半平面，系统稳定。平面的左半平面，系统稳定。(1)欧拉法： 稳定：(2)梯形法： 恒稳11zh( )(1)( )zX zhX z1z (

12、)(1)(1)( )22hhzX zX z2 2. .数值算法的选择原则数值算法的选择原则 Matlab提供了微分方程数值求解的一般方法,作为仿真算法的使用者,可不必考虑算法具体实现,而应关心各种方法在使用中会出现的问题,以及如何在仿真中恰当的选用这些方法. 一般,选用数值算法从以下几个方面考虑:(1)精度受算法和h影响截断误差+舍入误差=累计误差(2)计算速度受算法和h影响算法简单，速度就快些。(3)稳定性 受h影响，一般 h（2-3）系统最小时间系统最小时间4.3 数值算法中的数值算法中的“病态病态”问题问题 1 “病态”常微分方程例：(0)(1,0, 1)TXAXX其中123( ,)TT

13、Xx x x-211920A= 192120404040采用四阶龙格库塔法h=0.01时，计算时间长时，计算时间长h=0.04时，误差很大时，误差很大当当h0.05后，曲线发散振荡，数值不稳定，完全失去意义后，曲线发散振荡，数值不稳定，完全失去意义一般线性常微分方程组：一般线性常微分方程组：00( )( )( ),( )X tAX tBU tX tX的系数矩阵的系数矩阵A的特征值具有如下特征：的特征值具有如下特征：11Re()0max Re()min Re()iiii ni n 则称为则称为“病态病态”方程。方程。定义：定义：2 控制系统仿真中的“病态”问题1 病态系统中绝对值最大的特征值对应

14、于系统动态性能解中瞬态分量衰减最快的部分，它反映了系统的动态响应和系统的反应灵敏度。一般与系统中具有最小时间常数Tmin的环节有关，要求计算步长h取得很小。2 病态系统中绝对值最小的特征值对应于系统动态性能解中瞬态分量衰减最慢的部分，它决定了整个系统的动态过渡过程时间的长短。一般与系统中具有最小时间常数Tmax的环节有关，要求计算步长h取得很大。3 对于病态问题的仿真需要寻求更加合理的算法，以解决病态系统带来的选取计算步长与计算精度，计算时间之间的矛盾。3 “病态”系统的仿真方法 采用稳定性好，计算精度高的数值算法，并且允许计算步长能根据系统性能动态变化的情况在一定范围内作相应的变化，采用隐式

15、吉尔法101111kkkrk rkxa xa xa xhf该法已经证明对病态方程求解过程是数值稳定的。隐式吉尔法从理论上十分适应于病态系统 ，但需要解决好以下问题(1) 自启动 r阶多步算式无法自启动，需要用单步法求出前r步值(2) 预估迭代 迭代方法要求收敛性良好，否则在大步长时会造成数 值发散。(3) 变步长 初始阶段采用小步长，随后可逐步放大步长。 对不同精度要求的系统仿真，要考虑变阶次问题，即为减小每一步计算的截断误差，以提高精度，应选用较高的阶次，而当精度较低时，为减少工作量，则应选取较低的阶次。仿真时应根据估计误差 与给定的误差精度相比较改变步长或阶次来重新计算。 4.4 连续系统

16、状态方程的离散化上章所述的连续系统数学模型的离散化，上章所述的连续系统数学模型的离散化，是通过数值积分法实现的，尽管面向结是通过数值积分法实现的，尽管面向结构图的仿真方法是按环节给定参数，但构图的仿真方法是按环节给定参数，但是在计算时还是是在计算时还是按整个系统进行离散化按整个系统进行离散化，这就不便于引进非线性环节以进行非线这就不便于引进非线性环节以进行非线性系统的仿真。在本节，将介绍连续系性系统的仿真。在本节，将介绍连续系统离散模型的建立和仿真。统离散模型的建立和仿真。开 始输入开环传函分母、分子系数jica ,输入反馈函数 V、状态初始值 X0求开环状态方程系数阵 A,B,C求闭环状态方

17、程系数阵BvCAAb输入初始时间0t、终止时间ft、计算步长 h、输入幅值 r求龙格库塔法各次斜率4 , 3 , 2 , 1, jKj求11,kkyX输出数据、曲线1kkXXft到否？NY结 束 h改变时，叠代过程重复求解，费时繁琐不能对非线性环节单独考虑。连续系统离散化连续系统离散化思想：用差分方程描述连续系统的状态方程模型思想：用差分方程描述连续系统的状态方程模型（因为差分方程的主要特点就是方程中各变量由各相邻时刻的变（因为差分方程的主要特点就是方程中各变量由各相邻时刻的变 化量制约，这相当于递推方程）化量制约，这相当于递推方程）1 1、连续系统的离散化、连续系统的离散化DuCXyBuAX

18、X 设连续系统状态方程为设连续系统状态方程为00)(XtX其中其中为状态初始值为状态初始值.则由现代控制理论基础知，状态变量则由现代控制理论基础知，状态变量X(t)的解为的解为()00()( )kTAkTA kTkTeedXXButtAtAdee0)(0)()(uBXtX而而 t = ( k + 1 )T 时，可表示为时，可表示为()0(1) )() ()TTTkTekTedkTAAXXBu当系统输入当系统输入u(tu(t) )给定时，可求出系统离散化状态方程给定时，可求出系统离散化状态方程的解。一般，的解。一般， u(tu(t) )未知，通常采用两种方法近似处未知，通常采用两种方法近似处理：

19、理：(1)(1)令令u(kT+t)u(kT+t)u(kTu(kT) (0) (0t tT)T)相当在系统输入端加一个采样开关和零阶保持器相当在系统输入端加一个采样开关和零阶保持器 X X(k k+1)+1)T T) = ) = GXGX( (kTkT) + ) + HuHu( (kTkT) )G G = = e e A TA T，为，为t t = = T T 时的状态转移矩阵时的状态转移矩阵()0TTHedAB离散后的状态空间表达式为：离散后的状态空间表达式为：(1) ()()()()()x kTGx kTHu kTy kTCx kTDu kT、MatlabMatlab表示表示DuCXyBuA

20、XX已知连续系统状态方程为已知连续系统状态方程为在采样周期在采样周期T下离散后的状态空间表达可表示为：下离散后的状态空间表达可表示为：(1) ()()()()()x kTGx kTHu kTy kTCx kTDu kT()0G,TTTeHedAAB其中，在在MatlabMatlab中，若已知连续系统状态方程各阵模型参数中，若已知连续系统状态方程各阵模型参数 ( (A A、B B、C C、D D) ) 以及采样周期以及采样周期T T，则语句：，则语句： G G，H H = c2d ( = c2d (A A, ,B B, ,T)返回的矩阵返回的矩阵G G 、H H 就是所要求的就是所要求的( (

21、T T ) ) 、mm ( ( T T ) ) 。此外，此外， MatlabMatlab还提供了功能更强的求取连续系统离还提供了功能更强的求取连续系统离散化矩阵函数散化矩阵函数c2dm(),c2dm(),他容许调用时选用离散化变换他容许调用时选用离散化变换方式，并且得到的是标准的离散化状态方程。方式，并且得到的是标准的离散化状态方程。 G G, ,H H, ,C C, ,D D=c2dm (=c2dm (A A, ,B B, ,C C, ,D D, ,T T,选项选项) ) 表表 离散化变换方式选项离散化变换方式选项选选 项项说说 明明Zoh假设输入端加一个采样开关和零阶保持器假设输入端加一个采样开关和零阶保持器Foh假设输入端加一个采样开关和一阶保持器。假设输入端加一个采样开关和一阶保