连续系统的离散化方法课件

1、第四章第四章 连续系统的离散化方法连续系统的离散化方法4.1 常微分方程的数值解法常微分方程的数值解法4.1.1 数值求解的基本概念数值求解的基本概念已知一个一阶微分方程已知一个一阶微分方程 00( , )( )xf t xx tx 0 x10211,nntth tthtth12,nx xx10tth0()x th应用数值求解的思路是，从初值应用数值求解的思路是，从初值开始，在一系列时刻开始，在一系列时刻，求未知解，求未知解的近似解，的近似解，h是计算步长。若是计算步长。若x对对t的各阶导数都存在。则的各阶导数都存在。则x(t)在在用台劳级数表示为：用台劳级数表示为：2(2)( )00000(

2、)( )( )( )( )2!kkhhx thx thx txtxtk4.1.2 000,( )ttx tx10tth11( )x tx初始条件为初始条件为， 欧拉法的计算公式为欧拉法的计算公式为 1000( ,)xxhf tx其一般公式为 1( ,)kkkkxxhf tx0tf0f1fc0t1t23(2)(1)000000( ,)( ,)( ,)2!3!kknhhhRf txftxftxk称为截断误差称为截断误差 例例41 用欧拉法求下述微分方程的数值解。用欧拉法求下述微分方程的数值解。2(0)1xxx 解：因欧拉法的递推公式为解：因欧拉法的递推公式为 1( ,)kkkkxxhf tx所以所

3、以 2( ,)kkkf txx 则递推公式为：则递推公式为： 21kkkxxhx若取步长若取步长 0.1h 由由t=0开始计算可得开始计算可得 212223101 0.1 10.90.90.1 0.90.8190.8190.1 0.8190.75190.4627810 xxxx 上式的精确解是上式的精确解是 11xt数值计算与精确解的比较见表数值计算与精确解的比较见表 t00.10.20.31.0精确解x(t)10.90909090.83333330.76923070.510.90.8190.57190.46278104.1.3 龙格库塔法龙格库塔法2(2)0000()( )( )( )2!h

4、x thx thx txt100201021101122( ,)(,)()Kf txKf tbh xb K hxxh a Ka K2100000( ,)(),2hffxxhf txfttxxtx2K100( ,)Kf tx将将在点在点处作台劳级数展开，并取线性部分可得处作台劳级数展开，并取线性部分可得 20012100( ,)ffKf txbhb Khttxxtx将将 1010020012100( ,) ( ,)ffxxa hf txa h f txbhb hKttxxtx1K2K代入式代入式 比较比较 各项系数得各项系数得 122 122111,22aaa ba b待定系数个数超过方程个数，

5、必须先设定一个系数，然后即可求得其待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其参数。一般有以下几种取法参数。一般有以下几种取法：1、121210,1,2aabb则则1021002001( ,)(,)22xxhKKf txhhKf txK一般形式一般形式 12121( ,)(,)22kkkkkkxxhKKf txhhKf txK2、1212132,443aabb则则10121002001(3)4( ,)22(,)33hxxKKKf txKf th xhK 一般形式一般形式 112121(3)4( ,)22(,)33kkkkkkhxxKKKf txKf th xhK3、 12121

6、,12aabb则则 10121002001()2( ,)(,)hxxKKKf txKf th xhK 一般形式一般形式 112121()2( ,)(,)kkkkkkhxxKKKf txKf th xhK以上几种递推公式均称为二阶龙格库塔公式。是较典型的几以上几种递推公式均称为二阶龙格库塔公式。是较典型的几个常用算法。其中的方法个常用算法。其中的方法3又称为预估校正法，或梯形法。又称为预估校正法，或梯形法。 其意义如下：用欧拉法以斜率先求取一点，其意义如下：用欧拉法以斜率先求取一点， 再由此点求得另一斜率再由此点求得另一斜率101xxhK211001( ,)(,)Kf t xf th xhK0

