数学建模计算机模拟的技术实用教案

1、1 模拟基础模拟基础1.0 模拟的背景模拟的背景(bijng)、思路、思路应用领域：应用领域：第二次世界大战期间，第二次世界大战期间，J.V.Neumann等人将进行等人将进行的的“中子扩散中子扩散”的科研项目取名为的科研项目取名为“Monte-Carlo”运输系统模拟运输系统模拟摩天大楼安全疏散系统模拟摩天大楼安全疏散系统模拟国民经济发展模拟国民经济发展模拟人口增长系统模拟人口增长系统模拟供水系统模拟供水系统模拟第1页/共82页第一页，共82页。管理系统模拟管理系统模拟雷达系统模拟雷达系统模拟战争系统模拟战争系统模拟模拟思路：模拟思路：“模拟模拟”对系统抽象建模对系统抽象建模“试验试验”根据

2、模型设计算法，编程进行反复试验根据模型设计算法，编程进行反复试验“估计估计”根据试验数据根据试验数据“收集收集(shuj)”根据试验结果作出判断根据试验结果作出判断第2页/共82页第二页，共82页。1.1 1.1 模拟(mn)(mn)的基本知识 1.1.1 模拟的概念及作用 现实系统的数学或逻辑模型可能十分复杂，例如大多数具有随机因素的复杂系统，其中的一些随机性因素很难用准确的数学公式表述，从而也无法对整个系统采用解析法求解。模拟是处理这类实际问题的有力(yul)工具。第3页/共82页第三页，共82页。 模拟通常借助于计算机进行。计算机模拟：在已经建立的数学、逻辑模型的基础之上，通过(tngg

3、u)(tnggu)计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。第4页/共82页第四页，共82页。模拟的作用：对于很难用解析方法加以处理的问题，模拟是一种有效的技术；对建模过程中的假设进行鉴定，对理论研究的结论加以检验；对不同的实现方案进行多次模拟，按照既定(j dn)的目标函数对不同方案进行比较，从中选择最优方案。第5页/共82页第五页，共82页。 1.1.2 1.1.2 模拟的分类 通常，模拟时间是模拟的主要自变量。 设计正确的模拟时间推进机理：模拟过程中应根据系统的特性正确推进模拟时间，使系统中各要素与发生的事件保持(boch)(bo

4、ch)同步。第6页/共82页第六页，共82页。 推进模拟时间的基本方法： 下次事件法：将模拟时间由一个事件发生的时间点推进到紧接着的下一次事件发生的时间点。 固定时间步长法：模拟时间每次均以相等的固定步长向前推进，每到达一个新的模拟时间点需检查相应时间段内是否发生了事件。需根据实际(shj)问题合理设置模拟时间发生改变的步长。第7页/共82页第七页，共82页。 根据模拟过程(guchng)中因变量的变化情况进行分类： 1）离散型模拟：因变量在与事件时间有关的具体模拟时间点呈离散性变化。大多数系统（如排队服务系统）可采用离散型模拟。 时间推进方法：一般采用下次事件法 应当重点对系统状态可能发生改

5、变的事件进行描述，并确定这些事件之间的逻辑关系。第8页/共82页第八页，共82页。第9页/共82页第九页，共82页。 排队系统通常(tngchng)采用离散型模拟模型。其中，发生系统状态变化的事件有两个：一是有顾客到达；二是服务员完成服务。将最近发生上述两种事件之一的时刻设置为下次事件发生点，就可将服务过程描述为图2所示的模拟模型。第10页/共82页第十页，共82页。 2）连续型模拟：因变量随时间的改变呈连续性变化。在大多数计算机模拟过程中，按固定的步长推进模拟时间。 通常需建立一系列的由系统(xtng)状态变量组成的状态方程组，以描述状态变量与模拟时间的关系。 3）混合型模拟：因变量随时间的

6、推移而作连续性的变化并具有离散性的突变，如库存控制系统(xtng)。第11页/共82页第十一页，共82页。 1.1.3 模拟的方式 终态模拟：在规定的时间T内进行模拟运行，时间达到T时，模拟终止。其性能指标明显取决于系统的初始状态。 稳态模拟：随着模拟时间的推移，系统的性能逐渐趋于平稳。其目的是研究非终态系统长期运行条件下的稳态性能，模拟时间的长短取决于能否获得系统性能的优良(yuling)估计（可由模拟输出的精度确定）。第12页/共82页第十二页，共82页。 1.1.4 模拟的一般步骤 明确问题，建立模型。正确描述待研究问题，明确规定模拟的目的和任务，确定衡量系统性能或模拟输出结果的目标函数

