计算机模拟技术 27页

计算机模拟技术课程名：计算机模拟技术计算机模拟是在科学研究中常采用的一种技术，特别是在科学试验环节，利用计算机模拟非常有效。所谓计算机模拟就是用计算机来模仿真实的事物，用一个模型(物理的-实物模拟；数学的-计算机模拟)来模拟真实的系统，对系统的内部结构、外界影响、功能、行为等进行实验，通过实验使系统达到优良的性能，从而获得良好的经济效益和社会效益。计算机模拟方面的研究始于六十年代，早期的研究主要用于国防和军事领域(如航空航天、武器研制、核试验等)，以及自动控制等方面。随着计算机应用的普及，应用范围也在扩大，现在已遍及自然科学和社会科学的各个领域。在农业方面，我国从80年代开始进行作物生长发育模拟模型和生产管理系统的研究，目前有一定基础的：在小麦方面有北农大、中科院；棉花方面有中国农业大学、中国棉花所；水稻方面有江西农科院；在土壤水份、水资源及灌溉方面西北农业科技大学。目前影响较大的有比较成形的有江苏省农科院。目前的主要成果有：我国主要农作物栽培模拟优化决策系统RCSODS(水稻)和WCSODS(小麦-江苏省农科院)、MCSODS(玉米-河南省农科院)、CCSODS(棉花-中国农业大学)等。计算机模拟特别适合于实验条件苛刻、环境恶劣(如真空、高温、高压、有毒有害的场所)、试验周期长，花费大的场合。农作物的生产系统就很适合于计算机模拟：农作物的生产受各种条件的影响，不同作物、不同品种也有差异。比如，要想提高一种作物的产量，就先要作试验，通过试验了解这种作物的特性：抗旱性、耐寒性、对氮、磷、钾哪种肥更有效等。但农业的田间实验不能保证精度(除人为可控条件外，还有许多随机因素)、周期长(周期一年)，耗费大。可通过计算机模拟来实现：先建立这种作物生产系统的数学模型(依靠专业知识或试验数据。一般来说，诸如作物产量和农业环境的关系可用微分方程或其它方程来描述)，通过计算机模拟来找出这种作物的生长与农业环境相互作用的关系，以及各种条件之间的协迫情况。不仅可大大节省实验经费、加快研究进度(周期一年的实验结果几秒钟内即可得到)，这种模拟软件的开发还可与农业生产管理系统，决策系统相联系，实现对农作物生产的预测、分析、调控、设计的数字化和科学化。作为一门课程，不是研究某个特定系统的模拟问题，而是了解计算机模拟的一般过程、基本原则，掌握基础知识，掌握建模及动态模拟的一般方法。第一章计算机模拟概述1.1计算机模拟技术●研究对象在一个计算机模拟问题中，我们研究的对象是一个系统。系统：一些具有特定功能的、相互之间按一定规律联系着的实体的集合。如作物的生产系统可看作由作物、环境、技术、经济等要素构成的。各要素之间相互影响、相互联系，称为系统的相关性；一个系统是一个整体，整体内的各个部分不能分割，各因素之间必须相互协调，不能在任何一个环节出问题，才能使系统达到优良的状态，称为系统的完整性。●目标计算机模拟的目标是了解系统的各个实体之间的相互制约关系，从而使系统在预定的目标下达到最优和完善。如在作物生产系统中，怎样控制、实施各水、肥、栽培技术等，从而使产量最高，以获得最优的经济效益。●方法模拟的方法是先建立系统与环境相互作用的数学模型，用数学模型来类比、模仿现实系统(一个数学模型就是从数学上表达系统各因素之间的数量关系，或各因素之间协调的规则；从整个模拟过程来看就是一个算法，或一系列数据，这些数据综合描述一个系统过程或现象的重要行为\n)，然后在数学模型和对系统深刻了解的基础上，开发模拟软件，用影响系统目标的因素作为输入，通过计算机技术来表达系统各因素作用的状态。从数学的角度来看，模拟的过程就是对数学模型求解的过程，并把系统过程演示出来。●基础知识可见，对一个系统进行计算机模拟，(1)要对模拟对象有深刻的了解。如：一个公交车调试系统，要编制一个好的调度程序，必须先对现行系统作周密的调查，搞清哪些是影响调度程序的成分及实体，如现有车辆数、每车载客数、每趟车花费的时间、沿途乘客的密集程度、乘客的一般去向、乘客高峰期的人数等。只有经过周密的调查研究，才能形成一个完整的模拟系统；作物生长模拟系统中，也要搞清影响作物生长的各因素，具备丰富的专业知识，否则不能建立起精确的模型。在对一个系统的动态特性不完全清楚的情况下，有必要通过实验获取数据，以用于数学模型的建立。(2)要有数学知识。一般研究的系统较复杂，不能用简单的函数或方程来描述，要综合使用各种数学方法，才能使模型准确、可靠。(3)计算机知识和编程技巧。软件应完整地实现系统的数学描述，输出应直观、形象，如三维可视化输出等。软件的开发是项目的重要工作，也直接影响模拟的效果。软件的编写可使用任何编程语言，如C、VB、Java等。专门的语言如GPSS，GPSS是面向对象问题的离散事件的专用模拟语言，优其适用于排队系统。1961年，IBM公司发表GPSS的第一个版本，后来又有其它公司的各种版本。标准的版本有52个模块，每个模块用特定的名称和图形来表示其功能。现在一般常用的语言都有模拟库(已编好的用于实现模拟功能的函数)专门用于模拟软件的开发。1.2系统的分类可模拟的系统各种各样，不同类型的系统用不同的模型来描述。系统的分类方法很多，重要的分类方法是按系统的状态是否随时间变化来分：一个行为与时间有关的系统称为动态系统待研系统：静态系统：系统的行为与时间无关；用静态模型来描述，一般为数学方程、逻辑表达式等。如，电路的布尔表达式；电路中电压与电流的关系；系统的稳态解公式等动态系统：连续系统：系统的状态随时间连续变化；常用微分方程来描述；方程对所有时间点有效。如，卫星运行轨道，作物生长量等确定性系统：系统的输出完全由其输入来描述，即系统输入与输出按某种规则一一对应集中参数模型：用常微分方程来描述，即方程中的导数不是偏导数分布参数模型：用偏微分方程来描述，但一般用集中参数模型近拟表示随机系统：系统的输出是随机的，有规律的存在一族随机变量，且随机变量序列与时间有关(随机过程)。在确定性系统的模拟中使用随机变量的研究方法称为蒙特卡罗方法离散系统：系统状态的变化只在离散的时间发生；动态方程只在离散点上有效；一般为随机系统。如，库存问题；企业的管理系统等离散时间系统：时间步长固定；常用差分方程来描述离散事件系统：用事件来表示系统在时间间隔内的变化；常用概率模型来描述1.3建立数学模型对一个系统，确定其类型后有助于选择合适的方法建模，建模的一般步骤可分两个大的阶段：(1)实质内容模型阶段：首先对模拟对象进行调查(了解系统，搜索模拟所需的信息)、实验(\n参数估计等统计推断方法，确定参数及参数的敏感性)、分析(将信息分类、量化，确定描述系统的规则)，尽可能全面地掌握系统的基本特性、运动规律以及中间状态(最好是有对系统有深入了解的专业人员参与)，通过分析和逻辑推理，揭示系统内的规律。(2)形式数量阶段：在调查、实验、分析的基础上，进一步揭示系统内部的数量关系，并对其进行数学处理，即对系统用数学形式来描述：用变量描述系统状态，用各种数学方程定量表示各变量之间的相互联系，用递归方程描述系统状态的发展趋势。多数情况下，建立一个好的(能真实的描述系统，有代表性，能准确的模拟系统的数量信息，与实际系统较吻合)数学模型不易，优其是复杂多变的系统或系统本身的特性尚不完全清楚的情况(如对农作物开发新品种)。所以在数学模型建立后，必须进行模拟验证，如与真实系统相差较大，则要重新建模或修改模型。所以一个好的数学模型必须经过多次模拟，不断修改、完善，才能得到。一个数学模型是描述系统行为的一个算法或一系列方程，工程中对系统建模应先建立系统的需求规格说明，在制定模拟规划之前予以充分讨论，过程中需考虑以下因素：(1)模型中需考虑哪些有效的因素：一个系统可能有不同的行为(如一个作物生产系统中有作物的长势、产量、质量、病虫害等)，但不一定所有这些因素都要建立在模型中。模型中需要考虑的因素是：真正能简化系统的建模；对系统易于建模、测试和维护；使用较少的计算资源；对研究的系统有直接作用的。(2)建模细化到什么程度：根据系统的需求确定细化的程度，是只须建立一个简单的模型，还是要对系统行为精确描述。在模型的准确性和花费之间求得平衡。(3)与系统有相互作用的哪些外部环境考虑在建模中：如在作物生长模型中，气候、水、肥等。(4)在建模中采取什么技术：首先，是基于物理的方程，还是基于测试数据。如果基于物理的方程，是用微分方程，还是差分方程，考虑不考虑随机因素等。这个问题常由专业人员根据系统的专业知识来决定；如果基于试验数据，则建立经验方程。(5)建模时必须获得什么样的数据：如林木生长模型中，需要胸径、树高、材积等数据。另外还需考虑这些数据的方便的输出形式，以及模型需要的计算资源，如模型需要占用的内存、磁盘空间、消耗的CPU时间等。(6)建模和测试模型需多长时间，多少人力、物力、财力：随着模型复杂性的增加，成本也会增加。一般可从简单开始，随着应用逐步完善。(7)如何验证和确认模型：必须确认建立的模型被正确实现，模型所描述的行为与真实系统匹配到可接收的程度，才能有价值。那么以什么标准来衡量。现在已有验证、确认与签定(VV&A)技术来证明和验证模拟的精确性。以上问题应列在系统需求说明书中，并应用于最高层次：在实际的模拟问题中，可能将系统分解为子系统、组件，在对各组件、子系统建模时仍须遵从以上规则。初始计划编制评估测试分析设计实现需求计划编制部署图1-1迭代开发模型模拟步骤：(1)明确系统(2)建立模型(3)模型变换(4)软件设计开发(5)测试检验(6)评估、对结果的评价和分析一个模拟项目中各项工作的过程应该是一个迭代过程，如图1-1所示。下面通过实例来说明模拟过程。1.4应用举例1.4.\n1两物体追逐问题。设有一架歼击机追踪一架敌方轰炸机，假设两机相距10公里以内可实施攻击，且须在12分钟内完成追击任务，否则认为追击失败。设两机初始位置如图1-2所示。问题：对轰炸机的任一条特定航线，模拟歼击机的追击过程。(XB(0),YB(0))D(t)VF=200204060801001201401601802005025图1-2两机追踪模拟分析：在这个系统中，是对指定的轰炸机的一条航线而言，模拟歼击机的追击情况：歼击机按什么航线飞行，何时完成追击任务。能否歼灭敌机由在规定的时间内两机随时间变动的距离而定，所以实施功击的时间、距离是要输出的数据。在时刻t两机距离由歼击机在t时的位置决定，而t时的位置依赖于其速度和航向。为使模型简单，作如下简化：(1)设两机在同一平面飞行。于是三维问题转化为二维问题。(2)设歼击机的速度(必须考虑的因素)VF是常数(20km/mim)。变速须用微分方程描述，而常速即可用一般方程表达，求解更简单。(3)设歼击机的航向(必须考虑的因素)在△t(设1分钟)改变一次，而在1分钟以内操作不变。这样曲线运动即变为折线运动。(4)轰炸机航向(航线)可任意，比如现给定一条航线，如表1-1所示。表1-1轰炸机航线t012345678901112XB(t)809099108116125133141151160169179180YB(t)50484541353227212225293033(5)歼击机初始位置：XF(0)=0，YF(0)=0(初始条件)建立数学模型：设在t时刻两机位置为(XF(t),YF(t))、(XB(t),YB(t))，两机连线与水平线夹角为θ，则1分钟后歼击机的位置为：XF(t+1)=XF(t)+VFCosθ(1.1)YF(t+1)=YF(t)+VFSinθ(1.2)由(XF(t),YF(t))计算(XF(t+1),YF(t+1))时涉及角θ，从图中看出：Sinθ=[YB(t)-YF(t)]/D(t)(1.3)Cosθ=[XB(t)-XF(t)]/D(t)(1.4)而D(t)=SQR{[YB(t)-YF(t)]2+[XB(t)-XF(t)]2}(1.5)式(1.5)-(1.1)即是这个问题的数学模型。先算出两机之距离，不断判断是否在12分钟内到达可攻击的距离之内。程序流程图如图1-3所示。此例是连续系统，确定性模型，用一组方程来描述。下例是一个离散事件系统，是随机系统，用概率模型来描述。\n开始输入轰炸机航线数据t=0计算D(t)D(t)≤10输出t、D(t)输出目标逃脱信息t>12t=t+1计算新位置结束图1-3追击模拟流程图VB程序见实例。图1-4售货进程库存01234t300200100定货点1.4.2库存问题。商业部门为了合理利用有限的流动资金，每项商品都要在库存与销售之间求得平衡：库存量太大会增加管理费用、积压资金；库存量太小又可能造成缺货，也会造成销售损失、信誉损失。所以，当库存量不满足某一时段的顾客需求时，就要到厂家订货。这就需要采取一种策略：当库存量(比如布匹)降到P匹布时(称为定货点)则向厂家订货，定货量为Q匹(称为定货量)。假设现在有5种策略(方案)。从中选择一种使花费最少。表1-2库存策略策略PQ11251502125250315025041752505175300已知条件：(1)从发出订货单到收到货物需3天，即第i天订货，第i+3天收到。(2)每匹布的保管费为0.75元，缺货损失为1.8元/匹，订货费(包括手续费、采购差旅费及其他费用)为750元。(3)需求量为一个0-99之间的均匀分布随机数。(4)原始库存量为115匹，并设第一天没发出订货。\n分析：库存问题是商业上的一个重要问题在数学上有专门的模型研究，存储论也是运筹学的一个重要分支。库存问题用计算机模拟最合适，若通过销售来找最优方案，势必造成经济损失。这是典型的离散事件系统，用概率模型来描述。这里已给定概率模型：X~U［0,99］即密度函数为：p(x)=1/990≤x≤990x<0分布函数为：F(x)=x/990≤x≤99开始输入：预定到货日期、最初值、总费用、日期、原库存量等今天是否到货期产生随机需求量库存+Q需求量≤S(库存量)总费用+缺货损失算新库存量、保管费预期库存≤P预定到货Q，到货日期当日+3日期+1満180天输出结束图1-5库存问题模拟程序框图1x>99以下框图对每种策略模拟180天，选出效益最好的方案，如图1-5所示。C程序见实例。以上两个例子一个连续系统，状态随时间连续变化，用方程来描述；一个离散系统，到货和销售都按一些离散的步骤进行，存在随机性，只能用概率模型来描述。但解决问题的过程、分析方法类似：先分析系统，使对系统有充分的了解；建立数学模型，设定一些初始条件(如t=0时的状态，开始时的库存量)；系统状态的变化对应一组方程或一组规则，随着时是的变化，系统状态改变，当一个周期结束时，收集模拟过程的统计数据(即问题的解)。如果程序设计的好，就会使模拟非常逼真，就象一个真正的系统在演示一样。\n从理论上讲，任何问题都可用计算机模拟。计算机模拟具有经济、安全可靠、周期短等优点。对任何问题，只要建立起数学模型，改变参数值及变量值，就可模拟各种情况下的系统运行情况，从模拟输出的结果，可分析系统内各因素的权重及其制约关系，帮助决策者作出合理的决策，克服盲目性，使系统在实际运行过程中取得最好的效益。1.4.3传染病传播问题。假设某一地区有一种传染病正在流行，那么政府、医务部门就要采取措施来控制这种疾病的传播，要使采取的措施能够有效的控制传染速度，就必须先搞清被传染的人数跟哪些因素有关，被传染的人数是一个什么样的发展趋势，从而来有效的预测和控制，下面建立描述被传染人数的数学模型。传染病的传播涉及因素很多，不可能通过一两次简单的假设就能建好完善的数学模型。这里的作法是：先作出最简单的假设，看会得到什么样的结果，然后对不合理的地方再行修改，逐步得到较满意的模型。先讨论一个粗略的模型。模型I：假设：(1)每个病人在单位时间内传染其他人的人数为一个常数k0。(2)一人得病后，经久不愈，该人在传染期内不会死亡。记时刻t的得病人数为y(t)，开始模拟时有y0个传染病人，则在△t时间内增加的病人人数为y(t+△t)-y(t)=k0y(t)△t由导数定义得(在假设(1)、(2)下的数学模型)：dy/dt=k0y(t)(1.6)y(0)=y0初值问题(1.6)的解为：y(t)=y0ek0t(分析：)这个结果表明，得病人数将按指数形式无限增加，当t→∞时，y(t)→∞，显然与实际不符，说明上面的假设条件不合理。事实上，一个地区的总人数大致可视为常数(不考虑疾病传播期间出生的、迁出的、死亡的)，所以一个病人在单位时间内能传播的人数k0是在改变的：在初期，k0较大，随着病人的增多，健康人数的减少，被传染的机会也将减少，所以k0将变小。对上述模型进行改进：记t时刻的健康人数人S(t)，当总人数不变时，k0随S(t)减少而减少。模型II：假设：(1)总人数为常数n，t时刻的健康人数为S(t)，得病人数为Y(t)，则Y(t)+S(t)=n(1.7)(2)单位时间内一个病人传染的人数与当时的健康人数成正比，比例系数为k(医学上称为传染强度)(3)同模型I的假设(2)。由假设(2)，(1.6)式中的k0应为kS(t)，即：dy(t)/dt=kS(t)Y(t)(1.8)y(0)=y0将(1.7)式代入(1.8)式，得(上述假设下的数学模型)：dy(t)/dt=kY(t)(n–Y(t))(1.9)y(0)=y0初值问题(1.9)的解为：Y(t)=n/[1+(n/y0-1)e–knt](1.10)(分析：)由(1.10)式得：dy/dt=[kn2(n/y0-1)e–knt]/[1+(n/y0-1)e–knt]2(1.11)(1.10)式是被传染人数随时间变化的关系；(1.11)式是被传染病人的变化率与时间的关系，如图1-6和图1-7所示。\ndy/dt0tyny00t1t图1-6病人人数变化曲线图1-7病人变化曲线这个模型可预报疾病高峰到来的时间：令d2y/dt2=0，得极大值点：t1=ln(n/y0-1)/kn(1.12)由(12)式可知，当传染强度k或人数n较大时，t1变化较小，表明传染高峰到来的快，这与实际情况相吻合。但由(1.10)式知当t→∞时，y(t)→n，这意味着最终人人都被传染，这与实际不符，原因在于假设(3)不合理。再改进：可将人员分为三类：第一类为传染者(y)；第二类为易受传染者(s)，即这一类是非传染者，但能得病而成为传染者；第三类为除前两者之外的人(r)，包括患病死去的、痊愈后具有长期免疫力的、被隔离的。用y(t)、s(t)、r(t)分别表示这三类人数。模型III：假设：(1)总人数为常数n，则：y(t)+s(t)+r(t)=n(1.13)(2)同模型II的假设(2)(3)单位时间内病愈免疫的人数与当时病人的人数成正比，比例系数为m(恢复系数)由假设(3)，有dr/dt=my(t)(1.14)由于引入了r(t)，则模型II的方程(1.8)应改为：dy/dt=kS(t)y(t)–dr/dt(1.15)(1.15)式表示单位时间内病人人数的增加应等于被传染的人数减去病愈的人数。从(1.13)-(1.15)中消去dy，并设S(0)=S0，r(0)=r0得dS/dt=-kS(t)y(t)dr/dt=my(t)y(t)+s(t)+r(t)=n(1.16)S(0)=S0r(0)=r0模型(1.16)较好地描述了传染病传播问题。通过以上实例，对数学模型的建立就有了一个基本的思路，对计算机模拟技术、模拟的过程、问题的方法步骤也有了一个概括的了解。后面两章下分别对连续系统和离散系统讨论基本的模拟方法。习题：(1)对例1.4.1中的飞机追击问题采用极坐标(r,θ)，相应的方程为：dr/dt=VBCosθ-VFr·dθ/dt=-VBSinθ\nθ图1-8球摆其中r为两机距离，θ为两机连线的角度。编程模拟两机追击情况。(2)球摆：如图1-8所示，设绳长为L,夹角为θ，球质量为m，初始速度θ(0)=0，初始偏角θ0，确定摆动周期。提示：影响摆动周期的因素：①摆球重力；②摆球质量；③摆球尺寸相对于绳很小，故可视为质点建模；④绳子质量忽略；⑤磨擦力忽略；⑥空气阻力不考虑。建模：影响摆球运动的重力(使摆球回到中心位置的力)F的方向与绳子垂直；对绳子产生拉紧力的重力与绳子平行，它不影响摆球运动，忽略。F=-mgSinθ，将牛顿定律F=ma代入上式得：a=-gSinθ，其中a是摆球的切线加速度，由于L=θ"，故得运动方程：θ"=-L/g·Sinθθ(0)=θ0，θ´(0)=0\n第二章连续系统的计算机模拟技术连续系统的状态随时间连续的动状变化，这种变化依赖相关的因素，且遵从一定的规律，所以一般来说可用数学方程描述(代数方程、微分方程，此外还有传递函数、结构图及状态方程等)。问题：对数学模型怎样求解，求解过程也就是模拟的过程，即程序的算法。对前面的飞机追击的问题，其模型是一组代数方程，而且是显式的，所以，作法是：假设1分钟改变一次航向，每隔1分钟计算一次系统的状态，而且由模型知道了歼击机前一分钟的状态就可以算出后一分钟的位置。这样，通过迭代，就可求出问题的解。对一般的模拟问题，模型用微分方程表达(方程中除含有自变量外，还有导数或偏导数—分布参数型可通过一些变换转换为集中参数型)，而对一个常微分方程(只有线性的或几种特殊的能求出通解及特解，即解析解)一般来说求解析解是不可能的，(即使能求出解析解)在计算机上求解要用数值积分法：将连续的系统离散化，把微分方程转化为差分方程，通过迭代运算，求出问题的数值解。通过正确的控制步长、选择合理的算法即能达到要求的精度，这就是连续系统模拟的主要技术。数值积分法种类很多，本章介绍几种简单常用的方法。2.1欧拉(Euler)法f(t,y)f(tn,yn)y=0…tn-1tntn+1tn+2…tf(t,y(t))图2-1欧拉法的几何意义设连续系统的数学模型为dy/dt=f(t,y)(2.1)y(0)=y0为了模拟系统状态y随时间t的变化，需求解微分方程(2.1)的数值解(不是解析解)，为此，把模拟周期分为若干小区间，比如分为N个相等的小区间，如图2-1所示。只须在每个离散的点上计算系统的状态。在区间(tn,tn+1)上求积分，得y(tn+1)-y(tn)=积分的几何意义是小曲边梯形的面积。如果小区间取的足够小，则在区间(tn,tn+1)之间的f(t,y)可近似的看成常数f(tn,yn)，这样可用小矩形的面积代替小曲边梯形的面积，于是得在tn+1时的积分值为y(tn+1)≈yn+f(tn,yn)·h=yn+1(2.2)其中h=T/N，将(2.2)式写为差分方程形式，得yn+1=yn+f(tn,yn)·hn=0,1,2,…(2.3)这就是欧拉公式：它是一个递推的差分方程，任何一个新的数值解yn+1均基于前一个数值解yn以及导数值f(tn,yn)求得，只要给出初始条件y0及步长h，就可算出y1，由y1可算出y2，如此迭代计算y3，y4，…，直到满足所需计算的范围。欧拉法也叫折线法，特点是方法简单，计算量小，但计算精度也底。扩展I：改进的欧拉法(梯形法)为了提高计算的精度，可对欧拉法进行改进：从几何意义上看到，用小梯形的面积来代替小矩形的面积，必能提高精度。这样(2.3)式即可写成：yn+1=yn+h/2·[f(tn,yn)+f(tn+1,yn+1)]=yn+h/2·[fn+fn+1](2.4)但(2.4)式是隐式公式，即公式右端含有yn+1，而这是未知的待求量，故梯形法不能自行启动运算，要借助于其它算法：如用欧拉法算出ypn+1作为梯形法中ycn+1的预估值，再进行计算，这样，公式可写为：\n预估：ypn+1=yn+f(tn,yn)·h校正：ycn+1=yn+h/2·[f(tn,yn)+f(tn+1,ypn+1)](2.5)=yn+h/2·[fn+fpn+1](2.5)式也称为预估校正法，计算量稍大(需要附加的计算)，但精度要比欧拉法高，稳定性好。扩展II：多个输入的多阶系统在(2.1)式所表示的数学模型中，y是系统的状态，称之为状态变量，它随时间t而连续变化。(2.1)式是最简单的情况，只有一个状态变量，而它直接依赖时间而变化，在实际模拟中，状态变量可能不只一个，而影响系统状态的因素更多，这些因素是随时间变化的同，而系统状态的改变也正是由于这些因素随时间变化而变化。如在林木的生长系统中，要考察的状态可能有树高、胸径、材积等；在作物的生长系统中，考察的状态可能有株高、叶子的片数、叶子的宽度等，而影响系统状态的因素如水、肥、气候、土质、病虫害等，而这些因素都随时间而变化，如水份会随时间而被蒸发，肥会随时间被作物吸收，气候等自然环境更是随时间变化无常。所以与(2.1)式相应的数学模型应为：dyk/dt=f(x1(t),x2(t),…,xm(t),y1,y2,…,yq)k=1,2,…,q(2.6)yk(0)=yk,0其中yk为系统状态变量，xi为系统输入信号，这样的系统称为具有m个输入、q阶系统。(2.6)式是q个微分方程的方程组，q个初始条件。可以把(2.6)式表示成向量的形式，就和(2.1)式在形式上一致了：记Y为q个状态变量和向量，X为m个输入信号向量，即：Y=(y1(t),y2(t),…,yq(t))，X=(x1(t),x2(t),…,xm(t))则(2.6)式写为：dY/dt=f(X(t),Y)(2.6)’Y(0)=Y0对于(2.6)式所表达的数学模型，与(2.3)式相对应的欧拉公式为：yk,n+1=yk,n+f(x1(tn),x2(tn),…,xm(tn),y1,n,y2,n,…,yq,n)·h(2.7)k=1,2,…,q或写成向量的形式：Yn+1=Yn+f(Xn,Yn)·h(2.7)’例2.1设两种物质A和B合到一起产生化学反应，生成新物质C，假设1克A和1克B结合能产生2克C，形成C的速率与A和B的数量乘积成正比。同样C也可分解为A和B，C分解的速率正比于C的数量。问题：在给定A和B的数量后，模拟有多少C物质产生出来，以及达到稳定的时间。解：在任何时刻，设a，b，c分别是A，B，C的数量，则它们增加和减少的速度服从下列微分方程：da/dt=k2C-k1a·bdb/dt=k2C-k1a·b(2.8)dc/dt=2k1a·b-2k2C其中k1和k2是比例常数(实际问题中k1和k2会随温度、压力等发生变化，但在模拟过程中为简化模型，可视为常数)，是模型中的参数，-k1a·b中的负号是因为A和B是减少的过程。假设模拟从t=t0(一般取0)开始，使t以△t时间间隔增加(由步长△t和模拟周期可定出所分时间段数)，则相应的欧拉公式为(2.9)式。给出常数k1和k2值以及A、B、C的初始数量，即得迭代公式：an+1=an+(k2Cn-k1an·bn)·△tbn+1=bn+(k2Cn-k1an·bn)·△tcn+1=cn+(2k1an·bn-2k2Cn)·△t(2.9)a(0)=a0,b(0)=b0,c(0)=0k1=k1,0,k2=k2,0由(2.9)式，从t=0开始，由a0、b0、c0可算出a1、b1、c1，又可算出a2、b2、c2，(a2=a(2△\nt)…)，以△t为间隔，进行N=T/△次计算，就可算出周期T的系统状态，从而得出模拟结果。根据数学模型(2.9)就可设计编写模拟程序，可选用任何编程语言，在一些流行的语言中，对常用的模拟算法(比如后面要介绍的龙格-库塔(Runge-kutta)方法)都有相应的用于模拟计算的软件包。开始输入k1,k2,a0,b0,c0,T,△t,Nn=0n=n+1计算并保存an,bn,cnn≤N输出模拟结果结束图2-2例2-1模拟流程图publicclassRate{staticdoublek1,k2,h,t;staticdoubleA[]=newdouble[53];staticdoubleB[]=newdouble[53];staticdoubleC[]=newdouble[53];staticinti;publicstaticvoidmain(Stringargs[]){A[1]=100.0;B[1]=50.0;C[1]=0.0;t=0;h=0.1;k1=0.008;k2=0.002;for(i=1;i<53;i++){System.out.println(i+""+t+","+A[i]+","+B[i]+","+C[i]);strut(i);}}staticvoidstrut(inti){A[i+1]=A[i]+(k2*C[i]-k1*A[i]*B[i])*h;B[i+1]=B[i]+(k2*C[i]-k1*A[i]*B[i])*h;C[i+1]=C[i]+2.0*(k1*A[i]*B[i]-k2*C[i])*h;t=t+h;}}欧拉法是最简单的数值积分法，在介绍更好的数值积分法之前，先介绍几个关于数值积分的基本概念。1、单步法与多步法数值积分法是用递推公式求解，如果仅由前一时刻的数值yn就能算出后一时刻的数值yn+1，则称为单步法，反之，如果求yn+1时需用到yn，yn-1，yn-2，…等多个值，则称为多步法。如，欧拉法是单步法，而改进和欧拉法是多步法。单步法由递推公式自身就能启动运算(由初值即可算出y1，由y1可算出y2，…)，所以它是能自启动的算法；多步法在开始的时候要先用其它的方法计算该时刻前面的函数值，所以不能自启动运算。一般来说，由于多步法更能充分利用多个时刻的信息，所以模拟速度快，精度高，但计算量要大一些，后面还要介绍多步法的算法。2、显式与隐式在递推公式中，计算yn+1时所用到的数据均已算出的计算公式称为显式公式，相反，在算式中隐含未知量yn+1称为隐式公式。如欧拉法中递推公式是显式的，而梯形法中递推公式是隐式的。在隐式公式中，必须先用显式公式估计一个值，再用隐式公式迭代，即预估-校正法。隐式与显式相比，有明显的高精度和稳定性。3、误差在数值解法的过程中，每一步计算都会产生误差，误差的来源有两个方面：一是计算误差(即计算机计算本身的误差)，一个是公式误差(用差分方程代替微分方程)。分别称之为舍入误差和截断误差。\n舍入误差：计算机的字长是有限的，数字不能表示的完全精确，实际上计算的结果是用有限精度的有理数(如计算机中使用的浮点数)来近似无限精度的实数，所以在对动态系统模拟的过程中，舍入误差是不可避免的。通常舍入误差的大小与积分步长成反比，步长h越小，计算次数多则舍入误差大，但不能随意加大步长，否则将产生大的截断误差甚至影响系统稳定性。在给定运算序列的条件下，唯一可减少舍入误差的方法是增加数字表示的精度，如将单精度浮点数(有7位数的精度)表示改为双精度数(有15位数的精度)。实际的动态系统模拟时，多采用双精度数据类型，尽管这样变量占用的存储空间大，处理时间稍长，但对提高模拟的精度来讲是值得的。截断误差：是用差分方程代替微分方程产生的误差，即用数值解代替微分方程的精确解产生的误差，所以是公式本身的误差，一般用台劳级数来分析积分公式的精度：假设在tn时积分(精确值)已经算出，则用台劳级数可求得tn+1时的精确解：y(tn+1)=y(tn+h)=y(tn)+y’(tn)+h2/2!·y”(tn)+…+hr/r!·y(r)(tn)+O(hr+1)(2.10)如果一个数值解法是用前r+1项来近似的计算y(tn+1)，则后面的各项(记为)O(hr+1)是在这一步计算中引进的附加误差，称为这个算法的(局部)截断误差。误差O(hr+1)与hr+1同阶(h→0)，即计算的精度保特了r阶，此时称这个方法是r阶的精度。一个数值方程的阶数可视为衡量这一方法的精确度的重要标志，不同的数值解法有不同的精度。欧拉法只取(2.10)中的前两项作近似计算，所以是一阶精度的算法；梯形法相当于取(2.10)式中的前三项，故是二阶精度的方法。4、稳定性如果系统是稳定的(系统的状态随时间的推移逐步稳定在某个水平上)，则在模拟的迭代过程中，数值积分的解也应该是稳定的，但由于初始数据的误差及在迭代运算中的舍入误差会对后面的计算结果产生影响。当步长h选择不合理时，可能使模拟结果不稳定。对于欧拉法，可用下面的检验方程(其中λ为方程的特征根)：y’=λy(2.11)y(0)=y0(注意其精确解为y=y0eλt)来检验步长对数值解稳定性的影响：对(2.11)表示的方程，欧拉公式为：yn+1=yn+λynh=(1+λh)yn(2.12)y(0)=y0所以，要使数值解稳定，必须使：|1+λh|<1，解得|λh|<2或h<2T(T=1/λ是系统时间常数)。所以要保证欧拉方法计算的稳定性，步长h必须小于系统时间常数的2倍。在(2.11)中特征根λ在一定数量级范围变动，现令λ=1，来作一个具体的模拟，此时h应小于2，否则将不稳定。取h=0.01，h=0.05，h=1.0，h=1.9，h=2.0，h=2.1六个不同的步长，y0=1。模拟结果为：h=0.01：表现出较好精度h=0.05：虽然近似解接近精确解，但存在误差h=1.0：解在一步后趋于零，并一直保持，虽然稳定，但明显精度不高h=1.9：解振荡，幅度值逐衰减，并趋于稳定h=2.0：解不衰减，等幅振荡h=2.1：解的振荡幅度值递增，表明系统不稳定，数值解发散数值解图如图2-3所示。210-1-2210-1-2210-1-2210-1-2210-1-2210-1-2图2-3不同步长的数值解\n2.2龙格-库塔(Runge-Kutta)法R-K方法的基本思想是用台劳展开式的前几项来对微分方程求近似解。再以模型(2.1)为例：y'=f(t,y)(2.1)y(t0)=y0假设从t0开始，以h增长，h1=t0+h，y1=y(t0+h)，在t0附近展开成台劳级数，保留h2项，则有：y1=y0+f(t0,y0)·h+h2/2(δf/δy·dy/dt+δf/δt)(2.13)(此式括号中的导数是在(t0,y0)处的导数值)为了求(2.13)的解，假设(2.13)的解可写成如下的形式：y1=y0+(b1k1+b2k2)·h其中k1=f(t0,y0)(2.14)k2=f(t0+C2h,y0+a1k1h)(注意(2.13)式是关于函数y的导数的算式，而(2.14)式中k1和k2都是y在某点处的导数，相差的只是常量级的系数不同，问题正是要定出这些系数，从而得到数值解表达式，为此)对k2式右端的函数在(t0,y0)处展开台劳级数，保留h项，得：k2≈f(t0,y0)+(C2·δf/δt+a1k1·δf/δy·dy/dt)·h把k1、k2代入(2.14)中，得y1=y0+b1f(t0,y0)h+b2h[f(t0,y0)+(C2δf/δt+a1k1δf/δydy/dt)h](2.15)\n将(2.15)式与(2.14)式比较，得关于系数的方程：b1+b2=1b2C2=1/2(2.16)b2a1=1/2方程组(2.16)中四个未知量，三个方程，故有无穷多解，求出一个解，比如令b1=b2，得一个解：b1=b2=1/2a1=C2=1代入(2.14)，得如下一个公式：y1=y0+1/2(K1+K2)h其中k1=f(t0,y0)k2=f(t0+h,y0+k1h)这是从t0计算t1时刻的公式，写成一般的形式，得如下递推公式：yn+1=yn+h/2(K1+K2)其中k1=f(tn,yn)(2.17)k2=f(tn+h,yn+k1h)模型(2.1)的数值解公式(2.17)即称为R-K公式，这种数值解法即称为R-K方法。在(2.13)式中，只保留了h2项，故公式(2.17)的精度是2阶的，公式(2.17)称为二阶R-K公式。目前在实际模拟问题中，四阶R-K公式用的最为普遍。在推导四阶R-K公式时，在台劳公式中保留到h4项，推导过程与前面类似。一个四阶R-K公式可以是下面的形式：yn+1=yn+h/6(K1+2K2+2K3+K4)其中k1=f(tn,yn)k2=f(tn+h/2,yn+k1h/2)(2.18)k3=f(tn+h/2,yn+k2h/2)k4=f(tn+h,yn+k3h)问题：关于R-K方法还有下面一些问题需了解。(1)多种形式：方程组(2.16)有无穷多解，我们取了一个解，得到了公式(2.17)，也可以取其它的解，所以每一阶R-K公式都有多种形式。二阶R-K公式常用的除(2.17)式外，还有yn+1=yn+K2hk1=f(tn,yn)k2=f(tn+h/2,yn+k1h/2)(对应(2.16)的解b1=0,b2=1,C2=a1=1/2)。四阶R-K公式常用的除(2.18)式外，还有yn+1=yn+h/8(K1+3K2+3K3+K4)k1=f(tn,yn)k2=f(tn+h/3,yn+k1h/3)k3=f(tn+h2/3,yn+k1h/3+k2h)k4=f(tn+h,yn+hk1-hk2+hk3)(2)精度：也可推导3阶、5阶或更高阶的R-K公式，但在一般工程中，四阶公式就完全能达到要求的精度，而且四阶公式是最常用的，所以一般不使用更高阶的公式。另一个特殊情况：一阶R-K公式，只含有h的1次项，即：yn+1=yn+hf(tn,yn)，这就是欧拉公式，所以欧拉公式也可看面一阶R-K公式，其精度最低。R-K方法的精度取决于步长h\n及求解方法。一般来说，为达到同样的精度，四阶方法的步长可以比二阶方法的步长大10倍，而四阶方法每步的计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小。R-K方法可使用较大的步长也是其一个特点。(3)单步法：不论几阶的R-K公式都是单步法，在计算yn+1时只用到yn，可自启动，使用的存储量也小。另外，无论几阶的R-K公式，算式中均有两部分组成：一部分是上一步的结果yn，第二部分是h=tn–tn-1中对f(t,y)的积分，它是步长h乘以各点斜率的加权平均。如在四阶公式(2.18)中，取四点斜率k1、k2、k3、k4，对k2、k3各取两份，而k1和k4各取一份进行加权平均。(4)变步长的R-K方法：在前面讨论的这些数值方法中，都需要在模拟之前选择步长，且假设步长是固定的。这样作有一定的缺点：如果h选的太小，会增加计算量，增长模拟的时间；而太大又会达不到精度要求。另一种方法是在模拟过程中动态调整步长，这样，当状态变量变化缓慢时，步长可取大些，以减少计算时间；当状态变量变化快时，步长再选小些，以保证模拟的精度。这种步长的自动控制是根据每一步的误差大小来实现的，即根据模拟过程中每一步的局部误差大小来调整步长，使误差保持在规定的范围内。为了得到每一步的局部误差，一般采用两种不同阶次的递推公式：一般用低一阶的公式同时计算yn+1，并计算两者之差作为局部误差的估计值。由于附加的计算会使计算量加大，要使计算量最小，可选择R-K系数，使两个公式的ki相同，使中间结果对两个公式都适用。目前使用较多的一个四阶变步长的方法是R-K-Merson(墨森，1957年给出)法：其四阶公式为：yn+1=yn+h/6(K1+4K4+K5)k1=f(tn,yn)k2=f(tn+h/3,yn+k1/3)k3=f(tn+h/3,yn+h/6(k1+k2h))k4=f(tn+h/3,yn+h/8(k1+3k3)k5=f(tn+h,yn+h/2(k1-3k3+4k4)计算估计误差的三阶公式为：ÿn+1=yn+h/6(3K1-9K3+12K4)误差为：E=ÿn+1-yn+1=h/6(2K1-9K3+8K4-K5)在每一步计算后，计算误差，根据误差大小来调整步长，调整策略如下：1)当误差大于预先设定的最大允许误差Emax时，减少步长，一般将步长减半，并以新步长重新计算后再比较。2)当误差小于预先设定的最小误差Emin时，步长增加一倍，以新步长往下计算。3)如步长小于某一下限Hmin则不再减小，以免增加模拟时间，舍入误差过大。4)如果步长大于某一上限Hmax则不再增加步长。这种方法虽增加了计算量，但在多数情况下，假设具有同样的积分误差，变步长的模拟时间要小于固定步长的模拟时间。2.3线性多步法前已讨论过单步法和多步法：单步法中计算yn+1只用到yn和fn的值，而多步法计算yn+1时要用到前面多步的值。单步法计算简单，但多步法更能充分利用前面多步的信息，所以更能加快模拟速度、提高精度。一个线性公式对于输入因子来讲是一次的，所以无论公式推导还是计算，都较简单。在实际模拟工作中对非线性的数学模型一般要作线性化处理，以使求解简便。在线性多步法中以亚当姆斯(Adams)法使用最为普遍。下面介绍这种方法。亚当姆斯公式的一般形式为yn+1=yn+h[B–1f(tn+1,yn+1)+B0f(tn,yn)+…+Bkf(tn-k,yn-k)(2.19)各阶亚当姆斯公式的系数表如下：表2-1Adams系数表\n名称B–1B0B1B2B3一阶显式01000二阶显式03/2-1/200三阶显式023/12-16/125/120四阶显式055/24-59/2437/24-9/24一阶隐式10000二阶隐式1/21/2000三阶隐式5/128/12-1/1200四阶隐式9/2419/24-5/241/240由表2-1知，一阶显式Adams公式yn+1=yn+hf(tn,yn)就是欧拉公式。二阶隐式Adams公式yn+1=yn+h/2[f(tn+1,yn+1)+f(tn,yn)]就是改进的欧拉公式。多步法的缺点是不能自启动，所以，对于显式的Adams公式开始时要用单步法算出需要的值，然后才能用多步法迭代运算。对于隐式的Adams公式，一般先用相应的显式方法计算(预估)，再用隐式法计算(校正)，这种将显式与隐式合起来使用的方法，即为前面介绍的预估-校正法。显式公式计算简单，但隐式公式的稳定性好。实际的模拟工作中四阶Adams隐式公式用的最普遍，其预估-校正公式为：预估ypn+1=yn+h/24[55fn–59fn-1+37fn-2–9fn-3]校正ycn+1=yn+h/24[9fpn+1+19fn-5fn-1+fn-2]同阶的Adams法比R-K方法计算量小，比如四阶：隐式的Adams公式只须计算两次右函数，而R-K法要计算4次右函数，所以Adams法范用于实时模拟系统中。但使用Adams法时，如果f不连续，可能会造成较大的瞬时误差(实际上这种方法是基于插值算法)。本章以上介绍了几种微分方程的数值解法，除了以上介绍外，还有其它方法，对连续系统的模拟，先建立数学模型(其模型是微分方程的初值问题)，然后选择求解微分方程的计算方法。选择数值解法依据以下方面：1)精度要求：在数值解法中存在舍入误差和截断误差，不同的方法误差大小也有区别。一般来说，阶次越高，解就越精确。所以在前述方法中欧拉法精度最低，依次是梯形法、R-K方法、Merson法。此外，为减小误差，可选择较小的步长，但步长越小计算步数越大，计算速度也就越低。2)速度要求：不同的方法计算速度不同，比如同阶的Adams方法比R-K方法要快，在实际模拟中，要求更快的速度。所以计算速度也是衡量程序设计水平的一个重要标志。3)稳定性要求：每一种求解的数值方法都是通过某种离散化手续，将微分方程转化为差分方程(然后以代数方程的形式)来求解，而差分方程的求解会有计算误差，这种误差如果在迭代过程中恶性增长，就会“淹没”差分方程的“真解”，从而造成不稳定。因此，这也是在选择数值解法时应考虑的问题。例2.2\n第三章离散系统的计算机模拟前已讨论离散系统中状态的变化不连续，而是只在一些离散的点上发生。离散系统也分确定性和随机性两种类型：确定性系统在多次测试中对同样的初始条件和控制输入具有相同的响应，而随机系统在相同的输入和初始条件下每次运行行为不同；对于连续系统确定型更为普遍，而对于离散系统而言，确定型属于最优设计(比如前面库存的例子，确定型情况下是定货点、定货量、销售量均为确定)，多数具有随机性；实际上确定性是相对的，随机性是绝对的，自然界中没有系统不受随机因素的干扰(如前面的追击问题，随机的阵风会改变飞机的航向)。故对离散系统，只讨论随机模型，用概率分布来描述。模拟离散系统有两种用时模型：一种是固定时间步长模型，另一种是下一事件模型。前者在程序中产生一个时钟，这个时钟以固定时间步长更新，当一个时间步长结束时，检查系统是否有事件发生，有则处理，否则什么也不作，进入下一个时间步长。缺点：步长定的太小，则可能没有事件发生，浪费机时，使程序效率低下；步长太大，则在一个步长内可能有多个事件发生，而系统一般按一个事件处理，这样就漏掉了其它事件的差别，使模拟不准确；后者时间间隔的长短由两个事件出现的间隔而定，算法：系统中所有考虑的事件放在“事件表”中，由系统时钟选取一个最早发生的事件，并处理这个事件，然后由该事件发生的时间修改时钟。即系统不断地从事件表中选择事件、处理事件、修改时间，使模拟一步一步地向前推进。以后者用的较多。本章主要通过排队系统介绍离散系统模拟的基本方法。3.1排队系统排队系统是是常生活中经常遇到的现象，比如去食堂买饭要排队、去理发店理发要排队、去医院看病要排队等等。一般来说，当要求服务的数量超过服务机构(服务台、服务员)的容量，即到达的顾客不能立即得到服务，就会出现排队现象。排队系统是最典型、最重要的离散系统，很多问题都可归结为排队系统，从而用排队系统来模拟：如通信系统中电话占线问题、交通系统中车辆的堵塞和疏导、故障机器的停机待修、水库的存储调节、生产线的产品加工等等。在这些问题中顾客(不一定是人，要广范地理解)到达的时间和服务时间都是随机的，所以排队(队列可以是具体的排列，也可以无形的)不可避免。排队系统也称为随机服务理论。在排队系统中，要研究的问题是服务台的工作效率(或设备利用率)：比如在食堂买饭的排队系统中，应设几个窗口最为合理，如果增加服务窗口，就要增加投资或发生空闲浪费；如果卖饭窗口太少，排队现象就会严重，对顾客和社会都会带来不得的影响；在一个理发店中设置几个服务员才能即工作效率高，又不损失顾客等。因此，必须考虑在两者之间取得平衡。3.1.1排队系统的组成和特征一个排队系统都有三个基本组成部分：1)输入过程(即顾客到达的方式)2)排队规则3)服务机构不同的排队系统有不同的特征。1输入过程：有以下几种情况1)顾客的总体(称为顾客源)可能是有限的，也可能是无限的。比如上游河水流入水库可认为是无限总体，而一个工厂内停机待修的机器就是有限总体。2)顾客到达的方式可以是单个的，出可以是成批的。一般只讨论单个到达的情况。3)顾客相继到达的间隔时间可以是确定型的，也可以是随机型的。比如：在自动装配线上装配的机器部件必须按确定时间间隔到达装配点，但更多是随机型的。对于随机型，要知道单位时间内顾客到达数的概率分布。\n关于概率分布：前已说明，对于随机系统其模型不能用方程来描述，而是用概率分布描述。在排队系统中，首先要知道顾客到达的概率分布：如果不知道系统的理论分布，就必须先根据统计资料作出经验分布，然后按照统计学的方法(比如用假设检验的方法)确定服从哪种理论分布，并估计其中的参数。比如，在一个排队系统中统计一个月(30天)顾客到达的情况，每天有多少个顾客到达；然后根据源始数据统计出30天中有0个顾客到达的天数、1个顾客到达的天数、2个顾客到达的天数，…；根据这个统计资料可算出一天有k(k=1,2,…)个顾客到达的频率；用频率近似代替概率就得到的概率分布；然后将这种分布跟某种已知的理论分布作χ2检验；最后就可得到系统中顾客到达的概率分布。在排队系统中，在某些假设条件下，一个时间段中到达的顾客数服从泊松(Poisson)分布：记N(t)为在时间区间[0，t)内到达的顾客数，Pn(t)为在[0，t)内到达n个的概率。Pn(t,t+△t)表示在[t,t+△t)内有n个顾客到达的频率，即：Pn(t)=P{N(t)=n}Pn(t,t+△t)=P{N(t+△t)–N(t)=n}(△t>0,n≥0)若Pn(t,t+△t)满足下面三个条件，则说顾客的到达形成泊松流：1)在不相叠的时间区间内顾客的到达数相互独立(对无限总体是相互独立的，但有限总体未必，如一个工厂的机器在一个短时间内出现停机(顾客到达)的概率受已经待修数目的影响)。2)对充分小的△t，在[t,t+△t)内有一个顾客到达的频率与t无关，且均与时间长△t成正比，即P1(t,t+△t)=λ△t+o(△t)其中λ是常数，表示单位时间内1个顾客到达的概率；o(△t)是关于△t(当△t→0时)的高阶无穷小3)对充分小的△t，在内[t,t+△t)内有2个或2个以上到达的概率极小，可略。即：Pn(t,t+△t)=o(△t)可以证明，在以上假设下，内[0,t)内的顾客到达数N(t)服从泊松分布：Pn(t)=(λt)n/n!•e-λtt>0(3.1)N=1,2,…也可以证明：当N(t)服从P(λ)时，顾客相继到达的间隔时间服从指数分布，即密度函数为：f(t)=λe-λtt≥0其均值1/λ是平均到达的间隔时间，故在确定系统的概率分布时出可通过统计平均间隔时间来估计参数。2排队规则：有以下几种情况1)顾客到达时，如果所有服务台忙，则顾客可以离去，也可以等待，前者叫即时制，后者叫等待制。普通的电话呼叫属于前者，而登记市外长途电话的呼叫属于后者。对于等待制，为顾客服务的次序可有下面规则：先到先服务(FIFO)：即接次序接受服务后到先服务(LIFO)：如仓库中存放的没时间限制的物品；情报系统中最后到达的往往是最有价值的。随机服务有优先权的服务：如医院对病情严重的患者将优先治疗。2)从占用的空间来看，队列可以排在具体的处所；也可以抽象的。由于空间的限制或其它原因，有的系统要规定容量(队长)的最大限；有的是无限制的。3)从队列的数目看，有单列、多列(此时各列的顾客有的可以转移，有的不能)，前者较简单。3服务机构：从服务机构和工作模式上看有以下几种情况1)服务方式可以对单个顾客，也可以成批服务。如车站等候的顾客是成批服务。\n2)对顾客的服务时间可以是确定型，也可以是随机的。如果输入和服务时间都是确定型，就太简单了，这属于最优设计。多数是随机性的，对随机性的必须知道概率分布，概率分布可由以往的经验或专业知识得到。以上从三个方面介绍了不同的情况，每种情况的组合便构成一个排队模型，不同的排队模型有不同的模拟方法。在实际模拟一个排队问题时，首选要确定属于哪种类型，其中只有顾客到达的分布和服务时间的分布实测数据确定，其它因素在系统中均为给定。3.1.2排队系统的指标模拟一个排队系统的目的，是研究排队系统运行的效率，估计服务质量，确定系统参数的最优值。以决定系统的结构是否合理、研究设计改进措施等。这就必须确定用于衡量一个系统好坏的数量指标。不同的排队模型衡量标准及指标的算法不尽相同，下面讨论最简单的情况：单服务台、单队列、排队规则为先进先出、顾客到达的模式为泊松分布、服务时间为指数分布。对这一排队模型，常用下述指标衡量。1稳态平均延迟时间DD=Di/n其中Di为第i个顾客的延迟时间(即排队等待的时间)。故Di/n是n个顾客平均等待的时间，它是一个随机变量，在具体的模型中，通过建立差分方程，可求出其概率(差分方程的解)，称为瞬态解。一般来说，求瞬态解不易，即使求出也很难利用，因此，常用它的极限，称为稳态。稳态的物理含意是：当系统运行了无限长时间之后，初始(t=0)出发状态的概率分布的影响将消失，而系统的状态概率分布不再随时间变化。在实际应用的多数问题中，系统会很快趋于稳态，不需t→∞。2稳态平均滞留时间WW=Wi/n=(Di+Si)/n其中Wi为第i个顾客通过系统的滞留时间，它等于该顾客排队等待的时间Di和接收服务时间Si之和。在一些问题中，如机器故障中，无论等待修理或正在修理，都使工厂受到损失；而在购物、就诊等问题中顾客关心的是等待时间。3稳态平均队长QQ=Q(t)dt/T其中Q(t)为t时刻的队长，T为系统模拟时间。一般来说，Q越大说明服务率越低。排队成龙，是顾客最厌烦的。4系统中的态平均顾客数LL=L(t)dt/T=(Q(t)+S(t))dt/T其中L(t)为t时刻系统中的顾客数，它等于队列中的顾客数Q(t)与正在服务的顾客S(t)数之和。5服务台利用率ρρ==使用ρ可计算出以下指标：空闲率：1-ρ、Q=ρ2/(1-ρ)、L=ρ/(1-ρ)、W=L/、D=Q/\n以上是最简单的情况，在实际的排队系统中要复杂的多。下面就这种最简单的情况，举例通过手算来看一下模拟的过程。例3.1设：ti：第i个顾客到达的时间Ai：第i个顾客与前一顾客到达的时间间隔，即ti-ti-1Di：第i个顾客排队等待的时间Si：第i个顾客服务时间Ci：第i个顾客离去的时间，即ti+Di+Sibi：事件表中第i个事件发生的时间(此问题中的事件包括：顾客到达，顾客离去2种事件，bi是第i个任何一类事件发生的时间)，定义系统事件类型：I类：顾客到达，II类：顾客离去qi：第i个事件发生时的队长Zi：第i个事件发生时服务台状态，1为忙，0为闲些系统中，Ai和Si是随机变量，应先确定其分布(比如，根据专业知识或统计的方法确定分布类型，再由统计数据估计参数)，在这里给出一些样本值，根据这些值进行模拟：A1=15、A2=29、A3=24、A4=40、A5=22(分钟)S1=38、S2=36、S3=32、S4=23一般来说，模拟过程如下：(1)执行初始化操作：1)置初始时间t=t0，结束时间一般取t∞(此例中模拟150个时间单位)2)对事件表初始化，置系统初始事件的状态(2)推进模拟时钟：Time=t//修改时间While(TimeC1=53)b3=C1=53(<150，故处理：)Z3=0，D2=C1-t2=9，q3=q2-1=0，S2=36，故C2=C1+S2=53+36=89取下一事件：第3个顾客到达(t3=68