计算机仿真(编写)

3.1系统仿真概述(1)系统的概念系统有不同的定义，按照德谟克利特的理解，系统可以定义为：任何事物都是在联系中显现出来的，都是在系统中存在的，系统联系规定每一事物，而每一联系又能反映系统的联系的总貌。GGolden也给出了一个实用的定义：按照某些规律结合起来，互相作用、互相依存的所有实体的集合或总和”。比如，一个理发馆系统组成为：实体：服务员、顾客顾客：按某种规律到达，服务完毕后顾客离去服务员：根据顾客的要求，按一定的程序服务相互作用：顾客到达模式影响着服务员的工作忙闲状态，顾客排队状态服务员的多少和服务效率：影响着顾客接受服务的质量一个电动机调速系统：实体：电动机、测速元件、比较元件以及控制器。相互作用：实现按给定要求调节电动机的速度归纳以上定义，可以认为系统可由三部分组成：确定边界、输入、输出，而描述系统的“三要素”为：实体、属性、活动一一实体确定了系统的构成，也就确定了系统的边界；一一属性也称为描述变量，描述每一实体的特征：一一活动定义了系统内部实体之间的相互作用，从而确定了系统内部发生变化的过程。(2)模型的概念模型是实际系统本质的抽象与简化，模型主要适用于以下情况：(1)真实的系统尚未建立（2）可能会引起系统破坏或发生故障（3）难以保证每次试验的条件相同（4）试验时间太长或费用昂贵模型分为两大类一一物理模型，采用一定比例尺按照真实系统的“样子”制作一一沙盘模型一一数学模型，用数学表达式形式来描述系统的内在规律。定义如下集合结构：s=（T,x,a,Q,y,b,入） （3-1）T：时间基，描述系统变化的时间坐标T为整数则称为离散时间系统，T为实数则称为连续时间系统X：输入集，代表外部环境对系统的作用。X被定义为R",其中〃G/+,即X代表〃个实值的输入变量。输入段集，描述某个时间间隔内输入模式，是（XT）的子集。Q：内部状态集，是系统内部结构建模的核心。8：状态转移函数，定义系统内部状态是如何变化的。它是映射：3:QxA—Q其含义为：若系统在，。时刻处于状态q,并施加一个输入段刃：＜X,则b（q,仞）表示系统处于11状态。A：输出函数，它是映射：A:QxXxTY输出函数给出了一个输出段集。Y：输出段集，系统通过它作用于环境。系统模型水平可以分为三类：（1）行为水平一一亦称为输入/输出水平将系统视为一个“黑盒”，在输入信号的作用下，只对系统的输出进行测量；分解结构水平将系统看成若干个黑盒连接起来，定义每个黑盒的输入与输出，以及它们相互之间的连接关系;