7、，然后根据(gnj)系统的结构及作业规则，分析系统各状态变量之间的关系，以此为基础建立所研究的系统模型。第13页/共82页第十三页，共82页。 收集和整理数据资料。模拟技术的正确运用，往往由大量的输入数据作依靠。在随机模拟中，应认真分析具体收集到的随机性数据资料，确定系统中随机性因素的概率分布特性，以此为依据产生模拟过程所必需的抽样(chu yn)数据。 编制程序，模拟运行。 分析模拟输出结果：模拟结果的统计特性（样本均值、方差、置信区间等），灵敏性分析，选择最优方案。第14页/共82页第十四页，共82页。 注：模拟结果的统计分析模拟的输出结果是分布特征未知的随机变量，每次运行的结果仅仅是对该

8、随机变量所有(suyu)观察值总体的一次抽样，对总体的代表性很差，虽然可以增加模拟运行的时间从而增加抽样次数，但这些数据总是由一个“种子”经过一定的算法而获得的伪随机序列，它们是自相关的，并不能构成统计上独立的随机样本。第15页/共82页第十五页，共82页。 【问题】如图，一列火车(huch)(huch)从A A站开往B B站，某人每天赶往B B站上这趟火车(huch).(huch).AB火车火车(huch)运行运行方向方向某人某人(mu rn)他已了解到：他已了解到： 1）火车从火车从A站到站到B站的运行时间是均值为站的运行时间是均值为30分分钟，标准差为钟，标准差为2分钟的随机变量；分钟的

9、随机变量； 思考：思考：请研究他能否赶上这趟火车。请研究他能否赶上这趟火车。1.2 随机模拟案例：赶上火车的概率第16页/共82页第十六页，共82页。出发时刻出发时刻午后午后1:00 午后午后1:05 午后午后1:10频频 率率 0.7 0.2 0.1他到达他到达(dod)B (dod)B 站的时刻的频率站的时刻的频率分布为分布为时刻时刻 午后午后1:28 午后午后1:30 午后午后1:32 午后午后1:34频率频率 0.3 0.4 0.2 0.1他能否及时赶上他能否及时赶上(n shn)火车？火车？ 2）火车在下午大约）火车在下午大约1点离开点离开A站，离开时刻站，离开时刻的频率的频率(pn

10、l)分布如下：分布如下：第17页/共82页第十七页，共82页。明确明确(mngqu)问题：他能及时赶上火车的概率是多少？问题：他能及时赶上火车的概率是多少？i）分析法：用概率）分析法：用概率(gil)统计知识建立分析模型，统计知识建立分析模型，求解析解。（思考）求解析解。（思考）ii）模拟法：用概率统计知识建立模型，通过模拟求近）模拟法：用概率统计知识建立模型，通过模拟求近似解。似解。 即先建立模拟模型，然后通过计算机模拟得到问题的即先建立模拟模型，然后通过计算机模拟得到问题的近似解。在同样条件近似解。在同样条件(tiojin)下多次试验，计算他下多次试验，计算他能及时赶上火车的频率。能及时赶

11、上火车的频率。建模方向（思路）：建模方向（思路）：第18页/共82页第十八页，共82页。问题问题(wnt)分析：能及时赶上火车的充要条件是：分析：能及时赶上火车的充要条件是：其中其中 T1T1火车火车(huch)(huch)从从A A站出站出发的时刻；发的时刻；T2T2火车火车(huch)(huch)的的运行时间；运行时间；T T3 3他到达他到达B B站的时刻。站的时刻。是什么变量？是什么变量？如何模拟？如何模拟？213213TTTTTT即第19页/共82页第十九页，共82页。 ii）将午后）将午后1时记为时记为t=0，设火车，设火车(huch)运行时间运行时间T2服从正态分布：服从正态分布

12、： T2 N（30，22）。）。 基本假设：基本假设：i）假设）假设T1，T2，T3都是相互独立都是相互独立(dl)的随的随机变量机变量; 第20页/共82页第二十页，共82页。火车出发时刻火车出发时刻(shk)T1和人到达和人到达B站时刻站时刻(shk)T3的分布律分别为：的分布律分别为： T1（分） 0 5 10P(t) 0.7 0.2 0.1T3（分） 28 30 32 34P(t) 0.3 0.4 0.2 0.1建立模型：为了建立模型：为了(wi le)简化计算，将下午简化计算，将下午1点记为初点记为初始时刻。得到随机变量始时刻。得到随机变量T1和和T3的分布律如下：的分布律如下：能及

