控制系统CAD控制系统数字仿真

1、第4章 控制系统数字仿真数字仿真就是采用数学模型，在数字计算机上借助数值解法所进行的仿真实验。所谓数值解法，就是寻求y(t)在a,b区间内的一系列离散节点t1<t2<<tm<tm+1<上的近似值y1, y2, ，ym, ym+1, ，即求取ym+1 y(tm+1)。相邻节点的间距h=tm+1-tm称为步长，这里假定h为定值，即tm=t0+mh，m=0, 1, 2,。本章主要讲述数字仿真的基本理论和方法。4.1 数值积分法系统的动态特性通常用一阶微分方程组来描述，也即状态空间表达式。一般来说，只有极少数的微分方程能用到初等方法得其解析解（或用解析的方法得到精确解），

2、多数只能用近似数值求解。利用计算机求微分方程主要使用数值积分法，它是系统仿真的一个重要方法。在这里，我们主要研究一阶微分方程的形式，如： （1）， 求y(t) 解： 当t=tm+1，t0=tm时 （2）数值积分法时在已知初值的情况下，对f(t, y)进行近似积分，从而对y(t)进行数值求解的方法。下面介绍几种在数字仿真常用的数字积分法。1.欧拉法欧拉法又称为折线法，是最简单，也是最早的一种数值计算方法。对于式（2），如果积分间隔h=tm+1 - tm取得足够小，使得在tm与tm+1之间的f(t, y)可近似看做常数f(tm, ym)。这样式（2）可化为：即 （3）（3）式即为欧拉公式。欧拉公式

3、的几何解释：对于微分方程（1）的解y(t)看作是一条曲线，在任一步长内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长起点的斜率。基于上述的几何解释，我们从初始点 (t0, y0)出发向前推进(t1, y1)点，(t2, y2)点，t0t1t2tm+1t0y0y1y2ym+1y图中阴影部分即为误差。欧拉法的特点是：计算简单，但精度较低。例：用欧拉法求解初值问题（），h=0.1。解：因为则，该题解为：，将准确解y(tm)与近似解ym一起放入下表，可得：tm ym y(tm)tm ym y(tm) 0.1 1.1000 1.09540.2 1.1918 1.18320.3 1.2

4、774 1.26490.4 1.3582 1.34160.5 1.4351 1.4142 0.6 1.5090 1.48320.7 1.5803 1.59420.8 1.6498 1.61250.9 1.7178 1.67331.0 1.7848 1.7321由此表可以看出欧拉公式的精度很差。2.后退的欧拉法若用f(tm+1,ym+1)来代替f(tm,ym)，则（3）式可变为： （4）则（4）式称为后退的欧拉公式。后退的欧拉公式是隐式的（因为（4）式右边的ym+1是未知的），此时通常需要用迭代法求解，即： （5）t0t1t2tm+1t0y0y1y2ym+1y后退的欧拉公式的几何解释：在任一步长

5、内，用一段直线代替函数y(t)的曲线，此直线段得斜率等于该函数在步长终点的斜率。例：用后退的欧拉法求解初值问题（），h=0.1。解：后退的欧拉公式和欧拉公式的精度相同，都是一阶精度。3.梯形法比较欧拉公式和后退的欧拉公式可知，如果对这两种方法进行算术平均，即可大大消除主要误差，从而获得更大的精度，这种方法通常称为梯形法，其计算公式为： （6）同后退的欧拉公式一样，梯形公式也是隐式的。此时通常采用欧拉公式先预报一个，再将预报的代入（6）式进行校正，求出ym+1。梯形法的迭代公式为： （7）（7）式又被称为预估-校正公式。显然，梯形法要比欧拉法和后退的欧拉法精度更高，但计算量比欧拉法大。例：用梯形

6、法求解初值问题（），h=0.1。解：tm ym y(tm)tm ym y(tm) 0.1 1.0959 1.09540.2 1.1841 1.18320.3 1.2662 1.26490.4 1.3434 1.34160.5 1.4164 1.4142 0.6 1.4860 1.48320.7 1.5525 1.59420.8 1.6153 1.61250.9 1.6782 1.67331.0 1.7379 1.73214.龙格库塔法则：我们把称作tm,tm+1上的平均斜率，由此可见，只要给平均斜率提供一种算法，便可导出一种计算公式。当=0时，平均斜率为f(tm,ym)，为欧拉公式；当=1时，