7、x1K2K122KKK10012()2hxxhKxKK然后，从然后，从点开始，既不按该点斜率点开始，既不按该点斜率变化，也不按预估点斜率变化，也不按预估点斜率变化，而是取两者平均值变化，而是取两者平均值求得校正点，即：求得校正点，即：。0tf0f1f0t1t1001()2hxxff四阶龙格库塔法的计算公式为：四阶龙格库塔法的计算公式为：11234(22)6kkhxxKKKK1213243( ,)(,)22(,)22(,)kkkkkkkkKf txhhKf txKhhKf txKKf th xhK对于用状态方程表示的高阶线性系统对于用状态方程表示的高阶线性系统 XAXBUYCX0(0)XX其中状

8、态变量为其中状态变量为 12,nXx xx用四阶龙格库塔法时，有用四阶龙格库塔法时，有 ( ,)f t XAXBU计算公式为：计算公式为： 11234(22)6kkhXXKKKK其中，其中，1211322433( ,)( )(,)()()2222(,)()()2222(,)()()kkkkkkkkkkkkkkkkKf tXAXBU thhhhKf tXKA XKBU thhhhKf tXKA XKBU tKf th XhKA XhKBU th相应的输出为相应的输出为 11kkYCX按上式，取按上式，取0,1,2,kn不断递推，不断递推， 即可求得所需时刻即可求得所需时刻 11, ,nt tt(

9、 )kX t( )kYt的状态变量的状态变量和输出值和输出值hiK德国学者德国学者Felhberg对传统的龙格库塔法提出了改进，在每一个计算步长内对传统的龙格库塔法提出了改进，在每一个计算步长内对对f( )函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的步长为步长为，定义下面，定义下面6个个变量变量11(,),1,2,6;1,2,1iikijjkijKhf xK th iji下一步的状态变量可由下式求出：下一步的状态变量可由下式求出：611kkiiixxK 四阶四阶/五阶龙格库塔法系数表五阶龙格库塔法系数表iiji*i01

10、6/13525/2161/41/4003/83/329/326656/128251408/256512/131932/2197-7200/21977296/219728561/564302197/41041439/216-83680/513-845/4104-9/50-1/51/2-8/272-3544/25651859/4104-11/402/550这一方法又称为四阶这一方法又称为四阶/五阶龙格库塔法。五阶龙格库塔法。 4.1.4 微分方程数值解的微分方程数值解的MATLAB实现实现1、,23(, 0,)t xodeFun Tspan xoptions,45(, 0,)t xodeFun T

11、span xoptions该指令适用于一阶常微分方程组该指令适用于一阶常微分方程组 ( ,),1,2,iiixf t xin如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方程方可使用。程方可使用。2、输入参数、输入参数Fun为定义微分方程组为定义微分方程组 ( ,) ,1, 2 ,iiixftxinM函数文件名，可以在文件名加写函数文件名，可以在文件名加写，或用英文格式单引号界定文件名。，或用英文格式单引号界定文件名。3、在编辑调试窗口中编写一阶常微分方程组在编辑调试窗口中编写一阶常微分方程组 ( ,) ,1, 2

12、,iiixft xin的的M函数文件时，每个微分方程的格式必须与函数文件时，每个微分方程的格式必须与 ( ,)iiixf t x,it xix1,2,in一致，即等号一致，即等号严格以严格以“先自变量先自变量t，后函数，后函数”的固定顺序输入的固定顺序输入，表示微分方程的序数。表示微分方程的序数。左边为带求函数的一阶导数，右边函数的变量左边为带求函数的一阶导数，右边函数的变量4、输入参数、输入参数“Tspan”规定了常微分方程的自变量取值范围，规定了常微分方程的自变量取值范围，它以矩阵它以矩阵t0,tf的形式输入，表示自变量的形式输入，表示自变量0,tttf5、输入参数、输入参数x0表示初始条

13、件向量，表示初始条件向量， 0( 0)( 0)( 0)xx tx tx t微分方程组中的方程个数必须等于初始条件数，这是求微分方程微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解所必须的条件。特解所必须的条件。6、输入参数、输入参数options表示选项参数（包括表示选项参数（包括tol，trace），可缺省，），可缺省，即取默认值，即取默认值，tol是控制结果精度的选项对是控制结果精度的选项对ode23( )函数取函数取 ，对对ode45( )函数取函数取 。trace为输出形式控制变量，如果为输出形式控制变量，如果trace不为不为0，则，则 会将仿真中间结果逐步地由频幕显示出