13、时能及时(jsh)赶上火车的概率赶上火车的概率 p=PT3T1T2 第21页/共82页第二十一页，共82页。如果如果r为在为在(0,1)均匀分布的随机数，为了模拟随均匀分布的随机数，为了模拟随机变量机变量(su j bin lin)T1和和T3，可以通过如下，可以通过如下方法：方法：. 19 . 0,10, 9 . 07 . 0, 5, 7 . 00, 01rrrt. 0 . 19 . 0,34, 9 . 07 . 0,32, 7 . 03 . 030, 3 . 00,283rrrrt，则t1和t3可分别(fnbi)用来模拟随机变量T1和T3。第22页/共82页第二十二页，共82页。主要变量主

14、要变量(binling)说明：说明：n 模拟次数模拟次数k 临时变量临时变量(binling)，存储当前累计模拟次数，存储当前累计模拟次数count 存储赶上火车的次数存储赶上火车的次数两种不同风格两种不同风格(fngg)的算法的算法描述描述模拟算法(sun f)设计输入：输入：赶火车次数（天数）赶火车次数（天数）输出：输出：赶上火车的频率赶上火车的频率第23页/共82页第二十三页，共82页。第第 1 1 步步 输入模拟次数输入模拟次数n n第第 2 2 步步 k=1 k=1，count=0count=0第第 3 3 步步 当当k=nk T3, count=count+1 IF T1+T2 T

15、3, count=count+1，ENDEND第第 10 10 步步 执行第执行第3 3步步第第 11 11 步步 输出输出(shch)(shch)赶上火车频率赶上火车频率p=count/np=count/n算法算法1 1（分步骤（分步骤(bzhu)(bzhu)描述）描述）第24页/共82页第二十四页，共82页。i i）初始化：）初始化： 输入模拟次数输入模拟次数n;n; count=0; count=0;ii)ii)模拟模拟n n次次for i=1 for i=1 toto n, n,模拟随机变量模拟随机变量(su j bin (su j bin linlin)T1,T2,T3)T1,T2,

16、T3， 分别赋给分别赋给t1,t2,t3;t1,t2,t3;if t1+t2 t3,if t1+t2 t3,count=count+1count=count+1end if end if end forend forapp_prb=count / n;app_prb=count / n;算法算法2 2（伪代码（伪代码(di m)(di m)描述）描述）第25页/共82页第二十五页，共82页。n=input(输入模拟次数：输入模拟次数：);count=0;for i=1:n, rt1=rand; %模拟随机变量模拟随机变量(su j bin lin)t1（火车（火车从从A站出发的时刻）站出发的时

17、刻） if rt1=0.7 & rt10.9 T1=5; else T1=10; end T2=30+randn*2; %模拟随机变量模拟随机变量(su j bin lin)t2（火车的运行时间）（火车的运行时间） %模拟随机变量模拟随机变量(su j bin lin)t3（他到达（他到达B站的时站的时刻）刻） rt3=rand; if rt3=0.3 & rt3=0.7 & rt30.9 T3=32; else T3=34; end if T3 T1 + T2, count=count+1; endend%forprob=count/n 模拟程序第26页/共82页第二十六页，共82页。序号模

18、拟次数近似概率p110000.6280210000.6920310000.6530450000.6490550000.6260650000.6288模拟(mn)结果：每次模拟(mn)1000次或5000次第27页/共82页第二十七页，共82页。系统模拟注意事项：系统模拟注意事项：一次模拟结果毫无意义！一次模拟结果毫无意义！模拟是试验性的模拟是试验性的,是思维结果的验证。是思维结果的验证。必须必须(bx)进行足够多次的模拟，并对结果进行统进行足够多次的模拟，并对结果进行统计分析。计分析。 系统模拟特点：系统模拟是研究系统，特别是动态系统模拟特点：系统模拟是研究系统，特别是动态系统的重要方法，对于

19、：系统的重要方法，对于：结构复杂的系统；结构复杂的系统； 很难用解析方法求出变量关系的系统；很难用解析方法求出变量关系的系统；内部机理不明的内部机理不明的“黑箱黑箱”系统；系统； 为验证用其他方法建立的模型及结果为验证用其他方法建立的模型及结果,应是较好的应是较好的选择。选择。第28页/共82页第二十八页，共82页。1.3 随机变量的建模随机变量的建模利用理论分布，基于对问题的实际的、合理的假设，利用理论分布，基于对问题的实际的、合理的假设，选择适当的理论分布模拟随机变量选择适当的理论分布模拟随机变量优点是给出了各种理论结果出现的概率，便于进行数优点是给出了各种理论结果出现的概率，便于进行数学