状态结构水平不仅定义了系统的输入与输出，而且还定义了系统内部的状态集及状态转移函数。按照Or前的思路，模型可分类为卜表形式：表3.1模型分类模型描述变量的轨迹模型形式变量范围模型的时间集合连续离散空间连续变化模型偏微分方程连续时间模型空间不连续变化模型常微分方程差分方程✓离散时间模型离散(变化)模型有限状态机Q马尔可夫链9活动扫描V连续时间模型事件调度进程交互✓(3)仿真的概念1961年，G.W.Morgenthater,首次将仿真技术性地定义为“仿真意指在实际系统尚不存在的情况下对于系统或活动本质的实现”；1978年，K6rn定义“连续系统仿真”为“用能代表所研究的系统的模型作实验”；1982年，Spriet进一步将仿真的内涵加以扩充：所有支持模型建立与模型分析的活动即为仿真活动；1984年，Oren给出了仿真的基本概念框架“建模一实验一分析”，如图3-2所示，并定义“仿真是一种基于模型的活动”。其中，“系统、模型、仿真”三者之间的关系为：系统是研究的对象，模型是系统的抽象，仿真是对模型的实验。传统上意义上：“系统建模”属于系统辨识技术范畴；“仿真建模”即针对不同形式的系统模型研究其求解算法：“仿真实验”包括检验(Verification)：“仿真程序”的检验和致效(Validation)：将仿真结果与实际系统的行为进行比较现代仿真技术：将仿真活动扩展到上述三个方面，并将其统一到同一环境中。对于系统建模，它是：基本定律及系统辨识等方法的计算机程序化；用仿真方法确定实际系统的模型；基于模型库的结构化建模；采用面向对象建模(Object-OrientedModeling)方法，在类库的基础上实现模型拼合与重用。对于仿真建模，它包括：许多新算法和新软件；模型与实验分离技术，即模型的数据驱动(datadriven),而仿真问题又分为两部分：模型与实验：模型又分为两部分：参数模型和参数值而对于仿真实验，将实验框架与仿真运行控制区分开来；实验框架定义一组条件；输出函数的定义也与仿真模型分离开来。同时Oren描述了仿真概念的框架如图3-3所示。其中：“仿真问题描述”即为“仿真建模”；“行为产生”表达了“仿真实验”；“模型行为及其处理”就是“输出处理”fj为，“1模Mh■为及其处理模型行为《仿真数据)物选行为| 姑构行为行为处理：分析、:，小(4)仿真技术的应用仿真技术在系统设计中的应用情况有两种：-新系统设计：提供了强有力的工具。在可行性论证阶段，进行定量比较，为系统设计打下坚实的基础；在系统设计阶段，进行模型实验、模型简化并进行优化设计。-系统改造设计：涉及新的设备、部件或控制装置。利用仿真技术进行分系统实验，即一部分采用实际部件，另一部分采用模型，避免由于新的子系统的投入可能造成对原系统的破坏或影响，大大缩短开工周期，提高系统投入的一次成功率仿真技术在系统分析中的应用，主要体现在：在真实系统上进行试验：在真实系统上试验会破坏系统的正常运行：难以按预期的要求改变参数，或者得不到所需要的试验条件；很难保证每次的操作条件相同，难以对试验结果做出正确的判断：无法复原；试验时间太长、费用太大或者有危险等仿真在教育与训练中的应用，可制作成训练仿真系统：利用计算机并通过运动设备、操纵设备、显示设备、仪器仪表等复现所模拟的对象行为，并产生与之适应的环境，从而成为训练操纵、控制或管理这类对象的人员的系统。三大类：载体操纵型这是与运载工具有关的仿真系统，航空、航天、航海、地面运载工具，以训练驾驶员的操纵技术为主要目的。过程控制型用于训练各种工厂的运行操作人员如电厂、化工厂、核电站、电力网等搏奕决策型企业管理人员(厂长、经理)，交通管制人员(火车调度、航空管制、港口管制、城市交通指挥等)，军事指挥人员(空战、海战、电子战等)。仿真在产品开发及制造过程中的应用，主要体现在虚拟设计及制造(VirtualManufacturing)：是实际制造在计算机上的本质实现，是仿真技术以制造过程为对象的全方位的应用。(5)系统仿真的类型(i)根据模型的物理属性分类物理仿真：按照真实系统的物理性质构造系统的物理模型，并在物理模型上进行实验的过程称为物理仿真。物理仿真的优点是：直观、形象，也称为“模拟"。物理仿真的缺点是：模型改变困难,实验限制多，投资较大。数学仿真：对实际系统进行抽象，并将其特性用数学关系加以描述而得到系统的数学模型，对数学模型进行实验的过程称为数学仿真。计算机技术的发展为数学仿真创造了环境，亦称为计算机仿真数学仿真优点是：方便、灵活、经济。数学仿真缺点是：受限于系统建模技术，即系统数学模型不易建立。半实物仿真：即将数学模型与物理模型甚至实物联合起来进行实验。对系统中比较简单的部分或对其规律比较清楚的部分建立数学模型，并在计算机上加以实现对比较复杂的部分或对规律尚不卜分清楚的系统，其数学模型的建立比较困难，则采用物理模型或实物仿真时将两者连接起来完成整个系统的实验(ii)根据仿真计算机类型分类模拟计算机仿真：模拟计算机本质上是一种通用的电气装置，这是50—60年代普遍采用仿真设备。将系统数学模型在模拟机上加以实现并进行实验称为模拟机仿真。模拟机仿真是一种并行仿真，仿真时，代表模型的各部件是并发执行的。