14、来，否则将不会将仿真中间结果逐步地由频幕显示出来，否则将不 显示中间结果显示中间结果610310itit( )jixg t7、输出参数、输出参数t,x为微分方程组解函数的列表（为微分方程组解函数的列表（t和和x都是列矩阵），它都是列矩阵），它包含向量包含向量t各节点各节点 和与和与对应向量对应向量x的第的第j个分量值个分量值（即第（即第j个方程解），个方程解），i表示节点序列数。表示节点序列数。( )xg t()( )jjxgt8、输出参数、输出参数t,x缺省时，输出解函数的曲线，即函数缺省时，输出解函数的曲线，即函数及其各及其各的曲线。的曲线。阶导数阶导数求解微分方程的指令还有求解微分方程的

15、指令还有ode113（多步解法器），（多步解法器），ode15s（基于数字（基于数字微分公式的解法器），微分公式的解法器），ode23s（单步解法器），（单步解法器），ode23T（梯形规则的（梯形规则的一种自由插值实现），一种自由插值实现），ode23TB（二阶隐式龙格库塔公式）等。（二阶隐式龙格库塔公式）等。例例 41 求解常微分方程求解常微分方程 220txx，初始条件为初始条件为： (0)1;(0)1;0,10 xxt 解：方法解：方法1 把二阶微分方程化成两个一阶微分方程组：把二阶微分方程化成两个一阶微分方程组： 令令 12,xxxx则：则： 122212xxtxx 首先编制首先编制

16、M文件，并且函数名和文件，并且函数名和M文件名相同。文件名相同。 function xdot=wffc_1(t,x) %定义输入、输出变量和函数文件名定义输入、输出变量和函数文件名xdot=zeros(2,1); %明确明确xdot的维数的维数xdot (1)=x(1); % 第一个微分方程表示形式第一个微分方程表示形式xdot(2)=-x(1)+2+t2/pi; % 第二个微分方程表示形式第二个微分方程表示形式方法方法2 写出系统的状态方程写出系统的状态方程2010(2)101txx 首先编制首先编制M文件，并且函数名和文件，并且函数名和M文件名相同。文件名相同。function xdot=

17、wffc_1(t,x) %定义定义t,x,xdot和文件名和文件名xdot=0 1;-1 0*x+0;1*(2+t2/pi); %状态方程的表示形式状态方程的表示形式在命令窗口键入在命令窗口键入t,x=ode45(wffc_1,0,10,-1;1)，可得微分方程的数值解，可得微分方程的数值解，其前其前10组数据如下：组数据如下：t = 0 0.0167 0.0335 0.0502 0.0670 0.1507 0.2344 0.3182 0.4019 0.5964 -x = -1.0000 1.0000 -0.9828 1.0501 -0.9648 1.0999 -0.9460 1.1494 -

18、0.9263 1.1986 -0.8158 1.4395 -0.6856 1.6709 -0.5363 1.8917 -0.3691 2.1007 0.0830 2.5349 -22(1)0 xxxx(0)1;(0)1xx0,20t例例 42 求求Var der Pol微分方程微分方程，在初始条件，在初始条件下的数值下的数值。解解解：将二阶微分方程变换成一阶微分方程组解：将二阶微分方程变换成一阶微分方程组 12,xx xx 则则 12221212*(1)xxxxxx编制编制M文件，并且函数名和文件，并且函数名和M文件名相同。文件名相同。function xdot=vdpl(t,x)xdot=z

19、eros(2,1); 赋初值，并规定向量的维数。赋初值，并规定向量的维数。xdot(1)=x(2); 对第一个微分方程进行描述对第一个微分方程进行描述xdot(2)=2*(1-x(1)2)*x(2)-x(1); 对第二个微分方程进行描对第二个微分方程进行描述述在命令窗口键入在命令窗口键入t,x=ode45(vdpl,0,20,1;1)，可得微分方程的，可得微分方程的数值解，其前数值解，其前10组数据如下：组数据如下：t = 0 0.0502 0.1005 0.1507 0.2010 0.3136 0.4262 0.5389 0.65150.7484-x = 1.0000 1.0000 1.04