7、平均斜率为f(tm+1,ym+1)，为后退的欧拉公式；若采用tm和tm+1两点斜率值的算术平均作为平均斜率，则平均斜率为，此时为梯形公式（精度更高）。这给我们启示，如果设法在tm,tm+1区间内多预测几个点的斜率值，然后将它们进行加权平均作为平均斜率，则可构造出精度更高的计算公式，这就是龙格库塔的基本思想。（1）二阶龙格库塔法随意考查区间tm,tm+1区间内一点：，用tm和tm+p两个点的斜率值K1，K2线性组合得到平均斜率，则 （8）因为：所以：（8）又因为： （9）（根据泰勒公式）比较（8）式跟（9）式，可得： （10）满足（10）式的一族公式（8）统称为二阶龙格库塔公式。当p=1时，此时

8、，此时：即梯形公式是二阶龙格库塔的一个特例。另外给出一个常用的二阶龙格库塔公式：（）（2）三阶龙格库塔法为了进一步提高精度，取tm,tm+1区间内三个点，即tm，tm+p= tm +ph，tm+q= tm +qh，其中，0<pq1，用这三点的斜率值K1，K2，K3的线性组合得到平均斜率，此时： （11）经过推导可得： （12）满足（12）式的一族公式（11）统称为三阶龙格库塔公式。另外给出一个常用的三阶龙格库塔公式：（3）四阶龙格库塔法同理，取 tm,tm+1区间内四个点的斜率值进行线性组合来得到平均斜率，即为四阶龙格库塔法。因其计算很复杂，直接给出最常用的四阶龙格库塔公式： （13）例

9、：用四阶龙格库塔法求解初值问题（），h=0.2。解：tm ym(1) ym(2)y(tm)0.10.20.30.40.50.60.70.80.91.01.09591.18411.26621.34341.41641.48601.55251.61531.67821.73791.18321.34171.48331.61251.73211.18321.34161.48321.61251.7321四阶龙格库塔法的计算量比梯形法大将近一倍，但由于放大了步长，使得计算量几乎相同。由此可以看出选择算法的重要性。关于数值积分法的几点讨论：（1）单步法和多步法以上介绍的数值积分法都有一个特点，每次计算只需用到前一

10、次的计算结果，而不需要更前面的计算结果，我们称这类方法为单步法。优点：占用内存少能自启动（只要知道初值就可自行进行计算）容易实现变步长计算 多步法比单步法更精确（因为多步法利用的信息量更大）（2）显式和隐式 显式算法可直接利用递推公式求解下一步结果，而隐式算法需要采用迭代法。隐式算法计算过程复杂、计算速度慢，但通常精度较高。（3）数值积分的稳定性例：， 0t1.5其精确解为：取步长h=0.1，则欧拉法：y1.5=-1.09225×10-4四阶龙格库塔法：y1.5=3.95730×101精确解为：y(1.5)=9.54173×10-21可见此时数值积分的结果是不正确

11、的，或是不稳定的。原因：数值积分是一种近似积分，它在每一步的递推中将引入误差，所以若误差的积累越积越大，计算将出现不稳定。现在我们以下面的微分方程为例，来分析上面所述的各种数值积分法的稳定性， <0欧拉法则：所以：所以合理选择步长h，使其满足是保证欧拉法稳定的重要条件。梯形法则：所以：同理即梯形法在任何步长下都是稳定的。四阶龙格库塔法 合理选择步长h，使其满足可保证四阶龙格库塔法稳定。（4）步长与计算精度之间的关系h0h舍入误差截断误差采用数值积分法在计算机上进行仿真，主要包括两种误差：舍入误差与截断误差。这两种误差与步长之间的关系可用右图表示。从图中可以看出，两种误差对步长的要求是相互

