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

第四章

连续系统的离散化方法4.1常微分方程的数值解法

一.数值求解的基本概念设微分方程为则求解方程中函数x(t)问题的常微分方程初值问题所谓数值求解就是要在时间区间[a,b]中取若干离散点

求出微分方程在这些时刻的近似值取前两项近似：这种方法的几何意义就是把f(t,x)在区间[tk,tk+1]内的曲边面积用矩形面积近似代替。计算简单，计算量小，而且可以自启动。当h很小时，造成的误差是允许的。该算法具有一阶精度。取k=0,1,2,…N,从t0开始，逐点递推求解t1时的y1，t2时的y2…,直至tn时的yn，称之为欧拉递推公式。矩形面积1.欧拉法

欧拉法的特点：导出简单，几何意义明显，便于理解，能说明构造数值解法一般计算公式的基本思想。通常用它来说明有关的基本概念。例设系统方程为

用Euler法求其数值解（取步长,)递推公式为则已知方程的解析解为

精确解和解析解作比较：

误差在数量级，精度较差。t00.10.20.30.40.5……1.0精确解10.9090.8330.7690.6660.625……0.5数值解10.90.8190.7520.6590.6470.4632.龙格库塔法※

基本思想：取Taylor级数展开式前三项近似求解，并利用线性组合代替导数的求解。既可避免计算高阶导数，又可提高数值积分的精度，这就是Runge-Kutta法的基本思想。2.龙格库塔法※r为精度阶次，ai为待定系数，由精度确定；ki用下式表示线性组合等各阶导数不易计算，用下式中ki的线性组合代替1)当r=1时：与Taloy展开式相比较，可得a1=1,则上式成为欧拉递推公式2)当r=2时：将在点展成Taylor级数与台劳公式的二阶展开近似公式相比，可得以下关系：三个方程，四个未知数，解不唯一各个系数的几种取法——见书上。3)r=4时，四阶龙格库塔公式-最常用：仿真中遇到的大多数工程实际问题，四阶龙格库塔法以能满足精度要求，其截断误差o(h5)与h5同数量级。该法可以自启动。

4)、状态空间四阶龙格-库塔递推式若单输入单输出系统的状态空间表达式为：

在仿真中，对于n阶系统，状态方程可以写成一阶微分方程

根据四阶龙格-库塔公式,有T=tk时刻的xi值T=tk+h时刻的xi值

另

状态方程的四阶龙格-库塔公式如下：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求解解:令y1=x，y2=x’1、建立M-文件vdp.m如下：

functiondy=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’,[010],[1;1]);plot(T,Y(:,1),'-‘,T,Y(:,2))3、结果解

1、建立m-文件rigid.m如下：

functiondy=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、取t0=0，tf=12，输入命令：

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、结果如图图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.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时，误差很大当h>0.05后，曲线发散振荡，数值不稳定，完全失去意义系统矩阵的特征值差异较大一般线性常微分方程组：的系数矩阵A的特征值具有如下特征：则称为“病态”方程。定义：2控制系统仿真中的“病态”问题1病态系统中绝对值最大的特征值对应于系统动态性能解中瞬态分量衰减最快的部分，它反映了系统的动态响应和系统的反应灵敏度。一般与系统中具有最小时间常数Tmin的环节有关，要求计算步长h取得很小。2病态系统中绝对值最小的特征值对应于系统动态性能解中瞬态分量衰减最慢的部分，它决定了整个系统的动态过渡过程时间的长短。一般与系统中具有最小时间常数Tmax的环节有关，要求计算步长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=eAT，为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,’选项’)其中选项同上。