20、89 0.9437 1.0946 0.8762 1.1368 0.7995 1.1748 0.7158 1.2441 0.5157 1.2908 0.3156 1.3159 0.1311 1.3215 -0.0278 1.3131 -0.1431 -若在命令窗口输入若在命令窗口输入ode45(wffc_1,0,10,-1;1)，则可得到系统状，则可得到系统状态曲线，而不输出数据。如图态曲线，而不输出数据。如图44所示。所示。 该系统也可直接用该系统也可直接用SIMULINK进行仿真。首先建立进行仿真。首先建立SIMULINK仿真模型仿真模型如图如图45所示。所示。设置系统参数如下，将积分器初值

21、设置为系统要求的初值，如图设置系统参数如下，将积分器初值设置为系统要求的初值，如图46所示。所示。 XY示波器参数设置如图示波器参数设置如图47所示。所示。仿真时间参数设置如图仿真时间参数设置如图48所示。所示。 系统仿真曲线如图系统仿真曲线如图49和和410所示。所示。 4.2数值解法的稳定性及选择原则数值解法的稳定性及选择原则4.2.1 数值解法的稳定性数值解法的稳定性欧拉法欧拉法对对 0(0)xxxx按欧拉公式计算得按欧拉公式计算得 1(1),0,1,1kkkkxxh xhxkn这是一个一阶差分方程，这是一个一阶差分方程，Z变换后得变换后得 ( )(1)( )zX zhX z，其特征值为

22、，其特征值为 1zh 11zh11zh2h则该算法的稳定条件为则该算法的稳定条件为2h。算法的稳定域如图。算法的稳定域如图411所示。所示。 0-1-2Re()hIm()h不稳定域2.梯形法梯形法111()()22kkkkkkkhhxxffxxx( )(1)(1)( )22hhzX zX z1212hzh1212hzh1Re0 1z 可见只要原微分方程是稳定的（可见只要原微分方程是稳定的（），应用梯形公式进行数值计算时），应用梯形公式进行数值计算时）。因此梯形公式是恒稳公式。算法的稳定域如图）。因此梯形公式是恒稳公式。算法的稳定域如图412所示所示必然稳定（必然稳定（。梯形公式也是一种隐式公式

23、，所以隐式算法是一种恒稳算法。梯形公式也是一种隐式公式，所以隐式算法是一种恒稳算法。0-1-2Re()hIm()h不稳定域124.2.2 数值解法的选择原则数值解法的选择原则一般来说，选择数值计算方法从以下原则考虑：一般来说，选择数值计算方法从以下原则考虑：1、精度问题、精度问题 数值计算的精度主要受下面三类误差影响。截断误差、数值计算的精度主要受下面三类误差影响。截断误差、舍入误差和积累误差。舍入误差和积累误差。其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，

24、般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积累误差是上述两类误差积累的结果，它字长越长，舍入误差越小。积累误差是上述两类误差积累的结果，它同积分时间长短有关。一般积分步长越小，则积累误差越大（在一定同积分时间长短有关。一般积分步长越小，则积累误差越大（在一定的积分时间下）。所以，在一定的数值计算方法下，若从总误差考虑，的积分时间下）。所以，在一定的数值计算方法下，若从总误差考虑，必定有一个最佳步长值。必定有一个最佳步长值。2、计算速度问题、计算速度问题 计算速度决定于计算的步数以及每一步积分所需计算速度决定于计算的步数以及每一步积分所需的时间，而每一

25、步的计算时间同数值计算方法有关。例如四阶龙格的时间，而每一步的计算时间同数值计算方法有关。例如四阶龙格库塔法每一步需要求四个斜率值，花费较多的时间（但精度高）。在库塔法每一步需要求四个斜率值，花费较多的时间（但精度高）。在积分方法确定时，应在保证一定的计算精度下尽量能选用较大的步长积分方法确定时，应在保证一定的计算精度下尽量能选用较大的步长（必须满足数值稳定的条件），以缩短积分时间（有时往往选用变步（必须满足数值稳定的条件），以缩短积分时间（有时往往选用变步长的算法）。长的算法）。3、数值解的稳定性、数值解的稳定性 数值算法实际上就是将微分方程化成差分方程数值算法实际上就是将微分方程化成差分方

