连续系统的离散化方法

第四章 连续系统的离散化方法4.1 常微分方程的数值解法4.1.1 数值求解的基本概念已知一个一阶微分方程 0 0( , )( )x f t xx t x 0 x1 0 2 1 1, , , n nt t h t t h t t h 1 2, , , nx x x1 0t t h 0( )x t h应用数值求解的思路是，从初值开始，在一系列时刻，求未知解的近似解，h是计算步长。若x对t的各阶导数都存在。则x(t)在时的解用台劳级数表示为：2 (2) ( )0 0 0 0 0( ) ( ) ( ) ( ) ( )2! !k kh hx t h x t hx t x t x tk 4.1.2 欧拉法0 0 0, ( )t t x t x 1 0t t h 1 1( )x t x初始条件为计算时的， 欧拉法的计算公式为 1 0 0 0( , )x x hf t x 其一般公式为 1 ( , )k k k kx x hf t x 0 t f 0f 1f c 0t 1t2 3 (2) ( 1)0 0 0 0 0 0( , ) ( , ) ( , )2! 3! !k kn h h hR f t x f t x f t xk 称为截断误差 例41 用欧拉法求下述微分方程的数值解。2(0) 1x xx 解：因欧拉法的递推公式为 1 ( , )k k k kx x hf t x 所以 2( , )k k kf t x x则递推公式为： 21k k kx x hx 若取步长 0.1h由t=0开始计算可得 21 22 23 10 1 0.1 1 0.90.9 0.1 0.9 0.8190.819 0.1 0.819 0.75190.4627810 xxxx 上式的精确解是 11x t 数值计算与精确解的比较见表 t 0 0.1 0.2 0.3 1.0精确解x(t) 1 0.9090909 0.8333333 0.7692307 0.51 0.9 0.819 0.5719 0.4627810 4.1.3 龙格库塔法2 (2)0 0 0 0( ) ( ) ( ) ( )2!hx t h x t hx t x t 1 0 02 0 1 0 2 11 0 1 1 2 2( , )( , )( )K f t xK f t bh x b K hx x h a K a K 2 1 0 0 0 0 0( , ) ( ),2h f fx x hf t x f t t x xt x 2K 1 0 0( , )K f t x将在点处作台劳级数展开，并取线性部分可得 2 0 0 1 2 10 0( , ) f fK f t x b h b K ht t x xt x 将 1 0 1 0 0 2 0 0 1 2 10 0( , ) ( , ) f fx x a hf t x a h f t x bh b hKt t x xt x 1K 2K代入式 比较 各项系数得 1 2 2 1 2 21 11, ,2 2a a a b a b 待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其参数。一般有以下几种取法：1、 1 2 1 2 10, 1, 2a a b b 则1 0 21 0 02 0 0 1( , )( , )2 2x x hKK f t xh hK f t x K 一般形式 1 212 1( , )( , )2 2k k k kk kx x hKK f t xh hK f t x K 2、1 2 1 21 3 2, ,4 4 3a a b b 则1 0 1 21 0 02 0 0 1( 3 )4( , )2 2( , )3 3hx x K KK f t xK f t h x hK 一般形式 1 1 212 1( 3 )4( , )2 2( , )3 3k k k kk khx x K KK f t xK f t h x hK 3、 1 2 1 21, 12a a b b 则 1 0 1 21 0 02 0 0 1( )2( , )( , )hx x K KK f t xK f t h x hK 一般形式 1 1 212 1( )2( , )( , )k k k kk khx x K KK f t xK f t h x hK 以上几种递推公式均称为二阶龙格库塔公式。是较典型的几个常用算法。其中的方法3又称为预估校正法，或梯形法。 其意义如下：用欧拉法以斜率先求取一点， 再由此点求得另一斜率1 0 1x x hK 2 1 1 0 0 1( , ) ( , )K f t x f t h x hK 0 x 1K2K 1 22K KK 1 0 0 1 2( )2hx x hK x K K 然后，从点开始，既不按该点斜率变化，也不按预估点斜率变化，而是取两者平均值求得校正点，即： 。0 tf 0f 1f0t 1t1 0 0 1( )2hx x f f 四阶龙格库塔法的计算公式为：1 1 2 3 4( 2 2 )6k k hx x K K K K 12 13 24 3( , )( , )2 2( , )2 2( , )k kk kk kk kK f t x h hK f t x Kh hK f t x KK f t h x hK 对于用状态方程表示的高阶线性系统 X AX BUY CX 0(0)X X其中状态变量为 1 2, , nX x x x 用四阶龙格库塔法时，有 ( , )f t X AX BU 计算公式为： 1 1 2 3 4( 2 2 )6k k hX X K K K K 其中，12 1 13 2 24 3 3( , ) ( )( , ) ( ) ( )2 2 2 2( , ) ( ) ( )2 2 2 2( , ) ( ) ( )k k k kk k k kk k k kk k k kK f t X AX BU th h h hK f t X K A X K BU th h h hK f t X K A X K BU tK f t h X hK A X hK BU t h 相应的输出为 1 1k kY CX 按上式，取0,1,2, ,k n 不断递推， 即可求得所需时刻 1 1, , , nt t t ( )kX t()kYt的状态变量和输出值h iK德国学者Felhberg对传统的龙格库塔法提出了改进，在每一个计算步长内对f( )函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的步长为，定义下面6个变量11( , ), 1,2, ,6; 1,2, 1ii k ij j k ijK hf x K t h i j i 下一步的状态变量可由下式求出：61 1k k i iix x K 四阶/五阶龙格库塔法系数表i ji*i 0 16/135 25/2161/4 1/4 0 03/8 3/32 9/32 6656/12825 1408/2565 12/13 1932/2197 -7200/2197 7296/2197 28561/56430 2197/41041 439/216 -8 3680/513 -845/4104 -9/50 -1/51/2 -8/27 2 -3544/2565 1859/4104 -11/40 2/55 0这一方法又称为四阶/五阶龙格库塔法。 4.1.4 微分方程数值解的MATLAB实现1、 , 23( , , 0, )t x ode Fun Tspan x options , 45( , , 0, )t x ode Fun Tspan x options该指令适用于一阶常微分方程组 ( , ), 1,2,i i ix f t x i n 如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方程方可使用。2、输入参数 Fun为定义微分方程组 ( , ) , 1, 2 ,i i ix f t x i n M函数文件名，可以在文件名加写，或用英文格式单引号界定文件名。3、在编辑调试窗口中编写一阶常微分方程组 ( , ), 1,2,i i ix f t x i n 的M函数文件时，每个微分方程的格式必须与 ( , )i i ix f t x, it xix 1,2,i n 一致，即等号严格以“先自变量t，后函数”的固定顺序输入，表示微分方程的序数。左边为带求函数的一阶导数，右边函数的变量4、输入参数“Tspan”规定了常微分方程的自变量取值范围，它以矩阵t0,tf的形式输入，表示自变量 0,t t tf5、输入参数x0表示初始条件向量， 0 ( 0) ( 0) ( 0)x x t x t x t 微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解所必须的条件。6、输入参数options表示选项参数（包括tol，trace），可缺省，即取默认值，tol是控制结果精度的选项对ode23( )函数取 ，对ode45( )函数取 。trace为输出形式控制变量，如果trace不为0，则 会将仿真中间结果逐步地由频幕显示出来，否则将不 显示中间结果610 310 it it ( )j ix g t7、输出参数t,x为微分方程组解函数的列表（t和x都是列矩阵），它包含向量t各节点 和与对应向量x的第j个分量值（即第j个方程解），i表示节点序列数。( )x g t( )( )jjx g t8、输出参数t,x缺省时，输出解函数的曲线，即函数及其各的曲线。阶导数求解微分方程的指令还有ode113（多步解法器），ode15s（基于数字微分公式的解法器），ode23s（单步解法器），ode23T（梯形规则的一种自由插值实现），ode23TB（二阶隐式龙格库塔公式）等。 例 41 求解常微分方程 22 0tx x ，初始条件为： (0) 1; (0) 1; 0, 10 x x t 解：方法1 把二阶微分方程化成两个一阶微分方程组： 令 1 2,x xx x 则： 1 2 22 1 2x x tx x 首先编制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 写出系统的状态方程 20 1 0 (2 )1 0 1 tx x 首先编制M文件，并且函数名和M文件名相同。function xdot=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 -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 x x x x (0) 1; (0) 1x x 0, 20t例 42 求Var der Pol微分方程，在初始条件下的数值。解解：将二阶微分方程变换成一阶微分方程组 1 2,x x x x 则 1 2 22 1 2 12*(1 )x xx x x x 编制M文件，并且函数名和M文件名相同。function xdot=vdpl(t,x)xdot=zeros(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.0489 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所示。设置系统参数如下，将积分器初值设置为系统要求的初值，如图46所示。 XY示波器参数设置如图47所示。 仿真时间参数设置如图48所示。 系统仿真曲线如图49和410所示。 4.2数值解法的稳定性及选择原则4.2.1 数值解法的稳定性1.欧拉法对 0(0)x xx x 按欧拉公式计算得 1 (1 ) , 0,1, , 1k k k kx x h x h x k n 这是一个一阶差分方程，Z变换后得 ( ) (1 ) ( )zX z h X z ，其特征值为 1z h 1 1z h 1 1z h 2h 则该算法的稳定条件为2h 。算法的稳定域如图411所示。 0-1-2 Re( )h Im( )h不稳定域 2.梯形法1 1 1( ) ( )2 2k k k k k k kh hx x f f x x x ( )(1 ) (1 ) ( )2 2h hzX z X z 1 21 2hz h 1 21 2hz h 1Re 0 1z 可见只要原微分方程是稳定的（），应用梯形公式进行数值计算时）。因此梯形公式是恒稳公式。算法的稳定域如图412所示必然稳定（。梯形公式也是一种隐式公式，所以隐式算法是一种恒稳算法。 0-1-2 Re( )h Im( )h 不稳定域 12 4.2.2 数值解法的选择原则一般来说，选择数值计算方法从以下原则考虑：1、精度问题 数值计算的精度主要受下面三类误差影响。截断误差、舍入误差和积累误差。其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积累误差是上述两类误差积累的结果，它同积分时间长短有关。一般积分步长越小，则积累误差越大（在一定的积分时间下）。所以，在一定的数值计算方法下，若从总误差考虑，必定有一个最佳步长值。2、计算速度问题 计算速度决定于计算的步数以及每一步积分所需的时间，而每一步的计算时间同数值计算方法有关。例如四阶龙格库塔法每一步需要求四个斜率值，花费较多的时间（但精度高）。在积分方法确定时，应在保证一定的计算精度下尽量能选用较大的步长（必须满足数值稳定的条件），以缩短积分时间（有时往往选用变步长的算法）。3、数值解的稳定性 数值算法实际上就是将微分方程化成差分方程 进行求解，对一个稳定的微分方程组，经过变换得到的差分方程不一定稳定。不同的数值计算方法有不同的稳定域。从稳定性角度看，隐式比显式好。 4.3数值算法中的“病态”问题4.3.1 “病态”微分方程例43，已知系统的状态方程为 (0) 1 0 1 TX AX X 其中 21 19 2019 21 2040 40 40A ，试用MATLAB仿真 解：系统状态方程可写为： 1 1 2 32 1 2 33 1 2 321 19 2019 21 2040 40 40 x x x xx x x xx x x x 1x 2x 3x可建立系统的SIMULINK仿真模型如图413所示。其中第一个积分器输出为，第二个积分器输出为，第三个积分器输出为。 采用四阶龙格库塔法，选取固定步长的仿真方法，取h=0.01， 参数设置 由图可见，t=0.1s之后，曲线变化趋于平缓，若仍以h=0.01计算，会使计算时间拖得很长。 由结果可看出，误差很大，仿真有较大的失真。 因此，受数值稳定性限制，只能取小步长进行仿真，若系统是一个慢变过程，则计算速度大受影响。究其原因是状态方程的系数矩阵A的特征值差异较大。可求出系统的特征值如下：A=-21 19 -20;19 -21 20;40 -40 -40 A = -21 19 -20 19 -21 20 40 -40 -40 eig(A)ans = -2.0000 -40.0000 +40.0000i-40.0000 -40.0000i三个特征值为： 1 2 32, 40(1 ), 40(1 )j j 2 1Re( ) 10Re( ) ii而系统的特征值在实际系统中反映了动态过渡过程的作用。 各瞬态分量时间常数 “病态（stiff）”方程。表述如下：一般线性常微分方程组：0, (0)X AX Bu X X 的系数矩阵A的特征值i具有如下特征 11 ( ) 0max Re( ) min Re( )i i ii ni nRe （413）则称式（413）为病态方程，相应的系统称为病态系统。 4.3.2 “病态”系统的仿真方法采用自动变步长数值计算方法 对于例43， (0) 1 0 1 TX AX X 4.4 连续系统状态方程的离散化连续定常系统状态方程为： X AX BUY CX DX 0(0)X X连续系统状态方程解的一般形式为： 0 0 ( )0 tt t A ttX t e X t e U d A B TktkTt 1,0 若上式中令 Tkt 1采样时刻的状态 ( 1) ( 1) ( 1) ( ) ( )k TAT A k TkTX k T e X kT e BU d 令)0(, TddkT ( )0( 1) ( ) ( )TAT A TX k T e X kT e BU kT d kT Tk )1( kT采用零阶保持器，到期间保持器的输出为恒定值。且等于时刻的采样值, ( ) ( ) ( ) hU kT U kT U kT ( )0( 1) ( ) ( )TAT A TX k T e X kT e d BU kT 并记: )ATT A(T0G eH e d B 令 ddtTt , dtedtede T AtT tT T 000 )( AA 最后得: 0 ATT AtG eH e dtB 于是离散化的状态空间方程为：( 1) ( ) ( )( ) ( ) ( )X k T GX kT HU kTY kT CX kT DU kT c2d( )函数可以立即得出连续系统离散化的模型来，该函数的调用命令为： , 2 ( , , )G H c d A B T例44 已知连续系统的状态方程为： 2 1 1 70 1 0 20 2 1 31 2 4x x uy x ，用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 0 0.2003 1.1052H = 0.7473 0.19030.3355 则离散状态方程为： 1.2214 0.1162 0.1162 0.7473( 1) 0 0.9048 0 ( ) 0.1903 ( )0 0.2003 1.1052 0.3355X k X k U k MATLAB还提供了c2dm( )函数来作类似的变换，其调用格式为：G,H,Cd,Dd=c2dm(A,B,C,D,T,选项)它与c2d函数的区别在于，它可以允许用户自己选择变换方法。下面采用两种不同的方法进行离散化： 第一种： 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.3355 Cd = 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可看出用零阶保持器的离散化得出的结果和直接调用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.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.0000 D =-2.2204e-016由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。 下面介绍一种间接的变换方法，首先将要转换的连续传递函数转换成连续的状态方程模型，再对该连续状态方程离散化，然后将相应的离散状态方程再转换成传递函数，即得相应的脉冲传递函数。 可编制如下的程序，并命名为：tfc2d.m 可实现由连续传递函数表示的系统转换成脉冲传递函数表示的离散系统。function nd,dd=tfc2d(nc,dc,T) 定义函数其功能是将连续传递函数转换成脉冲传递函数A,B,C,D=tf2ss(nc,dc); 将连续传递函数转换成连续状态方程G,H=c2d(A,B,T); 将连续状态方程离散化得到离散状态方程nd,dd=ss2tf(G,H,C,D); 将离散状态方程转换为脉冲传递函数 可编制如下的程序，并命名为：tfd2c.m 可实现由脉冲传递函数表示的离散系统转换成连续传递函数表示的系统。function nc,dc=tfd2c(nd,dd,T) 定义函数其功能是将脉冲传递函数数转换成连续传递函A,B,C,D=tf2ss(nd,dd); 将脉冲传递函数转换成离散状态方程G,H=d2c(A,B,T); 将离散状态方程转换成连续状态方程nc,dc=ss2tf(G,H,C,D); 将连续状态方程转换为连续传递函数 例45 已知连续系统的传递函数为：，求其脉冲传递函数。23 212 86 3 2s ss s s 在MATLAB命令窗口出入：nc=1 12 8; dc=1 6 3 2; nd,dd=tfc2d(nc,dc,0.1)nd = 0 0.1255 -0.1546 0.0352dd = 1.0000 -2.5255 2.0758 -0.5488则离散脉冲传递函数为： 23 20.1255 0.1546 0.03522.5255 2.0758 0.5488s ss s s 采用下面的变换可得到原连续传递函数：nc,dc=tfd2c(nd,dd,0.1)nc = 0 1.0000 12.0000 8.0000dc = 1.0000 6.0000 3.0000 2.0000