20、分析和处理。但此方法仅限于十分简单的情况，当学分析和处理。但此方法仅限于十分简单的情况，当问题越复杂，数学处理变得越困难，并且丢失了试验问题越复杂，数学处理变得越困难，并且丢失了试验数据数据(shj)的信息。的信息。基于实际数据基于实际数据(shj)的频率做近似模拟的频率做近似模拟优点在于完全与观察数据优点在于完全与观察数据(shj)相符，并且随实际问相符，并且随实际问题的复杂程度增大不会产生更大的困难，仅增大工作题的复杂程度增大不会产生更大的困难，仅增大工作量而已。缺点是不便于进行数学分析，不得不依赖于量而已。缺点是不便于进行数学分析，不得不依赖于模拟得到的统计结果。模拟得到的统计结果。应用

21、常将两种模拟方法结合起来使用应用常将两种模拟方法结合起来使用第29页/共82页第二十九页，共82页。1.4 1.4 均匀分布随机变量均匀分布随机变量(su j bin lin)(su j bin lin)模模拟拟1.4.1 1.4.1 平方取中法平方取中法第30页/共82页第三十页，共82页。第31页/共82页第三十一页，共82页。1.4.2 线性同余法第32页/共82页第三十二页，共82页。线性同余法程序(chngx)：产生0c整数第33页/共82页第三十三页，共82页。第34页/共82页第三十四页，共82页。1.5 Matlab随机模拟函数(hnsh)：1.5.1 常见分布随机变量的模拟第

22、35页/共82页第三十五页，共82页。第36页/共82页第三十六页，共82页。第37页/共82页第三十七页，共82页。第38页/共82页第三十八页，共82页。第39页/共82页第三十九页，共82页。第40页/共82页第四十页，共82页。第41页/共82页第四十一页，共82页。第42页/共82页第四十二页，共82页。1.5.2 其它(qt)随机变量的模拟第43页/共82页第四十三页，共82页。第44页/共82页第四十四页，共82页。第45页/共82页第四十五页，共82页。2 Monte-Carlo2 Monte-Carlo模拟(mn)(mn)2.1 Monte-Carlo2.1 Monte-Ca

23、rlo原理原理(yunl)(yunl)第46页/共82页第四十六页，共82页。2.2 2.2 蒙特卡罗法应用(yngyng)(yngyng)2.2.1 求解求解(qi ji)非线性规划非线性规划第47页/共82页第四十七页，共82页。第48页/共82页第四十八页，共82页。第49页/共82页第四十九页，共82页。第50页/共82页第五十页，共82页。2.2.2 2.2.2 估算( sun)( sun)圆周率第51页/共82页第五十一页，共82页。第52页/共82页第五十二页，共82页。第53页/共82页第五十三页，共82页。第54页/共82页第五十四页，共82页。2.2.3 2.2.3 估算估

24、算( sun)( sun)定积分定积分第55页/共82页第五十五页，共82页。1 1）频率）频率(pnl)(pnl)法法第56页/共82页第五十六页，共82页。第57页/共82页第五十七页，共82页。第58页/共82页第五十八页，共82页。第59页/共82页第五十九页，共82页。2 2）平均值法）平均值法第60页/共82页第六十页，共82页。第61页/共82页第六十一页，共82页。第62页/共82页第六十二页，共82页。第63页/共82页第六十三页，共82页。第64页/共82页第六十四页，共82页。 例1：某港口提供有足够的泊位供船舶停靠，但现在仅有一个可供装卸的泊位，船舶先到则先进行装卸，如

25、果船舶得不到及时装卸而造成的滞期费为每小时100元。现要弄清该系统的性能，重点考察船舶进入该港后等待装卸的滞留时间(shjin)以及卸位（即装卸用的泊位）的利用率，从而进行经济效益分析。3 模拟模拟(mn)模型案例模型案例第65页/共82页第六十五页，共82页。到达间隔(小时)56 7 89101112131415161718到达船舶数(艘)13 6 791011111197654首先，对进入该港口的首先，对进入该港口的100艘船进行实际调查，记录其活动情况，得到这艘船进行实际调查，记录其活动情况，得到这100艘艘船到达港口的时间间隔船到达港口的时间间隔(jin g)和装卸时间的分布情况的频数

