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

第四章 连续系统的离散化方法 4.1 常微分方程的数值解法 00 ( , ) () dx f t x dt x t x 一 . 数值求解的基本概念 设微分方程为 则求解方程中函数 x(t)问题的常微分方程初值问题 所谓数值求解就是要在时间区间 a, b中取若干离散点 求出微分方程在这些时刻的近似值 0 1 2 Nx x x x 01 Na t t t b 2 ( 2 ) ( ) 0 0 0 0 0( ) ( ) ( ) ( ) ( )2 ! ! k khhx t h x t h x t x t x t k nnn tth 1 取前两项近似： 1 ( , )k k k kx h f t x x 这种方法的几何意义就是把 f(t,x)在 区间 tk,tk+1内的曲边面积用矩形面 积近似代替。计算简单，计算量小， 而且可以自启动。当 h很小时，造成 的误差是允许的。该算法具有一阶 精度。 取 k=0,1,2,N, 从 t0开始，逐点递推求解 t1时的 y1， t2时的 y2, 直至 tn时的 yn，称之为欧拉递推公式。 矩形面积 1. 欧拉法 欧拉法的特点： 导出简单，几何意 义明显，便于理解，能说明构造数 值解法一般计算公式的基本思想 。 通常用它来说明有关的基本概念。 例 设系统方程为 用 Euler法求其数值解（取步长 , ) 1.0h 2 , 0 1x x x 10 t 递推公式为 1 , 1 0 . 1n n n n n nx x h f t x x x 则 4628.01.01,0.1 7519.0819.01.01819.01.01,3.0 819.091.09.01.01,2.0 9.01.01,1.0 1,0 991010 2233 1122 0011 00 yyyt yyyt yyyt yyyt yt 已知方程的解析解为 精确解和解析解作比较： 误差在 数量级， 精度较差。 210 ty 1 1 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. 龙格库塔法 kktxff r为精度阶次， ai为待定系数，由精度确定； ki用下 式表示 线性组合 1 12 1 ( , ) , 2 , 3 , i i k k j j k f t b h x h b k i r 2 1 ()2! kkk k k t x k hx x h f f f f 1 1 r k k i i i x x h a k 2 ( 2 ) 0 0 0 0( ) ( ) ( ) ( )2 hx t h x t h x t x t 等各阶导数不易计算，用下式中 ki的线性组合代替 1)当 r=1时： 1 1 1 ,kkx x h a k 1 ( , ) ,kkk f t x 与 Taloy展开式相比较，可得 a1=1,则上式成为 11 ( , ) ,k k k k kx x h k x h f t x 欧拉递推公式 2)当 r=2时： 1 2 1 2 1 kk kk k f t , x k f t b h , x h b k 2 1 2 1kkk t xk f b h f b h k f 1 2 1kkf t b h , x h b k kkt ,x将 在点 展成 Taylor级数 与台劳公式的二阶展开近似公式相比，可得以下关系： 12 22 21 1 12 12 aa ba ab 三个方程，四个未知数，解不唯一 各个系数的几种取法 见书上。 3) r=4时，四阶龙格库塔公式 -最常用： 仿真中遇到的大多数工程实际问题，四阶龙格库塔法以能 满足精度要求，其截断误差 o(h5) 与 h5同数量级。该法可以 自启动。 1 1 2 3 4 1 21 32 43 22 6 22 22 kk kk kk kk kk h x x ( K K K K ) K f t , x hh K f t , x K hh K f t , x K K f t h , x hK 4)、状态空间四阶龙格 -库塔递推式 若单输入单输出系统的状态空间表达式为： X A X B U Y C X D U 在仿真中，对于 n阶系统，状态方程可以写成一阶微分方程 1 1 2 2 1 2 3( , , , , , ) ( 1 , 2 , , ) i i i in n i in x a x a x a x b u f t x x x x i n 根据四阶龙格 -库塔公式 ,有 1 1 2 3 4 1 1 1 2 2 2 1 1 1 1 3 1 1 2 2 4 1 1 3 3 ( 2 2 ) ( 1 , 2 , , ) 6 () ( ) ( ) ( ) 2 2 2 ( ) ( ) ( ) 2 2 2 ( ) ( kk i i i i i i k k k i i i in n i k kk i i i in n i i k kk i i i in n i i k kk i i i in n i h x x k k k k i n k a x a x a x b u t h h h k a x k a x k b u t h h h k a x k a x k b u t k a x hk a x hk ) ( ) ik b u t h T=tk时刻的 xi值 T=tk+h时刻的 xi值 另 1 1 1 1 2 1 2 1 1 2 1 1 2 1 2 2 2 2 3 1 3 2 3 3 4 1 4 2 4 4 , , , k k k T k k k T nn TT nn x x x x x x k k k k k k k k k k k k k + 1 k 1 xx KK 状态方程的四阶龙格 -库塔公式如下： 2 3 4 2 32 43 ( 2 2 ) 6 () ( ) ( ) 22 ( ) ( ) 22 ( ) ( ) k k k k h ut hh ut hh ut h u t h k + 1 k 1 1k k1 k k k + 1 k + 1 x x K K K K K A x B K A x K B K A x K B K A x K B y C x 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求解 例 1 经典的微分方程 2 2 2 2 1 0 01 01 d x dx ( x ) x dt dt x ( ) ; x ( ) 解 : 令 y1=x， y2=x 则微分方程变为一阶微分方程组： 12 2 2 1 2 1 12 21 0 2 0 0 y y y ( y ) y y y ( ) , y ( ) 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)0(,1)0(,0)0( 51.0 321 213 312 321 yyy yyy yyy yyy 解 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、 取 t0=0， tf=12，输入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+) 3、结果如图 图中， y1的图形为实线， y2的图形为“ *”线， y3的图形为“ +”线 . 4.2 数值算法的稳定性及求解原则 1.数值算法的稳定性 特征根在 s平面的左半平面，系统稳定。 (1)欧拉法： 稳定： (2)梯形法： 恒稳 11zh ( ) ( 1 ) ( )zX z h X z 1z ( ) ( 1 ) ( 1 ) ( ) 22 hhz X z X z 2.数值算法的选择原则 Matlab提供了微分方程数值求解的一般方法 ,作为仿真算 法的使用者 ,可不必考虑算法具体实现 ,而应关心各种方法在 使用中会出现的问题 ,以及如何在仿真中恰当的选用这些方法 . 一般 ,选用数值算法从以下几个方面考虑 : (1)精度 受算法和 h影响 截断误差 +舍入误差 =累计误差 (2)计算速度 受算法和 h影响 算法简单，速度就快些。 (3)稳定性 受 h影响，一般 h （ 2-3） 系统最小时间 4.3 数值算法中的 “ 病态 ” 问题 1 “病态”常微分方程 例： ( 0) (1 , 0 , 1 ) T X AX X 其中 1 2 3( , , )TTX x x x - 2 1 1 9 2 0 A = 1 9 2 1 2 0 4 0 4 0 4 0 采用四阶龙格库塔法 h=0.01时，计算时间长 h=0.04时，误差很大 当 h0.05后，曲线发散振荡，数值不稳定，完全失去意义 系统矩阵的特征值差异较大 一般线性常微分方程组： 00( ) ( ) ( ) , ( )X t A X t B U t X t X 的系数矩阵 A的特征值具有如下特征： 11 Re ( ) 0 m a x Re ( ) m in Re ( ) i iiinin 则称为“病态”方程。 定义： 2 控制系统仿真中的“病态”问题 1 病态系统中绝对值最大的特征值对应于系统动态性能 解中瞬态分量衰减最快的部分，它反映了系统的动态响 应和系统的反应灵敏度。一般与系统中具有最小时间常 数 Tmin的环节有关，要求计算步长 h取得很小。 2 病态系统中绝对值最小的特征值对应于系统动态性能 解中瞬态分量衰减最慢的部分，它决定了整个系统的动 态过渡过程时间的长短。一般与系统中具有最小时间常 数 Tmax的环节有关，要求计算步长 h取得很大。 3 对于病态问题的仿真需要寻求更加合理的算法，以解 决病态系统带来的选取计算步长与计算精度，计算时间 之间的矛盾。 3 “病态”系统的仿真方法 采用稳定性好，计算精度高的数值算法，并且允许计 算步长能根据系统性能动态变化的情况在一定范围内作相 应的变化，采用 隐式吉尔法 1 0 1 1 1 1k k k r k r kx a x a x a x h f 该法已经证明对病态方程求解过程是数值稳定的。 隐式吉尔法从理论上十分适应于病态系统 ，但需要解决好 以下问题 (1) 自启动 r阶多步算式无法自启动，需要用单步法求出前 r步值 (2) 预估迭代 迭代方法要求收敛性良好，否则在大步长时会造成数 值发散。 (3) 变步长 初始阶段采用小步长，随后可逐步放大步长。 对不同精度要求的系统仿真，要考虑变阶次问题， 即为减小每一步计算的截断误差，以提高精度，应选用 较高的阶次，而当精度较低时，为减少工作量，则应选 取较低的阶次。仿真时应根据估计误差 与给定的误差精 度相比较改变步长或阶次来重新计算。 4.4 连续系统状态方程的离散化 上章所述的连续系统数学模型的离散化， 是通过数值积分法实现的，尽管面向结 构图的仿真方法是按环节给定参数，但 是在计算时还是 按整个系统进行离散化 ， 这就不便于引进非线性环节以进行非线 性系统的仿真。在本节，将介绍连续系 统离散模型的建立和仿真。 开 始 输入开环传函分母、分子系数 ji ca , 输入反馈函数 V 、状态初始值 X 0 求开环状态方程系 数阵 A , B , C 求闭环状态方程系数阵 B v CAA b 输入初始时间 0 t 、终止时间 f t 、 计算步长 h 、输入幅值 r 求龙格 库塔法各次斜率 4,3,2,1, jK j 求 11 , kk yX 输出数据、曲线 1 kk XX f t 到否 ？ N Y 结 束 数值积分法叠代求解 h改变时，叠代过程重复求解，费时繁琐 不能对非线性环节单独考虑。 连续系统离散化 思想：用差分方程描述连续系统的状态方程模型 （因为差分方程的主要特点就是方程中各变量由各相邻时刻的变 化量制约，这相当于递推方程） 1、连续系统的离散化 DuCXy BuAXX 设连续系统状态方程为 00 )( XtX 其中 为状态初始值 . 则由现代控制理论基础知，状态变量 X(t)的解为 () 0 0 ( ) ( ) kT A k T A k Tk T e e d X X B u t tAtA dee 0 )( 0 )()( uBXtX 而 t = ( k + 1 )T 时，可表示为 ( ) 0 ( ( 1 ) ) ( ) ( ) T TTk T e k T e d k T AAX X B u 当系统输入 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 时的状态转移矩阵 () 0 T TH e d A B 离散后的状态空间表达式为： ( 1 ) ( ) ( ) ( ) ( ) ( ) x k T G x k T H u k T y k T C x k T D u k T 、 Matlab表示 DuCXy BuAXX 已知连续系统状态方程为 在采样周期 T下离散后的状态空间表达可表示为： ( 1 ) ( ) ( ) ( ) ( ) ( ) x k T G x k T H u k T y k T C x k T D u k T () 0 G, T TTe H e d AA B其 中 ， 在 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时进行离散化。 63 1 2 5 ( s )G ( s ) ( s ) ( s ) ( s ) 解：利用以下 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);