《第四章连续系统的离散化方法精选文档.ppt》由会员分享，可在线阅读，更多相关《第四章连续系统的离散化方法精选文档.ppt（38页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第四章连续系统的离散化方法本讲稿第一页，共三十八页其一般公式为 称为截断误差称为截断误差 例例41 用欧拉法求下述微分方程的数值解。用欧拉法求下述微分方程的数值解。本讲稿第二页，共三十八页解：因欧拉法的递推公式为解：因欧拉法的递推公式为 所以所以 则递推公式为：则递推公式为：若取步长若取步长 由由t=0开始计算可得开始计算可得 上式的精确解是上式的精确解是 数值计算与精确解的比较见表数值计算与精确解的比较见表 t00.10.20.31.0精确解x(t)10.90909090.83333330.76923070.510.90.8190.57190.4627810本讲稿第三页，共三十八页4.1.3
2、 龙格库塔法龙格库塔法将将在点在点处处作台作台劳级劳级数展开，并取数展开，并取线线性部分可得性部分可得 本讲稿第四页，共三十八页将将 代入式代入式 比较比较 各项系数得各项系数得 待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其参数。一般有以下几种取法参数。一般有以下几种取法：1、则则一般形式一般形式 本讲稿第五页，共三十八页2、则则 一般形式一般形式 3、则则 一般形式一般形式 以上几种递推公式均称为二阶龙格库塔公式。是较典型的几以上几种递推公式均称为二阶龙格库塔公式。是较典型的几个常用算法。其中的方法个常用算法。其中的
3、方法3又称为预估校正法，或梯形法。又称为预估校正法，或梯形法。本讲稿第六页，共三十八页其意义如下：用欧拉法以斜率先求取一点，其意义如下：用欧拉法以斜率先求取一点，再由此点求得另一斜率再由此点求得另一斜率然后，从然后，从点开始，既不按该点斜率点开始，既不按该点斜率变化，也不按预估点斜率变化，也不按预估点斜率变化，而是取两者平均值变化，而是取两者平均值求得校正点，即：求得校正点，即：。本讲稿第七页，共三十八页四阶龙格库塔法的计算公式为：四阶龙格库塔法的计算公式为：对于用状态方程表示的高阶线性系统对于用状态方程表示的高阶线性系统 其中状态变量为其中状态变量为 用四阶龙格库塔法时，有用四阶龙格库塔法时
4、，有 计算公式为：计算公式为：本讲稿第八页，共三十八页其中，其中，相应的输出为相应的输出为 按上式，取按上式，取不断不断递递推，推，即可求得所需时刻即可求得所需时刻 的状态变量的状态变量和输出值和输出值德国学者德国学者Felhberg对传统的龙格库塔法提出了改进，在每一个计算步长内对传统的龙格库塔法提出了改进，在每一个计算步长内对对f()函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的步长为步长为，定义下面，定义下面6个个变量变量下一步的状态变量可由下式求出：下一步的状态变量可由下式求出：本讲稿第九页，共三十八页 四阶四阶
5、/五阶龙格库塔法系数表五阶龙格库塔法系数表016/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、该指令适用于一阶常微分方程组该指令适用
6、于一阶常微分方程组 如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方程方可使用。程方可使用。2、输输入参数入参数为为定定义义微分方程微分方程组组 M函数文件名，可以在文件名加写函数文件名，可以在文件名加写，或用英文格式单引号界定文件名。，或用英文格式单引号界定文件名。3、在编辑调试窗口中编写一阶常微分方程组在编辑调试窗口中编写一阶常微分方程组 的的M函数文件时，每个微分方程的格式必须与函数文件时，每个微分方程的格式必须与 一致，即等号一致，即等号严格以严格以“先自变量先自变量t，后函数，后函数”的固定顺序输入的固定
7、顺序输入，表示微分方程的序数。表示微分方程的序数。左边为带求函数的一阶导数，右边函数的变量左边为带求函数的一阶导数，右边函数的变量4、输入参数、输入参数“Tspan”规定了常微分方程的自变量取值范围，规定了常微分方程的自变量取值范围，它以矩阵它以矩阵t0,tf的形式输入，表示自变量的形式输入，表示自变量5、输入参数、输入参数x0表示初始条件向量，表示初始条件向量，本讲稿第十一页，共三十八页微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解所必须的条件。所必须的条件。6、输入参数、输入参数options表示选项参数（包括表示选
8、项参数（包括tol，trace），可缺省，即取默），可缺省，即取默认值，认值，tol是控制结果精度的选项对是控制结果精度的选项对ode23()函数取函数取 ，对，对ode45()函数函数取取 。trace为输出形式控制变量，如果为输出形式控制变量，如果trace不为不为0，则，则 会将仿真中间结果逐步地由频幕显示出来，否则将不会将仿真中间结果逐步地由频幕显示出来，否则将不 显示中间结果显示中间结果7、输出参数、输出参数t,x为微分方程组解函数的列表（为微分方程组解函数的列表（t和和x都是列矩阵），它都是列矩阵），它包含向量包含向量t各节点各节点 和与和与对应向量对应向量x的第的第j个分量值个分
9、量值（即第（即第j个方程解），个方程解），i表示节点序列数。表示节点序列数。8、输出参数、输出参数t,x缺省时，输出解函数的曲线，即函数缺省时，输出解函数的曲线，即函数及其各及其各的曲线。的曲线。阶导数阶导数求解微分方程的指令还有求解微分方程的指令还有ode113（多步解法器），（多步解法器），ode15s（基于数字（基于数字微分公式的解法器），微分公式的解法器），ode23s（单步解法器），（单步解法器），ode23T（梯形规则的（梯形规则的一种自由插值实现），一种自由插值实现），ode23TB（二阶隐式龙格库塔公式）等。（二阶隐式龙格库塔公式）等。本讲稿第十二页，共三十八页例例 41 求解
10、常微分方程求解常微分方程，初始条件为初始条件为：解：方法解：方法1 把二阶微分方程化成两个一阶微分方程组：把二阶微分方程化成两个一阶微分方程组：令令 则：则：首先编制首先编制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 写出系统的状态
11、方程写出系统的状态方程本讲稿第十三页，共三十八页首先编制首先编制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
12、 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 -本讲稿第十四页，共三十八页例例 42 求求Var der Pol微分方程微分方程，在初始条件，在初始条件下的数值下的数值。解解解：将二阶微分方程变换成一阶微分方程组解：将二阶微分方程变换成一阶微分方程组 则则 本讲稿第十五页，共三十八页编制编制M文件，并且函数名和文件，并且函数
13、名和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.53
14、89 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进行仿真。首先建立进行仿真
15、。首先建立SIMULINK仿真模型仿真模型如图如图45所示。所示。设置系统参数如下，将积分器初值设置为系统要求的初值，如图设置系统参数如下，将积分器初值设置为系统要求的初值，如图46所示。所示。本讲稿第十七页，共三十八页XY示波器参数设置如图示波器参数设置如图47所示。所示。本讲稿第十八页，共三十八页仿真时间参数设置如图仿真时间参数设置如图48所示。所示。系统仿真曲线如图系统仿真曲线如图49和和410所示。所示。本讲稿第十九页，共三十八页4.2数值解法的稳定性及选择原则数值解法的稳定性及选择原则4.2.1 数值解法的稳定性数值解法的稳定性1.欧拉法欧拉法对对 按欧拉公式计算得按欧拉公式计算得
16、这是一个一阶差分方程，这是一个一阶差分方程，Z变换后得变换后得，其特征值为，其特征值为 则该则该算法的算法的稳稳定条件定条件为为。算法的。算法的稳稳定域如定域如图图411所示。所示。本讲稿第二十页，共三十八页2.梯形法梯形法可见只要原微分方程是稳定的（可见只要原微分方程是稳定的（），应用梯形公式进行数值计算时），应用梯形公式进行数值计算时）。因此梯形公式是恒稳公式。算法的稳定域如图）。因此梯形公式是恒稳公式。算法的稳定域如图412所示所示必然稳定（必然稳定（。梯形公式也是一种隐式公式，所以隐式算法是一种恒稳算法。梯形公式也是一种隐式公式，所以隐式算法是一种恒稳算法。本讲稿第二十一页，共三十八页
17、4.2.2 数值解法的选择原则数值解法的选择原则一般来说，选择数值计算方法从以下原则考虑：一般来说，选择数值计算方法从以下原则考虑：1、精度问题、精度问题 数值计算的精度主要受下面三类误差影响。截断误差、舍数值计算的精度主要受下面三类误差影响。截断误差、舍入误差和积累误差。入误差和积累误差。其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一般减小步长可减其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积累误差是上述
18、两类误差积累的结果，它同积分时间长短有关。一般积分步长越小，则累误差是上述两类误差积累的结果，它同积分时间长短有关。一般积分步长越小，则积累误差越大（在一定的积分时间下）。所以，在一定的数值计算方法下，若从总误积累误差越大（在一定的积分时间下）。所以，在一定的数值计算方法下，若从总误差考虑，必定有一个最佳步长值。差考虑，必定有一个最佳步长值。2、计算速度问题、计算速度问题 计算速度决定于计算的步数以及每一步积分所需的时间，计算速度决定于计算的步数以及每一步积分所需的时间，而每一步的计算时间同数值计算方法有关。例如四阶龙格库塔法每一步需而每一步的计算时间同数值计算方法有关。例如四阶龙格库塔法每一
19、步需要求四个斜率值，花费较多的时间（但精度高）。在积分方法确定时，应在要求四个斜率值，花费较多的时间（但精度高）。在积分方法确定时，应在保证一定的计算精度下尽量能选用较大的步长（必须满足数值稳定的条件），保证一定的计算精度下尽量能选用较大的步长（必须满足数值稳定的条件），以缩短积分时间（有时往往选用变步长的算法）。以缩短积分时间（有时往往选用变步长的算法）。3、数值解的稳定性、数值解的稳定性 数值算法实际上就是将微分方程化成差分方程进行求解，对数值算法实际上就是将微分方程化成差分方程进行求解，对一个稳定的微分方程组，经过变换得到的差分方程不一定稳定。不同的数值计算一个稳定的微分方程组，经过变换
20、得到的差分方程不一定稳定。不同的数值计算方法有不同的稳定域。从稳定性角度看，隐式比显式好。方法有不同的稳定域。从稳定性角度看，隐式比显式好。本讲稿第二十二页，共三十八页4.3数值算法中的数值算法中的“病态病态”问题问题4.3.1“病态病态”微分方程微分方程例例43，已知系统的状态方程为，已知系统的状态方程为 其中其中 ，试试用用MATLAB仿真仿真 解：系统状态方程可写为解：系统状态方程可写为：可建立系统的可建立系统的SIMULINK仿真模型如图仿真模型如图413所示。其中第一个积分器输出为所示。其中第一个积分器输出为，第二个积分器输出为，第二个积分器输出为，第三个积分器输出为，第三个积分器输
21、出为。本讲稿第二十三页，共三十八页本讲稿第二十四页，共三十八页采用四阶龙格库塔法，选取固定步长的仿真方法，取采用四阶龙格库塔法，选取固定步长的仿真方法，取h=0.01，参数设置参数设置 由图可见，由图可见，t=0.1s之后，曲线变化趋于平缓，若仍以之后，曲线变化趋于平缓，若仍以h=0.01计算，会使计算时间计算，会使计算时间拖得很长。拖得很长。本讲稿第二十五页，共三十八页由结果可看出，误差很大，仿真有较大的失真。由结果可看出，误差很大，仿真有较大的失真。本讲稿第二十六页，共三十八页因此，受数值稳定性限制，只能取小步长进行仿真，若系统是一个慢变过程，因此，受数值稳定性限制，只能取小步长进行仿真，
22、若系统是一个慢变过程，则计算速度大受影响。究其原因是状态方程的系数矩阵则计算速度大受影响。究其原因是状态方程的系数矩阵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三个特征值为：三个特征值为：而系而系统统的特征的特征值值在在实际实际系系统统中反映了中反映了动态过动态过渡渡过过程程的作用。的作用。各瞬态
23、分量时间常数各瞬态分量时间常数本讲稿第二十八页，共三十八页“病态（病态（stiff）”方程。表述如下：方程。表述如下：一般线性常微分方程组：一般线性常微分方程组：的系数矩的系数矩阵阵A的特征的特征值值具有如下特征具有如下特征（413）则称式（则称式（413）为病态方程，相应的系统称为病态系统。）为病态方程，相应的系统称为病态系统。本讲稿第二十九页，共三十八页4.3.2“病态病态”系统的仿真方法系统的仿真方法采用自动变步长数值计算方法采用自动变步长数值计算方法 对于例对于例43，本讲稿第三十页，共三十八页4.4 连续系统状态方程的离散化连续系统状态方程的离散化连续定常系统状态方程为：连续定常系统
24、状态方程为：连续系统状态方程解的一般形式为：连续系统状态方程解的一般形式为：若上式中令若上式中令 采样时刻的状态采样时刻的状态 令令 采用零阶保持器，采用零阶保持器，到到期间保持器的输出为恒定值。且等于期间保持器的输出为恒定值。且等于时刻的采样值时刻的采样值,并记并记:令令本讲稿第三十一页，共三十八页最后得最后得:于是离散化的状态空间方程为：于是离散化的状态空间方程为：c2d()函数可以立即得出连续系统离散化的模型来，该函数的调用命令为：函数可以立即得出连续系统离散化的模型来，该函数的调用命令为：例例44 已知连续系统的状态方程为：已知连续系统的状态方程为：，用，用MATLAB求求T=0.1时
25、，相应的离散状态方程。时，相应的离散状态方程。本讲稿第三十二页，共三十八页在在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 则离散状态方程为：本讲稿第三十三页，共三十八页MATLAB还提供了还提供了c2dm()函数来作类似的变换，其调用格式为：函数来作类似的变换，其调用格式为：G,H,Cd,Dd=c2dm(A,B,C,D,T,选项选项)它与它与c
26、2d函数的区别在于，它可以允许用户自己选择变换方法。函数的区别在于，它可以允许用户自己选择变换方法。下面采用两种不同的方法进行离散化：下面采用两种不同的方法进行离散化：第一种：第一种：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
27、=7.4854 1.9048 3.3584Cd=0.1111 0.2247 0.4152Dd=1.2364可看出用零阶保持器的离散化得出的结果和直接调用可看出用零阶保持器的离散化得出的结果和直接调用c2d()函数所得的结果相同。函数所得的结果相同。本讲稿第三十四页，共三十八页MATLAB的控制工具箱还给出了离散状态方程转化为连续状态方程的函数的控制工具箱还给出了离散状态方程转化为连续状态方程的函数d2c()，其调用格式为：，其调用格式为：A,B=d2c(G,H,T)如对上面采用如对上面采用c2d()得到的结果，采用如下的变换可得到原来的连续状态方程得到的结果，采用如下的变换可得到原来的连续状态
28、方程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
29、,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由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。本讲稿第三十六页，共三十八页下面介绍一种间接的变换方法，首先将要转换的连续传递函数转换成连续的状态方下面介绍一种间接的变换方法，首先将要转换的连续传递函数转换成连续的
30、状态方程模型，再对该连续状态方程离散化，然后将相应的离散状态方程再转换成传递函程模型，再对该连续状态方程离散化，然后将相应的离散状态方程再转换成传递函数，即得相应的脉冲传递函数。数，即得相应的脉冲传递函数。可编制如下的程序，并命名为：可编制如下的程序，并命名为：tfc2d.m 可实现由连续传递函数表示的系统转换可实现由连续传递函数表示的系统转换成脉冲传递函数表示的离散系统。成脉冲传递函数表示的离散系统。function nd,dd=tfc2d(nc,dc,T)定义函数其功能是将连续传递函数转换成脉冲传递定义函数其功能是将连续传递函数转换成脉冲传递函数函数A,B,C,D=tf2ss(nc,dc)
31、;将连续传递函数转换成连续状态方程将连续传递函数转换成连续状态方程G,H=c2d(A,B,T);将连续状态方程离散化得到离散状态方程将连续状态方程离散化得到离散状态方程nd,dd=ss2tf(G,H,C,D);将离散状态方程转换为脉冲传递函数将离散状态方程转换为脉冲传递函数 可编制如下的程序，并命名为：可编制如下的程序，并命名为：tfd2c.m 可实现由脉冲传递函数表示的离散系可实现由脉冲传递函数表示的离散系统转换成连续传递函数表示的系统。统转换成连续传递函数表示的系统。function nc,dc=tfd2c(nd,dd,T)定义函数其功能是将脉冲传递函数数转换成定义函数其功能是将脉冲传递函
32、数数转换成连续传递函连续传递函A,B,C,D=tf2ss(nd,dd);将脉冲传递函数转换成离散状态方程将脉冲传递函数转换成离散状态方程G,H=d2c(A,B,T);将离散状态方程转换成连续状态方程将离散状态方程转换成连续状态方程nc,dc=ss2tf(G,H,C,D);将连续状态方程转换为连续传递函数将连续状态方程转换为连续传递函数本讲稿第三十七页，共三十八页例例45 已知连续系统的传递函数为：已知连续系统的传递函数为：，求其脉冲传递函数。，求其脉冲传递函数。在在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则则离散脉冲离散脉冲传递传递函数函数为为：采用下面的变换可得到原连续传递函数：采用下面的变换可得到原连续传递函数：nc,dc=tfd2c(nd,dd,0.1)nc=0 1.0000 12.0000 8.0000dc=1.0000 6.0000 3.0000 2.0000本讲稿第三十八页，共三十八页