26、和累积频率分布。和装卸时间的分布情况的频数和累积频率分布。装卸时间（小时）910111213141516装卸船舶数（艘）2022191610832表1 100艘船到达(dod)港口的时间间隔频数表表2 100艘船装卸时间(shjin)频数表第66页/共82页第六十六页，共82页。到达间隔56789101112131415161718累积频率0.010.040.100.170.260.360.470.580.690.780.850.910.961.00装卸时间910111213141516累积频率0.200.420.610.770.870.950.981.00表3 船舶到达港口的时间(shjin

27、)间隔累积分布表表4 船舶(chunb)装卸时间累积分布表第67页/共82页第六十七页，共82页。为了比较准确地反映系统的性能，我们采用稳态模拟方式，可以考察系统运行360天的情况。假定港口每天24小时连续工作，因此模拟时间(shjin)T设置为8640小时。考虑到船舶到达港口的事件是影响整个系统状态的主要因素，可以将模拟事件设置为该事件的发生时刻。第68页/共82页第六十八页，共82页。 所用变量说明： t：一艘船到达港口的时间； tl：前一艘船驶离港口的时间； td：两艘船到达港口的时间间隔； ts：当前船舶装卸所需时间； tw：t时刻所有(suyu)已到达港口的船舶等待装卸时间总和； t

28、f：t时刻装卸位已空闲时间总和；第69页/共82页第六十九页，共82页。初始化数据(shj) t=tl=tf=tw=0产生船舶(chunb)到达间隔时间 td设置模拟时间 t=t+tdt8640打印输出产生装卸(zhungxi)服务时间tsttl统计装卸位空闲时间tf=tf+t-tl统计船舶等候时间tw=tw+tl-t设置船舶离港时间Tl=t+ts设置船舶离港时间Tl=tl+ts是是否否第70页/共82页第七十页，共82页。模拟次数到达船舶数到港船舶滞留总时间(小时)船舶平均滞留时间(小时)港口支付滞留费(元)装卸泊位空闲时间(小时)装卸泊位利用率1736748010.167480003539

29、6%272243696.0543690051594%3737910912.3690190042695%473657587.8257580038196%5741750110.1275010034196%673968409.2668400042295%表5 某港口对船舶服务(fw)的模拟结果第71页/共82页第七十一页，共82页。经过6次模拟计算，装卸泊位的平均利用率为95.3%，但到达的船舶的却因得不到及时的服务而造成平均每艘滞留9.295小时，每年到港船舶总滞留6843小时，因此港口每年约需支付68万元的滞留费，这显然是一笔不小的开支。为此，可修改模拟模型，考虑增设一个(y )装卸泊位，并重新

30、进行模拟运行，将这些数据提供给决策者一确定是否需要再新建一个(y )装卸泊位。第72页/共82页第七十二页，共82页。 练习：某公共汽车(qch)站每隔30分钟到达一辆汽车(qch)，但可能有0,3分钟的误差，此误差大小与前一辆汽车(qch)的运行无关。汽车(qch)最多容纳50名旅客，到达该汽车(qch)站时车内旅客人数服从20，50的均匀分布。到站下车的旅客人数服从3，7的均匀分布，每名旅客下车的时间服从1，7秒的均匀分布。旅客按每30分钟到达12个人的泊松分布到达汽车(qch)站，单列排队等车，先到先上，如果某位旅客未能上车，他将不再等候。旅客上车时间服从第73页/共82页第七十三页，共

31、82页。 4，12秒的均匀分布。上下车的规则是：先下后上，逐个上车，逐个下车。 假设每天共发25辆汽车，现在(xinzi)要求模拟一天汽车的运行情况，了解一天中在站内等候汽车的总人数、能上车及不能上车的人数、旅客排队时间分布情况、不能上车人数的分布情况等。第74页/共82页第七十四页，共82页。 例2：露天矿用电铲采掘矿石，然后用卡车运送到卸场。假设共有(n yu)m辆卡车装运，有n台（nm）电铲同时采掘供n辆卡车同时装卸；卸场有s个卸位（s=m），可供s辆卡车同时卸料。装运过程以班为单位，每个班开始工作时，m辆中的n辆卡车由n台电铲装车，其余mn辆卡车排队等候，每辆卡车装满矿石后驶离装料处，排于队列首位的待装卡车驶至空闲电铲前并调转车头接收装料，而刚才已装满的卡车驶至卸场时，若卸位空闲则调头卸完矿石并重新驶会采掘场，排在待装卡车的队尾。第75页/共82页第七十五页，共82页。为了充分发挥每台电铲、每辆卡车的效率，提高班产量(chnling)，就需要确定电铲、卡车及卸位之间的一个适当匹配数量。假定我们总共要模拟20个班运行的情况。已知每个班连续工作6小时，并且每个班开始工作时的初始状态均一致，