数学建模_计算机模拟技术1

模拟技术河南师范大学数学学院1模拟基础2Monte-Carlo模拟3模拟模型案例1模拟基础1.0模拟的背景、思路应用领域：•第二次世界大战期间，J.V.Neumann等人将进行的“中子扩散”的科研项目取名为“Monte-Carlo”•运输系统模拟•摩天大楼安全疏散系统模拟•国民经济发展模拟•人口增长系统模拟•供水系统模拟•管理系统模拟•雷达系统模拟•战争系统模拟模拟思路：•“模拟”－对系统抽象建模•“试验”－根据模型设计算法，编程进行反复试验•“估计”－根据试验数据•“收集”－根据试验结果作出判断1.1模拟的基本知识1.1.1模拟的概念及作用•现实系统的数学或逻辑模型可能十分复杂，例如大多数具有随机因素的复杂系统，其中的一些随机性因素很难用准确的数学公式表述，从而也无法对整个系统采用解析法求解。模拟是处理这类实际问题的有力工具。•模拟通常借助于计算机进行。计算机模拟：在已经建立的数学、逻辑模型的基础之上，通过计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。模拟的作用：•对于很难用解析方法加以处理的问题，模拟是一种有效的技术；•对建模过程中的假设进行鉴定，对理论研究的结论加以检验；•对不同的实现方案进行多次模拟，按照既定的目标函数对不同方案进行比较，从中选择最优方案。1.1.2模拟的分类通常，模拟时间是模拟的主要自变量。设计正确的模拟时间推进机理：模拟过程中应根据系统的特性正确推进模拟时间，使系统中各要素与发生的事件保持同步。推进模拟时间的基本方法：•下次事件法：将模拟时间由一个事件发生的时间点推进到紧接着的下一次事件发生的时间点。•固定时间步长法：模拟时间每次均以相等的固定步长向前推进，每到达一个新的模拟时间点需检查相应时间段内是否发生了事件。需根据实际问题合理设置模拟时间发生改变的步长。根据模拟过程中因变量的变化情况进行分类：1）离散型模拟：因变量在与事件时间有关的具体模拟时间点呈离散性变化。大多数系统（如排队服务系统）可采用离散型模拟。时间推进方法：一般采用下次事件法应当重点对系统状态可能发生改变的事件进行描述，并确定这些事件之间的逻辑关系。排队系统通常采用离散型模拟模型。其中，发生系统状态变化的事件有两个：一是有顾客到达；二是服务员完成服务。将最近发生上述两种事件之一的时刻设置为下次事件发生点。2）连续型模拟：因变量随时间的改变呈连续性变化。在大多数计算机模拟过程中，按固定的步长推进模拟时间。通常需建立一系列的由系统状态变量组成的状态方程组，以描述状态变量与模拟时间的关系。3）混合型模拟：因变量随时间的推移而作连续性的变化并具有离散性的突变，如库存控制系统。1.1.3模拟的方式终态模拟：在规定的时间T内进行模拟运行，时间达到T时，模拟终止。其性能指标明显取决于系统的初始状态。稳态模拟：随着模拟时间的推移，系统的性能逐渐趋于平稳。其目的是研究非终态系统长期运行条件下的稳态性能，模拟时间的长短取决于能否获得系统性能的优良估计（可由模拟输出的精度确定）。1.1.4模拟的一般步骤•明确问题，建立模型。正确描述待研究问题，明确规定模拟的目的和任务，确定衡量系统性能或模拟输出结果的目标函数，然后根据系统的结构及作业规则，分析系统各状态变量之间的关系，以此为基础建立所研究的系统模型。•收集和整理数据资料。模拟技术的正确运用，往往由大量的输入数据作依靠。在随机模拟中，应认真分析具体收集到的随机性数据资料，确定系统中随机性因素的概率分布特性，以此为依据产生模拟过程所必需的抽样数据。•编制程序，模拟运行。•分析模拟输出结果：模拟结果的统计特性（样本均值、方差、置信区间等），灵敏性分析，选择最优方案。注：模拟结果的统计分析模拟的输出结果是分布特征未知的随机变量，每次运行的结果仅仅是对该随机变量所有观察值总体的一次抽样，对总体的代表性很差，虽然可以增加模拟运行的时间从而增加抽样次数，但这些数据总是由一个“种子”经过一定的算法而获得的伪随机序列，它们是自相关的，并不能构成统计上独立的随机样本。【问题】如图，一列火车从A站开往B站，某人每天赶往B站上这趟火车.AB火车运行方向某人他已了解到：1）火车从A站到B站的运行时间是均值为30分钟，标准差为2分钟的随机变量；思考：请研究他能否赶上这趟火车。1.2随机模拟案例：赶上火车的概率出发时刻午后1:00午后1:05午后1:10频率0.70.20.1他到达B站的时刻的频率分布为时刻午后1:28午后1:30午后1:32午后1:34频率0.30.40.20.1他能否及时赶上火车？2）火车在下午大约1点离开A站，离开时刻的频率分布如下：明确问题：他能及时赶上火车的概率是多少？i）分析法：用概率统计知识建立分析模型，求解析解。（思考）ii）模拟法：用概率统计知识建立模型，通过模拟求近似解。即先建立模拟模型，然后通过计算机模拟得到问题的近似解。在同样条件下多次试验，计算他能及时赶上火车的频率。建模方向（思路）：问题分析：能及时赶上火车的充要条件是：其中T1—火车从A站出发的时刻；T2—火车的运行时间；T3—他到达B站的时刻。是什么变量？如何模拟？213213TTTTTT即ii）将午后1时记为t=0，设火车运行时间T2服从正态分布：T2~N（30，22）。基本假设：i）假设T1，T2，T3都是相互独立的随机变量;火车出发时刻T1和人到达B站时刻T3的分布律分别为：T1（分）0510P(t)0.70.20.1T3（分）28303234P(t)0.30.40.20.1建立模型：为了简化计算，将下午1点记为初始时刻。得到随机变量T1和T3的分布律如下：能及时赶上火车的概率p=P{T3＜T1＋T2}如果r为在(0,1)均匀分布的随机数，为了模拟随机变量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可分别用来模拟随机变量T1和T3。主要变量说明：n模拟次数k临时变量，存储当前累计模拟次数count存储赶上火车的次数两种不同风格的算法描述模拟算法设计输入：赶火车次数（天数）输出：赶上火车的频率第1步输入模拟次数n第2步k=1，count=0第3步当k=n，执行Step4－9，否则执行第11步第4步生成均匀分布随机数赋给r第5步由r及公式确定T1模拟火车出发时刻第6步生成均匀分布随机数赋给r；第7步由r及公式确定T3模拟人达到时刻第8步生成正态分布随机数T2模拟火车运行时间第9步IFT1+T2T3,count=count+1，END第10步执行第3步第11步输出赶上火车频率p=count/n算法1（分步骤描述）i）初始化：输入模拟次数n;count=0;ii)模拟n次fori=1ton,模拟随机变量T1,T2,T3，分别赋给t1,t2,t3;ift1+t2t3,count=count+1endifendforapp_prb=count/n;算法2（伪代码描述）n=input('输入模拟次数：');count=0;fori=1:n,rt1=rand;%模拟随机变量t1（火车从A站出发的时刻）ifrt10.7T1=0;elseifrt1=0.7&rt10.9T1=5;elseT1=10;endT2=30+randn*2;%模拟随机变量t2（火车的运行时间）%模拟随机变量t3（他到达B站的时刻）rt3=rand;ifrt30.3T3=28;elseifrt3=0.3&rt30.7T3=30;elseifrt3=0.7&rt30.9T3=32;elseT3=34;endifT3T1+T2,count=count+1;endend%forprob=count/n模拟程序序号模拟次数近似概率p110000.6280210000.6920310000.6530450000.6490550000.6260650000.6288模拟结果：每次模拟1000次或5000次系统模拟注意事项：一次模拟结果毫无意义！模拟是试验性的,是思维结果的验证。必须进行足够多次的模拟，并对结果进行统计分析。系统模拟特点：系统模拟是研究系统，特别是动态系统的重要方法，对于：结构复杂的系统；很难用解析方法求出变量关系的系统；内部机理不明的“黑箱”系统；为验证用其他方法建立的模型及结果,应是较好的选择。1.3随机变量的建模•利用理论分布，基于对问题的实际的、合理的假设，选择适当的理论分布模拟随机变量优点是给出了各种理论结果出现的概率，便于进行数学分析和处理。但此方法仅限于十分简单的情况，当问题越复杂，数学处理变得越困难，并且丢失了试验数据的信息。•基于实际数据的频率做近似模拟优点在于完全与观察数据相符，并且随实际问题的复杂程度增大不会产生更大的困难，仅增大工作量而已。缺点是不便于进行数学分析，不得不依赖于模拟得到的统计结果。•应用常将两种模拟方法结合起来使用1.4均匀分布随机变量模拟1.4.1平方取中法1.4.2线性同余法线性同余法程序：产生0－c整数1.5Matlab随机模拟函数：1.5.1常见分布随机变量的模拟1.5.2其它随机变量的模拟2Monte-Carlo模拟2.1Monte-Carlo原理2.2蒙特卡罗法应用2.2.1求解非线性规划框图初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1PMAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNjn?f(X)≥Q?YNX*=X,Q=f(X)k=k+1kMAXK?YN输出X,Q,停止YN在MATLAB软件包中编程，共需3个Ｍ文件:randlp.m,mylp.m,lpconst.m.主程序为randlp.m.%mylp.mfunctionz=mylp(x)%目标函数z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);%转化为求最小值问题%lpconst.mfunctionlpc=lpconst(x)%约束条件If3*x(1)+x(2)-10=0.5&3*x(1)+x(2)-10=-0.5%约束条件的误差为±0.5lpc=1;elselpc=0;end%randlp.mfunction[sol,r1,r2]=randlp(a,b,n)%随机模拟解非线性规划debug=1;a=0;%试验点下界b=10;%试验点上界n=1000;%试验点个数r1=unifrnd(a,b,n,1);%n1阶的［a,b］均匀分布随机数矩阵r2=unifrnd(a,b,n,1);sol=[r1(1)r2(1)];z0=inf;fori=1:nx1=r1(i);x2=r2(i);lpc=lpconst([x1x2])；iflpc==1z=mylp([x1x2])；ifzz0z0=z；sol=[x1x2]；endendend2.2.2估算圆周率2.2.3估算定积分1）频率法2）平均值法例1：某港口提供有足够的泊位供船舶停靠，但现在仅有一个可供装卸的泊位，船舶先到则先进行装卸，如果船舶得不到及时装卸而造成的滞期费为每小时100元。现要弄清该系统的性能，重点考察船舶进入该港后等待装卸的滞留时间以及卸位（即装卸用的泊位）的利用率，从而进行经济效益分析。3模拟模型案例到达间隔(小时)56789101112131415161718到达船舶数(艘)136791011111197654首先，对进入该港口的100艘船进行实际调查，记录其活动情况，得到这100艘船到达港口的时间间隔和装卸时间的分布情况的频数和累积频率分布。装卸时间（小时）910111213141516装卸船舶数（艘）2022191610832表1100艘船到达港口的时间间隔频数表表2100艘船装卸时间频数表到达间隔567891011121314151