26、程进行求解，对一个稳定的微分方程组，经过变换得到的差分方程不一进行求解，对一个稳定的微分方程组，经过变换得到的差分方程不一定稳定。不同的数值计算方法有不同的稳定域。从稳定性角度看，隐定稳定。不同的数值计算方法有不同的稳定域。从稳定性角度看，隐式比显式好。式比显式好。4.3数值算法中的数值算法中的“病态病态”问题问题4.3.1 “病态病态”微分方程微分方程例例43，已知系统的状态方程为，已知系统的状态方程为 (0)101TXAXX其中其中 211920192120404040A，试用，试用MATLAB仿真仿真 解：系统状态方程可写为解：系统状态方程可写为：11232123312321192019

27、2120404040 xxxxxxxxxxxx 1x2x3x可建立系统的可建立系统的SIMULINK仿真模型如图仿真模型如图413所示。其中第一个积分器输出为所示。其中第一个积分器输出为，第二个积分器输出为，第二个积分器输出为 ，第三个积分器输出为，第三个积分器输出为。采用四阶龙格库塔法，选取固定步长的仿真方法，取采用四阶龙格库塔法，选取固定步长的仿真方法，取h=0.01， 参数设置参数设置 由图可见，由图可见，t=0.1s之后，曲线变化趋于平缓，若仍以之后，曲线变化趋于平缓，若仍以h=0.01计算，会使计算时间计算，会使计算时间拖得很长。拖得很长。 由结果可看出，误差很大，仿真有较大的失真。

28、由结果可看出，误差很大，仿真有较大的失真。 因此，受数值稳定性限制，只能取小步长进行仿真，若系统是一个慢变过程，因此，受数值稳定性限制，只能取小步长进行仿真，若系统是一个慢变过程，则计算速度大受影响。究其原因是状态方程的系数矩阵则计算速度大受影响。究其原因是状态方程的系数矩阵A的特征值差异较大。的特征值差异较大。可求出系统的特征值如下：可求出系统的特征值如下：A=-21 19 -20;19 -21 20;40 -40 -40A = -21 19 -20 19 -21 20 40 -40 -40 eig(A)ans = -2.0000 -40.0000 +40.0000i-40.0000 -40

29、.0000i三个特征值为：三个特征值为： 1232,40(1),40(1)jj 21Re()10Re()ii而系统的特征值而系统的特征值在实际系统中反映了动态过渡过程在实际系统中反映了动态过渡过程的作用。的作用。 各瞬态分量时间常数各瞬态分量时间常数“病态（病态（stiff）”方程。表述如下：方程。表述如下：一般线性常微分方程组：一般线性常微分方程组：0,(0)XAXBuXX的系数矩阵的系数矩阵A的特征值的特征值i具有如下特征具有如下特征 11()0max Re()min Re()iiii ni nRe （413）则称式（则称式（413）为病态方程，相应的系统称为病态系统。）为病态方程，相应的

