民用爆炸物品库房建设项目可行性研究报告

1、第4章 控制系统数字仿真数字仿真就是采用数学模型，在数字计算机上借助数值解法所进行的仿真实验。所谓数值解法，就是寻求y(t)在a,b区间内的一系列离散节点t1t2tmtm+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）式即为欧拉公式。欧拉公式的几何解释：对于微分方程（1）的解y(t)

3、看作是一条曲线，在任一步长内，用一段直线代替函数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.2774 1.26490.4 1.3582

4、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后退的欧拉公式的几何解释：在任一步长内，用一段直线代替函数y(t)的曲线，此直

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

6、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时，平均斜率为f(tm+1,ym+1)，为后退

7、的欧拉公式；若采用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，其中，0pq1，用这三点的斜率值K1，K2，K3的线性组合得到平均斜率，此时： （11）经过推导可得： （12）满足（12）式的一族公式（11）统称为三阶龙格库塔公式。另外给出一个常用的三阶龙格库塔公式：（3）四阶龙格库塔法同理，取 tm,tm+1区间内四个点的斜率值进行线性组合来得到平均斜率，即为四阶龙格库塔法。因其计算很复杂，直接给出最常用的四阶龙格库塔公式： （13）例：用四阶龙格库塔法求解初值问题（），h=0.2。解

9、：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.0922510-4四阶龙格库塔法：y1.5=3.95730101精确解为：y(1.5)=9.5417310-21可见此时数值积分的结果是不正确的，或是不稳定的。原因：数值积分是一种近似积分，它在每一步的递推中将引入误差，所以若误

11、差的积累越积越大，计算将出现不稳定。现在我们以下面的微分方程为例，来分析上面所述的各种数值积分法的稳定性， 0表示有正的连接关系，wij0表示有负的连接关系；wii0，表示各环节i有自反馈。采用连接矩阵表示各环节的连接关系，使得对复杂系统结构的仿真变得简明方便。4.3 面向系统结构图的数字仿真4.3.1典型闭环系统的数字仿真1.典型闭环系统的结构形式典型闭环系统结构如下图所示：y-erG(s)v图中，V为系统的反馈系数，设其为常数。设图中开环传函为： （14）2.系统仿真模型与求解思路由（14）式可得能控标准型为：式中，其中，A、C阵中的，为（14）式首一化后的分子分母多项式系数，即， 。令

12、u=e则 u=e=r-v y令，则 求解思路：（1）欧拉公式：（2）梯形公式：（3）四阶龙格库塔公式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

13、;V=V0;n=n0;T0=t0;Tf=tf;h=h0;r=r0;形成开闭环系数矩阵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*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=

14、X+h*(K1+2*K2+2*K3+K4)/6; y=y;C*X; t=t;t(i)+h;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

15、=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;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; %

16、input(仿真时间Tf=);h=0.25; % input(计算步长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复杂连接的闭环系统数字仿真基本思路：与实际系统的结构图相对应，只需将各环节的参数及各环节间的连接方式输入计算机，仿真程序就能自动求出闭环系统的状态空间表达式，从而运

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

18、互关系表示出来，即用上节所讲的连接绝阵。举例说明：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）仿真程序框图与实现 系统模型参数矩阵输入处理

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

20、j wij 程序框图（参见书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

21、); else W(WIJ(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(nou

22、t);endt, yplot(t,y,t,e)4.4 连续系统的离散相似法4.4.1连续系统状态方程的离散化该状态方程的解为：当t=kT时， 当t=(k+1)T时，设t=-kT，即= t+kT，则d = dt 若采用零阶保持器，则 u(t+kT)= u(kT)，0t0?返回调用|u|s?MATLAB程序：function y=saturation(u,s) if(abs(u)=s) if(u0) y=s; else y=-s; end else y=u; end0-ssuy2.死区非线性数学描述为： 其中，s为死区值。NNYYy=u-sy=u+sy=0u0?返回调用|u|s?MATLAB程序：function y=dead(u,s) if(abs(u)=s) if(u0) y=u-s; else y=u+s; end else y=0; end0-ssuy3.滞环非线性数学描述为： 其中，s为滞环宽度；yb表示非线性环节的前一次