12、矛盾的，但两者之和有一个最小值，这样步长最好能选在最小值（当然实际要做到这一点很难）。4.2 控制系统的结构及其描述4.2.1控制系统常见的典型结构形式（1）单输入单输出开环控制结构y(t)u(t)r(t)控制器被控对象（2）单输入单输出前馈控制结构y(t)u(t)r(t)控制器被控对象补偿装置（3）单输入单输出闭环控制结构-y(t)u(t)r(t)控制器被控对象反馈环节凡是单输入单输出控制系统结构均可采用拓扑结构进行描述。拓扑结构不关心事物的细节，只讨论事物之间的相互关系，并通过图表示出来。通俗的讲就是如何把这些典型环节连接起来。-123456（4）多输入多输出控制结构用拓扑结构表示多输入多

13、输出系统，可对系统仿真带来极大方便。4.2.2控制系统的连接矩阵（1）单输入单输出系统-123456根据上图中的ui和yi之间的连接关系，可得：写成向量形式：即 U=WY+WOr, 其中，W称为系统连接矩阵，为n×n型，阵中各元素表示各环节之间的拓扑连接关系；WO为输入连接矩阵，为n×1（此例为一个输入），U=u1,u2,unT为各环节输入向量，Y =y1,y2,ynT为各环节输出向量。（2）多输入多输出系统r3r2r1-1234 35 2即 U=WY+WOr从W阵可以看出各环节之间的连接关系，即wij=0，环节i与环节j不相连；wij0，环节i与环节j相连，且wij>

14、;0表示有正的连接关系，wij<0表示有负的连接关系；wii0，表示各环节i有自反馈。采用连接矩阵表示各环节的连接关系，使得对复杂系统结构的仿真变得简明方便。4.3 面向系统结构图的数字仿真4.3.1典型闭环系统的数字仿真1.典型闭环系统的结构形式典型闭环系统结构如下图所示：y-erG(s)v图中，V为系统的反馈系数，设其为常数。设图中开环传函为： （14）2.系统仿真模型与求解思路由（14）式可得能控标准型为：式中，其中，A、C阵中的，为（14）式首一化后的分子分母多项式系数，即， 。令 u=e则 u=e=r-v y令，则 求解思路：（1）欧拉公式：（2）梯形公式：（3）四阶龙格库塔公

