燕山大学自动化专业,控制系统仿真及MATLAB语言,第四章 连续系统的离散化方法,4.1,常微分方程的数值解法,一,.,数值求解的基本概念,设微分方程为,则求解方程中函数,x(t,),问题的常微分方程初值问题,所谓数值求解就是要在时间区间,a,b,中取若干离散点,求出微分方程在这些时刻的近似值,取前两项近似：,这种方法的几何意义就是把,f(t,x,),在区间,t,k,t,k+1,内的曲边面积用矩形面积近似代替。计算简单，计算量小，而且可以自启动。当,h,很小时，造成的误差是允许的。该算法具有一阶精度。,取,k=0,1,2,N,从,t,0,开始，逐点递推求解,t,1,时的,y,1,，,t,2,时的,y,2,直至,t,n,时的,y,n,，称之为欧拉递推公式。,矩形面积,1.,欧拉法,欧拉法的特点：,导出简单，几何意义明显，便于理解，能说明构造数值解法一般计算公式的基本思想,。,通常用它来说明有关的基本概念。,例,设系统方程为,用,Euler,法求其数值解（取步长,),递推公式为,则,已知方程的解析解为 精确解和解析解作比较：,误差在 数量级，精度较差。,t,0,0.1,0.2,0.3,0.4,0.5,1.0,精确解,1,0.909,0.833,0.769,0.666,0.625,0.5,数值解,1,0.9,0.819,0.752,0.659,0.647,0.463,2.,龙格库塔法,基本思想：,取,Taylor,级数展开式前三项近似求解，并利用线性组合代替导数的求解。,既可避免计算高阶导数，又可提高数值积分的精度，这就是,Runge-Kutta,法的基本思想。,2.,龙格库塔法,r,为精度阶次，,a,i,为待定系数，由精度确定；,k,i,用下式表示,线性组合,等各阶导数不易计算，用下式中,ki,的线性组合代替,1),当,r=1,时：,与,Taloy,展开式相比较，可得,a,1,=1,则上式成为,欧拉递推公式,2),当,r=2,时：,将 在点 展成,Taylor,级数,与台劳公式的二阶展开近似公式相比，可得以下关系：,三个方程，四个未知数，解不唯一,各个系数的几种取法,见书上。,3)r=4,时，四阶龙格库塔公式,-,最常用：,仿真中遇到的大多数工程实际问题，四阶龙格库塔法以能满足精度要求，其截断误差,o(h,5,),与,h,5,同数量级。该法可以自启动。,4),、状态空间四阶龙格,-,库塔递推式,若单输入单输出系统的状态空间表达式为：,在仿真中，对于,n,阶系统，状态方程可以写成一阶微分方程,根据四阶龙格,-,库塔公式,有,T=,t,k,时刻的,x,i,值,T=,t,k,+h,时刻的,x,i,值,另,状态方程的四阶龙格,-,库塔公式如下：,RK,法的特点：,1,需要存储的数据少，占用的存储空间少；,2,只需知道初值，即可启动递推公式进行计,算，可自启动；,3,容易实现变步长运算。,4,每积分一步需要计算多次右函数，计算量,大。,基于龙格库塔法，,MATLAB,提供了求常微分方程数值解,的函数，一般调用格式为：,t,x=ode23(xfun,t0,tf,x0),t,x=ode45(xfun,t0,tf,x0),常微分方,程函数名,起始,时间,终止,时间,初始状,态向量,输入,输出,4/5,阶龙格,-,库塔算法,2/3,阶龙格,-,库塔算法,3.,常微分方程,Matlab,求解,解,:,令,y,1,=x,，,y,2,=x,1,、建立,M-,文件,vdp.m,如下：,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,Y,(:,2),3,、结果,解,1,、建立,m-,文件,rigid.m,如下：,function,dy,=,rigid(t,y,),dy,=zeros(3,1);,dy(1)=y(2)*y(3);,dy(2)=-y(1)*y(3);,dy(3)=-0.51*y(1)*y(2);,2,、,取,t,0,=0,，,t,f,=12,，输入命令：,T,Y=ode45(rigid,0 12,0 1 1);,plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3,、结果如图,图中，,y,1,的图形为实线，,y,2,的图形为“*”线，,y,3,的图形为“,+”,线,.,4.2,数值算法的稳定性及求解原则,1,.,数值算法的稳定性,特征根在,s,平面的左半平面，系统稳定。,(1),欧拉法：,稳定：,(2),梯形法：,恒稳,2,.,数值算法的选择原则,Matlab,提供了微分方程数值求解的一般方法,作为仿真算法的使用者,可不必考虑算法具体实现,而应关心各种方法在使用中会出现的问题,以及如何在仿真中恰当的选用这些方法,.,一般,选用数值算法从以下几个方面考虑,:,(1),精度,受算法和,h,影响,截断误差,+,舍入误差,=,累计误差,(2),计算速度,受算法和,h,影响,算法简单，速度就快些。,(3),稳定性,受,h,影响，一般,h,（,2-3,）,系统最小时间,4.3,数值算法中的“病态”问题,1,“,病态”常微分方程,例：,其中,采用四阶龙格库塔法,h=0.01,时，计算时间长,h=0.04,时，误差很大,当,h0.05,后，曲线发散振荡，数值不稳定，完全失去意义,系统矩阵的特征值差异较大,一般线性常微分方程组：,的系数矩阵,A,的特征值具有如下特征：,则称为“病态”方程。,定义：,2,控制系统仿真中的,“,病态”问题,1,病态系统中绝对值最大的特征值对应于系统动态性能解中瞬态分量衰减最快的部分，它反映了系统的动态响应和系统的反应灵敏度。一般与系统中具有最小时间常数,T,min,的环节有关，要求计算步长,h,取得很小。,2,病态系统中绝对值最小的特征值对应于系统动态性能解中瞬态分量衰减最慢的部分，它决定了整个系统的动态过渡过程时间的长短。一般与系统中具有最小时间常数,T,max,的环节有关，要求计算步长,h,取得很大。,3,对于病态问题的仿真需要寻求更加合理的算法，以解决病态系统带来的选取计算步长与计算精度，计算时间之间的矛盾。,3,“,病态”系统的仿真方法,采用稳定性好，计算精度高的数值算法，并且允许计算步长能根据系统性能动态变化的情况在一定范围内作相应的变化，采用,隐式吉尔法,该法已经证明对病态方程求解过程是数值稳定的。,隐式吉尔法从理论上十分适应于病态系统，但需要解决好以下问题,(1),自启动,r,阶多步算式无法自启动，需要用单步法求出前,r,步值,(2),预估迭代 迭代方法要求收敛性良好，否则在大步长时会造成数 值发散。,(3),变步长 初始阶段采用小步长，随后可逐步放大步长。,对不同精度要求的系统仿真，要考虑变阶次问题，即为减小每一步计算的截断误差，以提高精度，应选用较高的阶次，而当精度较低时，为减少工作量，则应选取较低的阶次。仿真时应根据估计误差 与给定的误差精度相比较改变步长或阶次来重新计算。,4.4,连续系统状态方程的离散化,上章所述的连续系统数学模型的离散化，是通过数值积分法实现的，尽管面向结构图的仿真方法是按环节给定参数，但是在计算时还是,按整个系统进行离散化,，这就不便于引进非线性环节以进行非线性系统的仿真。在本节，将介绍连续系统离散模型的建立和仿真。,数值积分法叠代求解,h,改变时，叠代过程重复求解，费时繁琐,不能对非线性环节单独考虑。,连续系统离散化,思想：用差分方程描述连续系统的状态方程模型,（因为差分方程的主要特点就是方程中各变量由各相邻时刻的变,化量制约，这相当于递推方程）,1,、连续系统的离散化,设连续系统状态方程为,其中,为状态初始值,.,则由现代控制理论基础知，状态变量,X,(,t,),的解为,而,t,=(,k,+1),T,时，可表示为,当系统输入,u(t,),给定时，可求出系统离散化状态方程的解。一般，,u(t,),未知，通常采用两种方法近似处理：,(1),令,u(kT+t),u(kT,)(0,t,T),相当在系统输入端加一个采样开关和零阶保持器,X,(,k,+1),T,)=,GX,(,kT,)+,Hu,(,kT,),G,=,e,A T,，为,t,=,T,时的状态转移矩阵,离散后的状态空间表达式为：,、,Matlab,表示,已知连续系统状态方程为,在采样周期,T,下离散后的状态空间表达可表示为：,在,Matlab,中，若已知连续系统状态方程各阵模型参数,(,A,、,B,、,C,、,D,),以及采样周期,T,，则语句：,G,，,H,=c2d(,A,B,T,),返回的矩阵,G,、,H,就是所要求的,(,T,),、,m,(,T,),。,此外，,Matlab,还提供了功能更强的求取连续系统离散化矩阵函数,c2dm(),他容许调用时选用离散化变换方式，并且得到的是标准的离散化状态方程。,G,H,C,D,=c2dm(,A,B,C,D,T,选项,),表,离散化变换方式选项,选 项,说 明,Zoh,假设输入端加一个采样开关和零阶保持器,Foh,假设输入端加一个采样开关和一阶保持器。,Tustin,采用双线性变换（,Tustin,算法）方法,Prewarp,采用改进的,Tustin,变换方法,Matched,采用,SISO,系统的零极点匹配法,2.,离散函数的连续化,在,MATLAB,中也提供了从离散化系统转换为连续系统各系数矩阵求取的功能函数，其调用格式分别如下,A,B,=d2c(,G,H,T,),或,A,B,C,D,=d2cm(,G,H,C,D,T,选项,),其中选项同上。,例,对连续系统在采样周期,T=0.1,时进行离散化。,解：利用以下,Matlab,命令，对系统进行离散化,k=6;,Z=-3;,P=-1,；,-2,；,-5,；,T=0.1;,A,B,C,D=zp2ss(Z,P,k);,G1,H1=c2d(A,B,T);,G2,H2,C2,D2=c2dm(A,B,C,D,T,zoh);,G3,H3,C3,D3=c2dm(A,B,C,D,T,foh);,3,、脉冲传递函数求解,连续 离散,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);,离散 连续,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);,