数字计算机仿真：将系统数学模型用计算机程序加以实现，通过运行程序来得到数学模型的解，从而达到系统仿真的目的。早期的数字计算机仿真则是一种串行仿真，因为计算机只有一个中央处理器(CPU),计算机指令只能逐条执行。数字模拟混合仿真：为了发挥模拟计算机并行计算和数字计算机强大的存贮记忆及控制功能，以实现大型复杂系统的高速仿真，将系统模型分为两部分，其中一部分放在模拟计算机上运行，另一部分放在数字计算机上运行，两个计算机之间利用模/数和数/模转换装置交换信息。(iii)根据仿真时钟与实际时钟的比例关系分类实际动态系统的时间基称为实际时钟，系统仿真时模型所采用的时钟称为仿真时钟，以此将仿真分为：实时仿真：即仿真时钟与实际时钟完全一致，模型仿真的速度与实际系统运行的速度相同当被仿真的系统中存在物理模型或实物时，必须进行实时仿真.亚实时仿真：即仿真时钟慢于实际时钟，模型仿真的速度慢于实际系统运行的速度，也称为离线仿真。超实时仿真：即仿真时钟快于实际时钟，模型仿真的速度快于实际系统运行的速度。(iiii)根据系统模型的特性分类连续系统仿真连续系统是指系统状态随时间连续变化的系统，分为：-集中参数系统模型，一般用常微分方程(组)描述。-分布参数系统模型，一般用偏微分方程(组)描述。其中，离散时间变化模型中的差分模型归为连续系统仿真范畴。•离散事件系统仿真离散事件系统是指在某些随机时间点上系统状态发生离散变化的系统。与连续系统的主要区别在于：-状态变化发生在随机时间点上这种引起状态变化的行为称为“事件”，因而这类系统是由事件驱动的；“事件”往往发生在随机时间点上，亦称为随机事件，因般都具有随机特性；一系统的状态变量往往是离散变化的；系统的动态特性很难用人们所熟悉的数学方程形式描述研究与分析的主要目标是系统行为的统计性能而不是行点轨迹。(6)系统仿真的一般步骤一般情况下，系统仿真可分为以下步骤，如图3-3。建模与形式化：确定模型的边界，模型进行形式化处理；仿真建模：选择合适的算法，确定算法的稳定性、计算精度、计算速度：程序设计：将仿真模型用计算机能执行的程序来描述，程序中要包括仿真实验的要求仿真运行参数、控制参数、输出要求；模型校验：程序调试，检验所选仿真算法的合理，检验模型计算的正确性(Verification)仿真运行：对模型进行实验仿真结果分析：对系统性能作出评价，对模型可信性进行检验(Validation),只有可信的模型才能作为仿真的基础2连续系统仿真2.1连续系统模型描述连续系统是指系统状态变化在时间上是连续的，可以用方程式(常微分方程、偏微分方程、差分方程)描述系统模型。按照式(3T),连续系统数学模型可定义如下集合结构：S=(T,X,Q,0,y,6,九)典型的形式有：.常微分方程一一输入/输出水平dny dn-'y dy cTh d^udtn df-x dt dfi dr' (其中〃为系统的阶次，q(i=0,L2,…为系统的结构参数，与(7=QL2…，力为输入函数的结构参数，它们均为实常数。.传递函数一一输入/输出水平若系统的初始条件为零，对(3-2)式两边取拉氏变换后稍加整理：VC-SjG(s)= =- U(5) < J (33)Za—s」六。(2)式称为系统的传递函数。.状态空间描述——状态结构水平系统内部模型即为状态空间模型。状态空间描述的一般形式为：

状态方程X=状态方程X=AX+BU(3-4)输出方程Y=CX (3-5)式3-4为状态方程，式3-5为输出方程。其中A是n*n维矩阵，B是n*l维矩阵，C是l*n维矩阵。对形如式3-2的单数入单输出的n阶系统，引进n个内部状态变量xl,x2,…,xn,易于将其转换为上述形式的状态方程。作用函数U为单输入u,输出变量Y为单输出。2.2模型结构变换连续系统仿真要将这个系统的模型在计算机上实现出来，首先耍把系统的各种描述形式转换成内部模型状态空间模型，将其称为模型结构变换。.输入/输出水平模型到内部模型的变换一连续系统，它的数学模型如(3-6)式所示dydtn+adydtn+a\dn-xydtn-'今引进〃个状态变量:则有dnydn则有dnydn-xydn~2y= =-a, -a2 -dtn1dtn-x2dl■--any+u(t)=~axxn-出工〃一i anx\+〃«)将上述〃个一阶微分方程写成矩阵形式可得 0 0 1 0 •0 __ 尤0 0 1 •0x}0X2x20x=—・・・・・・・・••・・・*+u*,0 0 0 ••1_xn_Xn1~an-a„-\一q-2…~a\_(3-7)(3-8)外部模型变换到内部模型不唯一，所以仿真模型也不唯一。一个系统有多种实现，最小实现的充要条件是(4B、。为完全能控且完全能观测。

2系统状态初始值变换如果系统是非零初始条件，那么从外部模型变换到内部内部模型还必须考虑如何将给定的初始条件转变为相应的状态变量的初始值。若系统是由如下一般形式的〃阶微分方程来描述：dnyd'^ydydnud'^uduTOC\o"1-5"\h\z0n——+ax——H \-a.—+a„y=cn i-c, rH \-c. l-c„w"df1 zdt°df1 dr't dt"系统初始条件为：=/),〃"(4))=〃，，(i=L2…,〃T)。伴随方程法：一阶微分方程组的状态变量记为为(i=1,2,•••,〃)，如果它们满足如下关系：%]=aoy-cQu (8)Xj=xj+]-ajy+CjU (9)x„=-an+cnu (10)»=5(石+，〃) (11)该状态方程与原方程等价。伴随方程法显式地表示了状态变量与原输入/输出变量及其高阶导数之间的关系，因而易于进行初始值的转换。这样得到状态方程及输出方程：JA一q 1 0 …0—a2 0状态方程及输出方程：JA一q 1 0 …0—a2 0 1 …0—00…1-an 0 0 …0C=—[1,0…0]即其中=CX+Du其中6,2—。0。25=D=。/。0设为=1,初值转换方程:伴随方程有多种形式，因而得到的状态方程也不唯一。那么，实现这种初值转换的条件是什么呢？考虑转换后得到的系统状态空间模型为:X=AX+BUY=CX即假定〃的〃阶导数项的系数4=0,已知系统的初始条件为:y(o),夕(o),…严，(0)〃(0),«(0),…〃"(0)则为了由上述初始值求出状态变量Xi,x2---xn的初始值，可列出以下方程：y(/)=Cx(t)y(t)=Cx(t)=CAx(t)+CBu(t)y(z)=Cx(t)=CAx(t)+CBu(t)=CA2x(t)+CABu(t)+CBu(t)于是可得下列矩阵方程y(n=6>X(r)+Tu(r) (⑶其中y(O=[y(O即)…u(r)=[«(/)«(/)…"""(nF

0=0=CCA为n〃方阵CAn1-00…oCB0…oT=CABCBCAn-2BCAn3BCB0由(16)式可得： x(/)=^1[y(r)-Tu(/)] (14)即，若 夕।存在，则可由(14)式求出x(t)的初始值。由控制理论可知，9是(A、B、。的能观判别阵，若(A、B、C)是完全能观的，则0非奇异。这就是说，由高阶微分方程输入/输出变量初始值转变为状态初始值的条件是：内部模型(A、B、C)是完全能观的。3.3经典的连续系统仿真建模方法学其主要内容包括：连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。3.3.1离散化原理及要求问题：数字计算机在数值及时间上的离散性用于处理被仿真系统数值及时间上的连续性。连续系统仿真，从本质上，对原连续系统从时间、数值两个方面对原系统进行离散化并选择合适的数值计算方法来近似积分运算，使离散模型近似于原连续模型。相似原理:设系统模型为：7=f(y,u,t),其中〃为输入变量，为系统变量；令仿真间间隔为力，离散化后的输入变量为6(乙)，系统变量为和.)，其中，”表示£=〃瓦如果 虱，”四“),y(tn)^y(tn)即，(乙)=力(乙)一〃(乙)=0,ey(tn)=y(tn)-y(tn)^O(对所有〃…)则可认为两模型等价。原连续模型炉=f(y,»,/),©Ai仿真模型j=/(女日,/“)虱t")图三-1相似原理

对仿真建模方法三个基本要求：(1)稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。(2)准确性：有不同的准确性评价准则，最基本的准则是：绝对误差准则：Cy(tn)|二位(乙)一V(乙)相对误差准则：ey'<tn)- )S0其中3规定精度的误差量。(3)快速性：若第〃步计算对应的系统时间间隔为4=，”+]—乙，计算机由y(f“)计算小用)需要的时间为Tn,若T^hn称为实时仿真，7〃〈力〃称为超实时仿真Tphn称为亚实时仿真，对应于离线仿真对y=/(歹，〃/)，已知系统变量歹的初始条件>«0)=No，要求V随时间变化的过程》(7),为初值问题，需要用到积分算法。其计算过程:由初始点y(，o)=y0的/(仙M))y(t)=九+1f(t,y)dt由欧拉法得：三％+&•/("，％)对任意时刻1加1K)兄+1=y(tn+i)=yn+(乙+i-乙)•/(。，K)截断误差正比于42^f(t.y)图3.2数值积分法原理常用的梯形法表示： 儿+|=y(tn+,)=y„+^h[f(tn,yj+f(tn+l,y^f(t.y)图3.2数值积分法原理是隐函数形式。预报可采用欧拉法估计初值，然后用梯形法校正：校正公式：y»+\)-yn+)"(如州)+/(3/2)]

预报公式：=yn+"9,北）反复迭代，直到满足卜：；；-v；+i|工£。经典的数值积分法可分为两类：单步法与多步法3.2龙格库塔法.龙格-库塔法基本原理对歹“用）=y（fn）+1'/"，y）dt若令：y”=y（fn）Q"=则有歹（乙+1）三v"+i=V"+Q”Q“的数值求解称作“右端函数”计算问题。在乙）附近展开台劳级数，只保留好项，则有:%=汽+/4,为“+；《?+务1/(15)2oyatat1(15)假设这个解可以写成如下形式：乃=乂）+（。］占+a2k2）h其中 匕=/。0，歹0） ◎=/&+纳，No+&后/）对22式右端的函数展成台劳级数，保留"项，可得：&三/4，氏）+（“普+仇左奈）|,J代入，则有：%=%+4/（柘比）+。%=%+4/（柘比）+。(16)将（16）式与（15）式进行比较，可得：+a2=1,a2bl=1/2,a2b2=1/2于是有四个未知数4，av4，台2,但只有三个方程，因此有无穷多个解。1若限定为=%，则 =7J4=仇=1

计算公式:必计算公式:必=y（）+g（K+与）其中/=/4，%），攵2=/«o+A，必）+左/）h若写成一般递推形式，即为：歹«”+1）三兄+1=匕,+5（/+%2）其中 勺=/（%%）,k2=f（t〃+h,y〃+k［h）3截断误差正比于力，称为二阶龙格库塔法（简称RK-2）。截断误差正比于A的四阶龙格一库塔法（简称RK-4）公式：V（。+1）三/+1=尤+；（K+2左2+2&+左4）

O其中：匕=/（/"，"）k2=f^n+-»yn+-k^&=/（，"+《,X,+京2）左4=/（乙+〃，北+麻3）.龙格一库塔法的特点：（1）形式多样性例：q,a2,",与非唯一解，可以得到许多种龙格一库塔公式：yn+i=yn+k2h,_zh〃，、其中 々=/（%y”） &=/（。+］，尤+］匕）各种龙格--库塔法可以写成如下一般形式：yn+i=yn+hXc-ki/=1i-l其中 号=/⑥+仙，yn+h工/kj）i=L2,…，sj=\式中各系数满足以下关系。1=0/-Iaj— b-i=2,3,•••,sJj=i±c,=iLi=ls称为级数，表示每步计算右端函数/的最少次数。可以证明，1阶公式至少要计算一次，2阶公式^n=2； 4阶公式品=4:依此类推。单步法在计算y”+i时只用到y”，而不直接用V”_1，丁"-2等项。其优点为：存储量减小，可以自启动。可变步长步长〃在整个计算中并不要求固定，可以根据精度要求改变但是在一步中，为计算若干个系数左，则必须用同一个步长力。速度与精度四阶方法的〃可以比二阶方法的〃大10倍，每步计算量仅比二阶方法大一倍，高于四阶的方法由于每步计算量将增加较多，而精度提高不快。3.实时龙格一库塔法实时仿真：要求仿真模型的运行速度往往与实际系统运行的速度保持一致。一般的数值积分法难以满足实时仿真的要求，这不仅仅是因为由这些方法所得到的模型的执行速度较慢，而且这些方法的机理不符合实时仿真的特点。考虑系统 生=/(y，«(/))dt几+】=yn+-(^i+质)RK-2公式如下： “1%2=九+1，〃”+1，4+%〃)一个计算步内分两子步：J时刻：利用当前的%7计算h 计算一次右端函数f需///2

5班々时刻：应计算3尽管此时为+1々已经得到，但%+1则无法得到。（若对与+1也进行预报一一加大仿真误差）。仿真执行延迟h/2一一输出要迟后半个计算步距y〃+i=yn+她实时2实时2阶龙格一库塔法:h后2=/。"+1/2，〃〃+1/2，%+/%）时刻：应计算耳，利用当前的un,“需要A/2；%场々时刻，应计算七此时+1/2已经得到，unh/2也可得到，42的计算就不会引入新的误差。计算一次右端函数/需要〃/2,可实时输出yn采入2）计克为+i并输出3.3线性多步法.线性多步法基本原理基本原理：利用一个多项式去匹配变量若干已知值和它们的导数值。设：，'〃+1 ，，。+1时刻的北，N”+i，,,•，%+"］和/"，"+i'…，9”+"1已知,则该方法可分两步：预报:山歹〃，，〃+1,,♦, 和歹"，9?”+4-1来tlMyn+k,yn+k

校正：若歹"+*也已知，由它们来计算/"+k采用的多项式具有以下形式(勿阶)：产.〃7 (t_/、 m〃⑺这4个一=力尸< ,=0hJ,=om(t_八1 m"⑺=TXU—,— =-iZu/I(=1yhJ,=1其中：4是待定系数，T=tn+k~f.h在，“十斤时刻，2=o,可得到：yn+k=yn+k- (17))=2夕“)=》》=-4,4c丁=-疙即”i=\ \ n J(1)预报公式令m或T,从(18)式得到如下方程组：叫+d，+""••+/若+2/+4*+…+2%以+3町+32御+…+3”力**)=2夕“)=》》=-4,4c丁=-疙即”i=\ \ n J(1)预报公式令m或T,从(18)式得到如下方程组：叫+d，+""••+/若+2/+4*+…+2%以+3町+32御+…+3”力***d：+2d：+…+md](j=l,2 k) (18)二＞〃+"1«=yn+k~2i-yn+k—3=一明+k-l4P+2x +•--+mx2〃jd：=-hyn+k_24P+2x3d；+•--+mx3W_1以=-hyn+k.1 2 tn 〃+£—J■■■11111■-4]r -12 22 23 … 2"'4 匕小213 32 3, … 3"4 匕+一1kk2K…km4T=yn01 2 3 ••md. ~hy..k Jw+i-1012x23x22…mx2m",d —hyA+l n+k-2012x33x22…wx3m-'4+2 -忱*3012xk3xk2…mxk,n~l将其写成矩阵形式：Vpdp-Zp （20）其中上标。表示预报。其解为： dp=（VP）-/ZPd -hy(19)(21)由于为常数阵，其逆存在，z向量中的各元素为已知值，因而d向量的各元素值可计算得到，从而由北+&=4，9”+人=一!4,得到下一时刻的预报值。缺点：只有d（），4是所需要的，其它元素的计算成为多余，得不到4+无，夕”+斤与%，y〃+1，…，歹"+I和夕”，夕〃+1，…，夕〃+"1的显式表达式(22)yn+k=eod(22)定义辅助变量(23)定义辅助变量(23)(24)此式可改写为(y0)‘0P=e。(24)向量。”的元素可划分为两个组 0。=[。3a2,…(25)0'T1a；02x3球03x3204x3305x34_用_0TOC\o"1-5"\h\z1110 012 31 1例：A=3,则(24)式为:1 22 32例：A=3,则(24)式为:1 23 33 3 3x221 24 34 4 4x231 25 35 5 5x24可计算得到： *=(一18,9,10二一9,一18,—3)。'只依赖于k,即先前y〃+4-/和/〃+«_/的个数，而与它们的数值无关。这样，可以预先求解(24)式得到。丹+斤=或"=(0P[ZP (26)从而得到片+*的显式表达式：(27)(2)校正公式预报公式一一显式公式，未包括夕”+人校正：对该预报值应进行校正，即先预报得到歹“‘，然后再用此值推出/”+*由乂，/+1，…，和yn+k，以及%J〃+i，…，％+"i来预报yn+t,可令炉2k~l,从(2)式得到如下方程组:

d:+d：+d；i+d；=y…4+2/+4琥+-一+2〃％：=兄+"24+3/+32/+...+39：=小3■

.4p+24+…+机〃：=-w.+id：+2*2d；+---+Wx2'e%=-hyn+k_2d；+2x3以+…+加x3W-'dcm=-hyn+k3(29)(30)将其写成矩阵形式： Vcdc=Z(29)(30)其中上标c表示校正，可得dc=(vf)-/zc1000 •■■ 0d°yn+kI111-••14122223 •••2md?yn+k-2133233 •..3m"3yn+k—31k左2k3•••kma—yn(31)0123 .•• mdk+i-吼+1012x23x2?•••mx2"idk+2一勿…-2012x33x32-••mx3m-1•一勿…*_d,”_定义：e]=（0,1,0,…,0）为（m+1）xl的列向量，上标7表示转置。将e；左乘（30）式可得:

-hyn+k=e；dc=e](Vc)-Zc (32)定义(33)可改写为仅。)>=e,(34)/二麻M：,…&话埠…闱7(35)例：k=3-111010101av01012012-22x22x320012)33x2°3x32b：0012444x2,4x3、b：0012555x2’5x3’A.0同样，。’只依赖于A,即先前夕“+*_/和的个数，而与它们的数值无关。这样(36)1A (36)1A £二一£(工。)——〃工耳"+3〃7=0从而(37)2.线性多步法误差分析为了便于分析，对预报公式和校正公式，定义统一的表达式：TOC\o"1-5"\h\zk kE 肥…=0 (38)/=0 /=0%=0，显式预报/〃+上，笈0=°，显式预报歹〃+左A=A=••,=A=0时称为后向差分公式(BDF)a。,Bo同时均不等于0时为隐式校正公式，k称为公式的阶次。假设变量各时间的精确值已经得到，将其代入(38)式，可得：Z叼心一切)一应尸M-讷) (39)i=O 1=0在乙+左附近，将每个函数展开成泰勒级数:

yQn+k—曲)-y^n+k)-与y(fn+k)+粤产“〃+%)+…yqn+k一诙)=/«〃+*)—卓»«〃+〃)+粤y^n+k)+…(40)对所有〃i=0,1,2,…,公，将(40)式代入(39)式，合并同类项，可得Cj(g+〃c+9(心)+…+g严(-)+.••(41)其中Co=CXq+。1+CC-,+…+cxkG=(-1)'{=(&+2a2h hkaj++电+。卢 1-/3k)}C2=(-l)2{=(q+2辿+…+EaJ+%(仪+2A4 1-k/3k))C3=(-1)'{+@+23a2♦ 卜/4)++(4+2%2T 卜/月)}■*

■G=(-1广也(6+2川见+…+%")+$(4+2*,民+…+亡0)}(42)如G),£,G,g均为0,则称为p阶公式t卬«诙)-应月M*-%)=〃"Cj"(a)+…(43)1=0 /=0预报与校正公式的截断误差系数可以用(44)式来计算，即以=(一，"{,(3+22+…+/％)+呆片+29+…+上口)}(44)3.稳定性分析仿真方法选择的基本要求：仿真计算不改变原系统的绝对稳定性。dy!dt=jLty,=a+且Re^=a<0原系统是稳定的。观察欧拉法仿真递推公式/+i=y"+勿，儿)故有 y„+i=yn+丽(45)%(炉0,1,2,…)为它的一个仿真解，设北+£“为其准确解，即TOC\o"1-5"\h\z(%+1+£〃+1)=(兄+£〃)+睢(y〃+J) (46)用(46)式减去(28)式，可得：£”+i=£”+hjLlEn即 £”+1-(1+ =0 (47)特征方程为z—(1+"?)=0 (48)显然，为了使扰动序列&不随〃增加而增长，必须要求：\\+^lh\<\ (49)我们称它所对应的域就是该算法的稳定域：K211/a|,即分小于等于系统时间常数的两倍。确定数值积分法稳定域的一般方法测试方程：dyldt=juy,〃=a+场数值积分公式片+1=M砧)•、” (50)其中p(6〃)是一个关于〃4高阶多项式函数，则只有当I2(〃口)|<1时，算法才稳定。3.4离散事件系统仿真基础3.4.1基本概念离散事件系统：系统中的状态只是在离散时间点上发生变化，而且这些离散时间点一般是不确定的。

例如：单人理发馆系统，设上午9:00开门，下午5:00关门顾客到达时间一般是随机的，为每个顾客服务的时间长度也是随机的。系统的状态：服务台的状态（忙或闲）、顾客排队等待的队长。状态量的变化也只能在离散的随机时间点上发生。在离散系统仿真中，有以下术语需要定义：（1）实体分为两大类：临时实体及永久实体临时实体：在系统中只存在一段时间的实体。这类实体由系统外部到达系统，通过系统，最终离开系统。永久实体：永久驻留在系统中的实体。只要系统处于活动状态，这些实体就存在，或者说，永久实体是系统处于活动的必要条件。临时实体按一定规律不断地到达（产生），在永久实体作用下通过系统，最后离开系统，整个系统呈现出动态过程。（2）事件引起系统状态发生变化的行为。从某种意义上说，这类系统是由事件来驱动的。顾客到达事服务开始事服务结束事图3.4事件、活动、进程三者关系示事件表：实现对系统中的事件进行管理，表中记录每顾客到达事服务开始事服务结束事图3.4事件、活动、进程三者关系示系统事件：系统中固有事件，“程序事件”，用于控制仿真进程。（3）活动用于表示两个可以区分的事件之间的过程，它标志着系统状态的转移。顾客的到达事件与该顾客开始接受服务事件之间可称为一个活动（4）进程进程由若干个事件及若干活动组成，一个进程描述了它所包括的事件及活动间的相互逻辑关系及时序关系。（5）仿真钟离散事件动态系统的状态本来就只在离散时间点上发生变化，因而不需要进行离散化处理。由于引起状态变化的事件发生时间的随机性，仿真钟的推进步长则完全是随机的两个相邻发生的事件之间系统状态不会发生任何变化，因而仿真钟可以跨过这些“不活动”周此仿真钟的推进呈现跳跃性，推进速度具有随机性。时间控制部件是必不可少的，以便按一定规律来控制仿真钟的推进。（6）统计计数器某一次仿真运行得到的状态变化过程只不过是随机过程的•次取样，它们只有在统计意义下才有参考价值。 在仿真模型中，需要有一个统计计数部件，以便统计系统中的有关变量。离散事件系统仿真研究的•般步骤：（1）系统建模（2）确定仿真算法（3）建立仿真模型（4）设计仿真程序（5）仿真结果分析4.2离散事件系统仿真策略离散事件系统仿真策略主要分三种：事件调度法、活动扫描、进程交互法。1、主要概念（1）成分（Component）：相应于系统中的实体，用于构造模型中的各个部分，可分为两大类：主动成分（Active-typeComponent）：可以主动产生活动的成分，如排队系统中的顾客，它的到达将产生排队活动或服务活动。被动成分（Passive-typeComponent）：本身不能激发活动，只有在主动成分作用下才产生状态变化。（2）描述变量：成分状态、属性的描述。（3）成分间的相互关系：描述成分之间相互影响的规则。在一个模型中，主动成分对被动成分可能产生作用，而主动成分之间也可能产生作用。c=｛a,,a2,…，a”）成分集合，a,是第/个成分分量（1<i4a=｛。,.a：；…,a.｝主动成分子集，a,是第j个主动成分分量（i<JW相，相<〃）。a、…,a,） 被动成分子集，a.是第〃个被动成分分量（14左<一个模型中， n=m+lS为所有成分的状态变量，值域为P=｛p/,P-f…，pn] 参数（属性）集合。成分。的状态下一发生变化的时刻，值域为｛H1,8）｝Da（s）成分a在状态变量值为s时的条件是否满足Da（S）-true,表示满足，d*（s）=false表示不满足。TIME为模型仿真钟的值，值域为怵】6}。2、事件调度法（EventScheduling）事件调度法基本思想：统状态的变化，按时间顺序确定并执行每个里件发生时有关的逻辑关系。・r事件均放在事件在中。模型中有个时间控制成分，该成分从事件表中选择具有最早发生时间的事件，并将仿真钟修改到该里件发生的时间，再调用,j该事件相应的事件处理模块，该力事件产生为■策略的非形式描述：成分集合={ava2,--an}.主动成分集CA={ava2,---am}被动成分集CP={am+vam+2,--an/描述变量：描述每一主动成分aeC\的变量，a的状态$a值域§a，$a下一变化时刻的时间变量L描述每•被动成分aeCp的变量，a的状态Sa，值域Sa（被动成分的状态变化只有在主动成分作用下才能发生，其发生时间由主动成分来确定，因而不需要时间变量。）描述所有成分的属性的变量：参数集合仁山，Pz,…,P，］成分间的相互关系•每个主动成分aeC］的影响受主在a作用下其状态变化的描述，称为事件处理流程；•各成分处理的优先级，即同时发生时的处理顺序（解结规则）。在事件调度法中，一般主动成分也同时具有被动成分属性，以便接受其它主动成分的作用。事件调度法策略描述执行初始化操作，包括：置初始时间t=to,结束时间tg—te事件表初始化，置系统初始事件成分状态初始化S=((Sa、，ta、)L，(Sa，/a)>Sa%)%”S/w+l操作事件表，包括取出具有s)=min{(a|aeCA}事件记录修改事件表推进仿真钟TIME=t(s)While(TIME<=ro0)则执行Case根据事件类型/7=1执行第1类事件处理程序*(*第/.类事件处理程序对成分的状态变化进行建模，而且要进行统计计算)户2执行第2类事件处理程序2=m执行第m类事件处理程序endcase取出具有=加比"a\aECa}事件记录**(**若具有=加a£C//事件记录有若干个，则按解结规则处理)置仿真时间TIME=t(s)endwhile***(***该策略描述未包括对仿真结果的分析等)2活动扫描法(ActivityScanning)羽件调度法中仿真钟的推进仅仅依据CA/准则，问该事件发生的任何条件的测试必须在该事件处理程序内部去处♦甲如果条件满足,该事件发生，否则,则推迟或取消该事件发生。从本质上来说，事件调度法是一种“预定事件发生时间”的策略。这样，仿真模型中必须预该策略对于活动持续时间确定性较强■可以是■某种分布的随机变量)的系统是比较方便的。,与事件的发生不仅与时间有关，而且与其它条件行关，即只有满足某些条件时才会发生。在点则表现出来了.由于这类系统的活动持—的不确定性，因而无法预定活动的开始或终止时间。活动扫描法的基本思想是：用活动的观点建模。系统由成分组成，而成分包含着活动，这些活动的发生必须满足某些条件；每个主动成分均有个相应的活动，例程：仿真过程中,活动的发生时间也做为■之一，—之其它条件具有更高的优先♦.设。/S)表示成分a在系统状态S下的条件是否满足(o/s)=5ue则表示满足，々/se则表示不满足)，J表示成分a的状态下一发生变化的时刻，活动扫描法每一步要对系统中所有主动成分进行扫描，当：GW仿真钟当前值TIME,且(ii)Da(S)=true^,执行该成分。的活动子例程。所有主动成分扫描一遍后，则又按同样顺序继续进行扫描，直到仿真结束。显然，活动扫描法由于包括了对事件发生时间的扫描，因而它也具有事件调度法的功能。实现措施：(1)设置系统仿真钟TIME与成分仿真钟系统仿真钟表示系统仿真进程的推进时间，而成分仿真钟则记录该成分的活动发生时刻，两者的关系可能有三种情况：ta>TIME表示该活动在将来某一时刻可能发生；,a=TIME表示该活动如果条件满足则应立即发生；<TIME表示该活动按预定时间早应发生，但因条件未满足，到目前为止实际上仍未发生，当前是否发生，则只要判断其发生的条件。(2)设置条件处理模块该模块用于测定Da(S)的值及系统仿真钟与成分仿真钟之间的关系，记FUTURE(S)={^Va>TIME},PRESENT(S)={口'a=TIME},PAST(S)={ a<TIME}该模块将满足下列条件：(i)&£PRESENT(S)uPAST(S),且(ii)0a(S)=true的成分置于可激活(ACTIVATABLE)的成分集合中，即JaePRESENT(S)uPAST(S)ACTIVATABLE^=[\Da(S)=true此时，系统仿真钟不推进，仅仅处理成分活动，包括修改成分仿真钟。如果可激活的成分集合为空，则将系统仿真钟推进到下一•最早发生的活动生成时刻，即TIME=min('aFUTURE(S)策略描述：执行初始化操作，包括：置初始时间t=t0,结束时间忆二。设置主动成分的仿真钟力(i=l,2,…,m)ui成分状态初始化：S=((S/),•••(5,t),5，…S)设置系统仿真钟TIME=t"While(TIME<C),则执行扫描for户最高优先数到最低优先数将优先数为J的成分置成iIf(ta<TIME且Da(S)=true)执行活动子例程i退出，重新开始扫描end//'end/brTIME=min(^aFUTURE(S))endwhile策略核心：建立活动子例程模型，包括此活动发生引起的状态变化（自身的），对其它成分的状态产生的作用等等，条件处理模块：策略实现的本质，它相应于事件调度法中的定时模块。4、进程交互法(ProcessInteractive)进程交互法采用进程（Process）描述系统，它将模型中的主动成分所发生的事件及活动按时J序进行组♦.从