30、系统称为病态系统。4.3.2 “病态病态”系统的仿真方法系统的仿真方法采用自动变步长数值计算方法采用自动变步长数值计算方法 对于例对于例43， (0)101TXAXX4.4 连续系统状态方程的离散化连续系统状态方程的离散化连续定常系统状态方程为：连续定常系统状态方程为：XAXBUYCXDX0(0)XX连续系统状态方程解的一般形式为：连续系统状态方程解的一般形式为： 00()0tt tA ttX teX teUdABTktkTt1,0若上式中令若上式中令 Tkt1采样时刻的状态采样时刻的状态 (1)(1)(1) ()( )kTATA kTkTX kTeX kTeBUd令令)0(,TddkT ()

31、0(1) ()()TATA TX kTeX kTeBU kTdkTTk)1(kT采用零阶保持器，采用零阶保持器，到到期间保持器的输出为恒定值。且等于期间保持器的输出为恒定值。且等于时刻的采样值时刻的采样值,()()()hU kTUkTU kT()0(1) ()()TATA TX kTeX kTed BU kT并记并记: )ATTA(T0GeHed B令令ddtTt,dtedtedeTAtTtTT000)(AA最后得最后得:0ATTAtGeHe dtB于是离散化的状态空间方程为：于是离散化的状态空间方程为：(1) ()()()()()X kTGX kTHU kTY kTCX kTDU kTc2d

32、( )函数可以立即得出连续系统离散化的模型来，该函数的调用命令为：函数可以立即得出连续系统离散化的模型来，该函数的调用命令为：,2 ( , , )GHc d A B T例例44 已知连续系统的状态方程为：已知连续系统的状态方程为：211701020213124xxuyx ，用，用MATLAB求求T=0.1时，相应的离散状态方程。时，相应的离散状态方程。在在MATLAB命令窗口输入：命令窗口输入： A=2 -1 -1;0 -1 0;0 2 1; B=7;2;3; C=1 2 4; D=0; G,H=c2d(A,B,0.1)G = 1.2214 -0.1162 -0.1162 0 0.9048 0

33、 0 0.2003 1.1052H = 0.7473 0.19030.3355 则离散状态方程为：1.22140.11620.11620.7473(1)00.90480( )0.1903( )00.20031.10520.3355X kX kU kMATLAB还提供了还提供了c2dm( )函数来作类似的变换，其调用格式为：函数来作类似的变换，其调用格式为：G,H,Cd,Dd=c2dm(A,B,C,D,T,选项选项)它与它与c2d函数的区别在于，它可以允许用户自己选择变换方法。函数的区别在于，它可以允许用户自己选择变换方法。下面采用两种不同的方法进行离散化：下面采用两种不同的方法进行离散化： 第

34、一种：第一种： G,H,Cd,Dd=c2dm(A,B,C,D,0.1,zoh)G = 1.2214 -0.1162 -0.1162 0 0.9048 0 0 0.2003 1.1052H = 0.7473 0.1903 0.3355Cd = 1 2 4Dd = 0第二种第二种G,H,Cd,Dd=c2dm(A,B,C,D,0.1,tustin)G = 1.2222 -0.1170 -0.1170 0 0.9048 0 0 0.2005 1.1053H = 7.4854 1.9048 3.3584Cd = 0.1111 0.2247 0.4152Dd =1.2364可看出用零阶保持器的离散化得出的

35、结果和直接调用可看出用零阶保持器的离散化得出的结果和直接调用c2d ( )函数所得的结果相同。函数所得的结果相同。 MATLAB的控制工具箱还给出了离散状态方程转化为连续状态方程的函的控制工具箱还给出了离散状态方程转化为连续状态方程的函数数d2c( )，其调用格式为：，其调用格式为： A,B=d2c(G,H,T)如对上面采用如对上面采用c2d( )得到的结果，采用如下的变换可得到原来的连续状态方程得到的结果，采用如下的变换可得到原来的连续状态方程A,B=d2c(G,H,0.1)A = 2.0000 -1.0000 -1.0000 0 -1.0000 0 0 2.0000 1.0000B = 7

36、.0000 2.0000 3.0000除了除了d2c( ) 函数外，函数外，MATLAB还提供了还提供了d2cm( )函数，其调用格式为：函数，其调用格式为： A, B, C, D=d2cm(G, H, Cd, Dd, T,选项选项)其中转换方法的选项意义如表其中转换方法的选项意义如表42所述。所述。对上面采用对上面采用Tustin变换得到的结果，采用如下的变换可得到原来的连变换得到的结果，采用如下的变换可得到原来的连续状态方程续状态方程A,B,C,D=d2cm(G,H,Cd,Dd,0.1,tustin)A = 2.0000 -1.0000 -1.0000 0 -1.0000 0 0 2.0000 1.0000B = 7.0000 2.0000 3.0000C = 1.0000 2.0000 4.0000D =-2.2204e-016由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。 下面介绍一种间接