15、式K1，K2，K3，K4为状态变量X在t=tk时的四组斜率，由此可求出tk+1时刻的状态变量Xk+1，从而求得tk+1时刻的输出：yk+1=CXk+1。所以只要给出了状态变量的初值，就可以由上面的递推公式，求出各时刻的状态变量Xk和输出yk。系统结构图开环传函开环状态空间表达式闭环状态空间表达式采用数值积分法（如四阶龙格库塔法）进行仿真3.仿真程序框图与实现（1）程序框图（参见书P133页图4.9）（2）程序 输入数据a=a0,a1,.,an;b=b0,b1,.,bm; X0=x0,x1,.,xn;V=V0;n=n0;T0=t0;Tf=tf;h=h0;r=r0;形成开闭环系数矩阵b=b/a(1

16、);a=a/a(1);A=a(2:n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(1,n-1),1'm1=length(b);C=fliplr(b),zeros(1,n-m1);A1=A-B*V*C;X=X0;y=0;t=T0;运算求解N=round(Tf -T0)/h);for i=1:N K1=A1*X+B*r; K2=A1*(X+h*K1/2)+B*r; K3=A1*(X+h*K2/2)+B*r; K4=A1*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y;C*X; t=t;t(i)+h;

17、end输出结果t,yplot(t,y)例题（书P135）程序如下：程序1：k=1;a=conv(1 0 0,conv(0.25 1,0.25 1);b=2*k, k;X0=0 0 0 0'V=1;n=4;T0=0;Tf=300;h=0.3;r=1;B=b/a(1);a=a/a(1);A=a(2:n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(1,n-1),1'm1=length(b);C=fliplr(b),zeros(1,n-m1)A1=A-B*V*CX=X0;y=0;t=T0;N=round(Tf-T0)/h);for i=

18、1:N K1=A1*X+B*r; K2=A1*(X+h*K1/2)+B*r; K3=A1*(X+h*K2/2)+B*r; K4=A1*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y,C*X; t=t,t(i)+h;endplot(t,y)程序2：r=1;num1=2 1;den1=0.0625 0.5 1 0 0;num2=1;den2=1;num,den=feedback(num1,den1,num2,den2);A,B,C,D=tf2ss(num,den)Tf=100; % input('仿真时间Tf=');h=0.25; % in

19、put('计算步长h=');X=zeros(length(A),1);y=0;t=0;for i=1:Tf/h K1=A*X+B*r; K2=A*(X+h*K1/2)+B*r; K3=A*(X+h*K2/2)+B*r; K4=A*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y;C*X; t=t;t(i)+h;endplot(t,y)4.3.2复杂连接的闭环系统数字仿真基本思路：与实际系统的结构图相对应，只需将各环节的参数及各环节间的连接方式输入计算机，仿真程序就能自动求出闭环系统的状态空间表达式，从而运用数值积分法进行求解。（1）典型环

20、节比例环节： 积分环节： 比例积分环节： 一阶惯性环节： 一阶超前滞后环节：二阶振荡环节： 对于二阶振荡环节， yuyu-yu可见，我们通常遇到的都是一阶环节，即使是二阶振荡环节，也可用两个一阶环节等效连接得到，因此在本章中，我们认为所遇到的环节均为一阶环节。典型的一阶环节（通用）yiui用此格式可描述所有环节的输入输出关系。设输入向量为U=u1,u2,unT输出向量为Y =y1,y2,ynT则系统中所有环节的输入输出关系可表示为：(A+BS)Y=(C+DS)U （15）其中，（2）仿真求解把系统中各环节描述出来后，还需要进行数值积分求解，此时需要把各环节的相互关系表示出来，即用上节所讲的连接

21、绝阵。举例说明：yr-123 3 2即 U=WY+WOr （16）, 连接矩阵及模型参数确定之后，求取闭环系统的状态空间表达式（仿真模型）。将（16）式代入（15）式，得：设Q=B - DW，且逆矩阵存在，H=CW A，则为不增加求解的复杂性，应使DW0=0，则设，则求解一阶微分方程组的矩阵形式，就可采用数值积分法进行求解。应注意两点： 为保证Q-1的存在，应按照bi0的原则确定每个典型环节，即避免以纯比例、纯微分环节作为典型环节。 为保证DWO=0，应使含有微分项系数（即 di0）的环节不直接与外加输入相连。（3）仿真程序框图与实现 系统模型参数矩阵输入处理一般将系统中各典型环节的系数整理成

22、如下矩阵的形式（假设系统由n个典型环节组成）：，A=diag(P(:,1)或 a=a1,a2,anA=diag(a)然后由程序自动转换生成A、B、C、D阵（diag函数）。 连接矩阵输入方法系统连接矩阵W和输入连接矩阵WO中大部分元素为零，非零元素较少，为加快输入速度，通常采用只输入非零元素的方式来构成WIJ阵（由W和WO阵中的非零元素及对应的环节号所构成的矩阵），然后由程序自动转换成W和WO。具体作法是建立非零元素矩阵WIJ阵（m×3）i j wij其中m为W和WO中非零元素的个数，i为被作用环节，j为作用环节（外加输入用0来表示），wij为对应的连接系数i j wij 程序框图（

23、参见书P141页图4.16） 程序实现P=a1,b1,c1,d1; a2,b2,c2,d2; an,bn,cn,dnWIJ=i, i, wij ; ; Y0=y10, y20,., yn0 ;L=L0 ;n=n0 ;T0=t0 ;Tf =tf ;h=h0 ;r=r0 ;nout=nput;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if WIJ(k,2)=0 W0(WIJ(k,1)=WIJ(k,3); else W(W

24、IJ(k,1),WIJ(k,2)=WIJ(k,3); endendQ=B-D*W;Q1=inv(Q);H=C*W-A;A1=Q1*H;B1=Q1*C*W0;Y=Y0;y=Y(nout);t=T0;e=0;N=round(Tf-T0)/h*L);for i=1:Nfor j=1:L K1=A1*Y+B1*r; K2=A1*(Y+h*K1/2)+B1*r; K3=A1*(Y+h*K2/2)+B1*r; K4=A1*(Y+h*K3)+B1*r; Y=Y+h*(K1+2*K2+2*K3+K4)/6;en y=y;Y(nout); t=t;t(i)+h*L; e=e;r-Y(nout);endt, yp

25、lot(t,y,t,e)4.4 连续系统的离散相似法4.4.1连续系统状态方程的离散化该状态方程的解为：当t=kT时， 当t=KT时，设t=-kT，即= t+kT，则d = dt 若采用零阶保持器，则 u(t+kT)= u(kT)，0t<T其中， 若采用一阶保持器，则其中，y(k+1)T=CX(k+1)T+Du(k+1)T所以若采样零阶保持器，则若采样一阶保持器，则4.4.2典型环节的离散化对阶次较高的系统，求取并非易事，再加上它是按连续系统进行离散化的，不适用于非线性系统。所谓按典型环节离散化，就是将系统分成若干个典型环节，然后针对每个典型环节进行离散化。这种方法具有改变参数方便，易于

26、引进非线性环节等特点。设某典型环节的传函为：设 （17）则 （18）由（17）式得： （19）由（18）式可得： （20） 则由（19）和（20）式可得：按离散化步骤，可得：其中，对于典型环节来说，只要知道 A、B、C、D值并代入上式就可得离散化系数阵。（1）积分环节。 （A=0，D=0）（2）比例积分环节。（A=0，D0）（3）惯性环节。（A0，D=0）（4）一阶超前滞后环节。（A0，D0）4.4.3离散相似法的框图与实现（1）程序框图（参见书P149页图4.18）（2）各环节输入作用Uk=WYk+Wor式中，Uk=u1k,u2k,unkT为各环节的输入量，n为输入环节数，Yk=y1k,y2

27、k,ynkT为各环节的输出量，r为外加参考输入。 （3）程序实现for i=1:n （A=P(:,1)，B=P(:,2)，C=P(:,3)，D=P(:,4)） if(A(i)=0) FI(i)=1; FIM(i)=h*C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1; FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); end else FI(i)=exp(-h*A(i)/B(i); FIM(i)=(1-FI(i)*C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;

28、FID(i)=0; if(D(i)=0) FIC(i)=1- （A(i)*D(i) ）/ （C(i)*B(i)）; FID(i)=D(i)/B(i); end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;t=T0:h*L:Tf;N=length(t);for k=1:N-1 for p=1:L Ub=Uk; Uk=W*Y+W0*r; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Y=FIC'.*X+FID'.*Uf; end y=y,Y

29、(nout);endplot(t,y) 4.5 非线性系统的数字仿真4.5.1典型非线性环节及其仿真实现1.饱和非线性0-ssuy数学描述为：其中，s为饱和值。NNYYy=sy=-sy=uu>0?返回调用|u|s?MATLAB程序：function y=saturation(u,s) if(abs(u)>=s) if(u>0) y=s; else y=-s; end else y=u; end0-ssuy2.死区非线性数学描述为： 其中，s为死区值。NNYYy=u-sy=u+sy=0u>0?返回调用|u|s?MATLAB程序：function y=dead(u,s) i

30、f(abs(u)>=s) if(u>0) y=u-s; else y=u+s; end else y=0; end0-ssuy3.滞环非线性数学描述为： 其中，s为滞环宽度；yb表示非线性环节的前一次输出值，其值需要判断u、y的变化率大小。YNNYNYu-syb?Ny=u-sy=yb返回调用u>ub?Yu<ub?u+syb?y=u+syb=yMATLAB程序：function y,yb=backlash(yb,u,ub,s) if(u>ub) if(u-s)>=yb) y=u-s; else y=yb; end else if(u<ub) if(u+s

31、)<=yb) y=u+s; else y=yb; end else y=yb; end end yb=y;0-ssuy4.继电非线性数学描述为： NNYYy=sy=-s返回调用u>0?u<0?MATLAB程序：function y=sign(u,s) if(u>0) y=s; else if(u<0) y=-s; end4.5.2非线性特征的判断本章的程序只考虑饱和、喝死、滞环三种典型非线性特性，并规定非线性环节不单独作为一个仿真环节，而是根据系统结构情况，把非线性环节添加到相应的线性典型环节之前或之后，因此，对于有非线性环节的系统并不增加典型环节数。这样，对于一

32、个典型环节除了要给出参数a、b、c、d外，还要给出其附加的非线性环节标志和参数，以便在程序中识别并执行。在输入数据时，设立非线性标志向量Z和参数向量S：非线性标志向量Z：Z=z1,z2,znZ中的每一个元素对应一个环节的非线性特征，n为环节序号，为此约定如下：zi=0 典型环节前后无非线性环节zi=1 典型环节前有饱和非线性环节，应修正ui zi=2 典型环节前有死区非线性环节，应修正uizi=3 典型环节前有滞环非线性环节，应修正uizi=4 典型环节后有饱和非线性环节，应修正yizi=5 典型环节后有死区非线性环节，应修正yizi=6 典型环节后有滞环非线性环节，应修正yi参数向量S：S=

33、s1,s2,snS用来存储非线性环节的参数（1）求得ui后，由zi判断各环节输入端是否有非线性环节，若有，则根据标志位确定类型，转相应子程序，修正ui（2）求得yi后，由zi判断各环节输出端是否有非线性环节，若有，则根据标志位确定类型，转相应子程序，修正yi（3）各种非线性环节的实线程序以函数文件的形式存储起来，供主程序调用。4.5.3程序框图及仿真实现1.程序框图（参见书P156页图4.27）2.程序实现Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=zeros(n,1);%设立一个向量单元Ubb存储非线性环节的前一次输出值t=T0:h*L:Tf;N=lengt

34、h(t);for k=1:N-1 for p=1:L Ub=Uk; %设立一个向量单元Ub存储非线性环节的前一次输入值 Uk=W*Y+W0*r; for i=1:n if(Z(i)=0) if(Z(i)=1) Uk(i)=saturation(Uk(i),S(i); end if(Z(i)=2) Uk(i)=dead(Uk(i),S(i); end if(Z(i)=3) Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i); end end end Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+

35、FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; for i=1:n if(Z(i)=0) if(Z(i)=4) Y(i)=saturation(Y(i),S(i); end if(Z(i)=5) Y(i)=dead(Y(i),S(i); end if(Z(i)=6) Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i); end end end end y=y,Y(nout);endplot(t,y)4.6 计算机控制系统的数字仿真4.6.1 采样控制系统的仿真方法1.差分方程递推求解法已知，即，则所以，则所

36、以又因为，所以若分别得到控制器的传递函数D(z)和被控对象的传递函数G(z)时，此时可通过求解下面一组高阶差分方程进行仿真求解。写成向量形式：（1）程序框图（参见书P162图4.33）（2）程序实现clear allR=r0;T=0.1;T0=t0;Tf=tf;A=a1,a2,an;B= b0,b1,bm;C= c1,c2,cl;D= d0,d1,dr;E=zeros(r+1 ,1); % length(D);U1=zeros(L,1); % length(C);U2=zeros(m+1,1); % length(B);Y=zeros(n,1); % length(A);yk=0;y=0;t=

37、0;ek=0;uk=0;N=length(T0:Tf/T);for k=1:N ukb= uk;ek=R-yk; E=ek;E(1:r); uk=-C*U1+D*E; U1=ukb;U1(1:L-1); U2=uk;U2(1:m); yk=-A*Y+B*U2; Y=yk;Y(1:n-1); y=y,yk; t=t,k*T;endplot(t,y)采用差分方程递推求解的结果无截断误差。从理论说算法时精确的，只要D(z)和G(z)已知即可，缺点是当系统复杂时，G(z)难以求取，且不便考虑非线性环节。2.连续部分按环节离散化当系统连续部分较复杂时，不用化简求取G(z)，可按照连续系统典型环节离散化仿真方法，将连续部分中各环节离散化处理后，与采样部分一并进行仿真。（1）程序框图（参见书P163图4.35）（2）程序实现clear allr=2;P=0 1 1 0;1 1 1 0;G=2.72,-1;F=0.717;W=0 0;1 0;W0=1 0'h=0.01;T0=0;Tf=30;Tm=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);n=length(A);m1=length(G);m2=length(F);f