1、计算机仿真、随机模拟,基本内容：,1、计算机仿真的基本概念,2、随机数的产生,3、时间步长法,4、事件步长法,6、实验作业,5、Monte Carlo方法,一、计算机仿真的基本概念,仿真：,就是将所研究的对象用其它手段加以模仿的一种活动。,实物仿真,非实物仿真,如曹冲称象、军事演习、风洞试验等,计算机仿真：,是一种非实物仿真方法，通过建立数学模型、编制计算机程序实现对真实系统的模拟。从而了解系统随时间变化的行为或特性, 以评价或预测一个系统的行为效果,为决策提供信息的一种方法.它是解决较复杂的实际问题的一条有效途径。,鼠疫的检测和预报，三峡的安全、生态，公交车的调度,航空管理, 经营投资，道路
2、的修建, 通信网络服务，电梯系统服务等等。,计算机仿真举例：,难以用数学公式表示的系统，或者没有建立和求解数学模型的有效方法 虽然可以用解析的方法解决问题，但数学的分析与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法 希望能在较短的时间内观察到系统发展的全过程，以估计某些参数对系统行为的影响 难以在实际环境中进行实验和观察时，计算机仿真是唯一可行的方法，例如太空飞行的研究 需要对系统或过程进行长期运行比较，从大量方案中寻找最优方案,适用计算机仿真解决的问题：,便于重复进行试验，便于控制参数，时间短，代价小；,计算机仿真的优点：,可以在真实系统建立起来之前，预测其行为效果，从而可以从不同
3、结构或不同参数的模型的结果比较之中，选择最佳模型；,对于缺少解析表示的系统，或虽有解析表示但无法精确求解的系统，可以通过仿真获得系统运行的数值结果；,对于随机性系统，可以通过大量的重复试验，获得其平均意义上的特性指标。,仿真常用术语:,系统：,一些具有特定功能相互之间以一定的规律联系着的物体所组成的总体,系统的对象、系统的组成元素都可以称为实体.,实体：,属性是对实体特征的描述，可以是文字型、数字型或逻辑型,属性:,实体在一段时间内持续进行的操作或活动,活动:,例如：银行系统，销存系统,改变系统状态的瞬间变化的事情.,事件:,为了限制所研究问题涉及的范围, 用系统边界把所研究的系统与影响系统的
4、环境区分开来.,系统边界：,事件表一般是一个有序的记录列，每个记录包括事件发生时间、事件类型等一些内容,事件表:,系统的状态是指在某一时刻实体及其属性值的集合,状态:,表示仿真时间的变量,仿真时钟:,仿 真 研 究 步 骤,系统分析,模型构造,模型运行,输出结果,明确问题和提出总体方案: 把被仿真系统的内容表达清楚； 弄清仿真的目的、系统的边界； 确定问题的目标函数和可控变量； 找出系统的实体、属性和活动等。,系统分析,建立模型；选择合适的仿真方法（如时间步长法、事件表法等）；确定系统的初始状态；设计整个系统的仿真流程图。 收集数据； 编写程序、程序验证； 模型确认。,模型构造,运行：确定具体
5、的运行方案,如初始条件、参数、步长、重复次数等,然后输入数据,运行程序。 改进：将得出的仿真结果与实际系统比较,进一步分析和改进模型,直到符合实际系统的要求及精度为止。,模型的运行与改进,设计出结构清晰的仿真结果输出。包括提供文件的清单,记录重要的中间结果等。输出格式要有利于用户了解整个仿真过程 ,分析和使用仿真结果.,设计格式输出仿真结果,(库存问题)某电动车行的仓库管理人员采取一种简单的订货策略，当库存量降低到P辆电动车时就向厂家订货，每次订货Q辆，如果某一天的需求量超过了库存量，商店就有销售损失和信誉损失，但如果库存量过多，会导致资金积压和保管费增加。若现在已有如下表所示的两种库存策略，
6、试比较选择一种策略以使总费用最少。,计算机仿真举例：,这个问题的已知条件是： (1)从发出订货到收到货物需隔3天。 (2)每辆电动车保管费为0.50元/天，每辆电动车的缺货损失为1.60元/天，每次的订货费为75元。 (3)每天电动车需求量是0到99之间均匀分布的随机数。 (4)原始库存为110辆，并假设第一天没有发出订货。,分析：这一问题用解析法讨论比较麻烦，但用计算机按天仿真仓库货物的变动情况却很方便。我们以30天为例，依次对这两种方案进行仿真，最后比较各方案的总费用，从而就可以做出决策。,计算机仿真时的工作流程是早上到货、全天销售、晚上订货，以一天为时间步长进行仿真。,首先检查这一天是否
7、为预定到货日期，如果是，则原有库存量加Q，并把预定到货量清为零；如果不是，则库存量不变。,（程序附后）,接着仿真随机需求量，这可用计算机语言中的随机函数得到。若库存量大于需求量，则新的库存量减去需求量；反之，则新库存量变为零，并且要在总费用上加缺货损失。,然后检查实际库存量加上预定到货量是否小于重新订货点P，如果是，则需要重新订货，这时就加一次订货费。如此重复运行30天，即可得所需费用总值。,由此比较这两种方案的总费用，可以得到最好的方案。,1均匀分布的随机数及其产生,对随机现象进行模拟，实质上是要给出随机变量的模拟也就是说利用计算机随机地产生一系列数值，它们的出现服从一定的概率分布,则称这些
8、数值为随机数。,最常用的是在(0,1)区间内均匀分布的随机数，也就是我们得到的这组数值可以看作是(0,l)区间内均匀分布的随机变量的一组独立的样本值,其它分布的随机数可利用均匀分布的随机数产生,二、随机数的产生,产生模拟随机数的计算机命令,在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：,1、产生m*n阶(0,1)均匀分布的随机数矩阵：rand (m, n)产生一个(0,1)均匀分布的随机数：rand,2、产生m*n阶(a, b)均匀分布U(a, b)的随机数矩阵：unifrnd (a, b, m, n)产生一个(a, b)均匀分布的随机数：unifrnd(a, b),ran
9、dn(m,n)： m*n阶N(0,1)标准正态分布随机数矩阵,3、正态分布,To Matlab（rnd）,若连续型随机变量X 的概率密度函数为其中 0为常数，则称X 服从参数为 的指数分布。,指数分布的期望值为,电子元件的寿命，电话的通话时间，微生物的寿命，随机服务系统中的服务时间都服从指数分布。,指数分布在排队论、可靠性分析中有广泛应用。,注意：Matlab中，产生参数为 的指数分布的 命令为exprnd( ).,例 顾客到达某商店的间隔时间服从参数为0.1的指数分布,指数分布的均值为1/0.1=10。指两个顾客到达商店的平均间隔时间是10个单位时间.即平均10个单位时间到达1个顾客. 顾客
10、到达的间隔时间可用exprnd(10)模拟。,设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为其中 0为常数，则称X服从参数为 的泊松分布。,容器内的细菌数，十字路口的交通事故，寻呼台的寻呼次数及每天到商店购买商品的顾客数等都服从泊松分布。,泊松分布的期望值为,指数分布与泊松分布的关系：,(1)指两个顾客到达商店的平均间隔时间是10个单位时间.即平均10个单位时间到达1个顾客.(2)指一个单位时间内平均到达0.1个顾客,例 (1)顾客到达某商店的间隔时间服从参数为0.1的指数分布 (2)该商店在单位时间内到达的顾客数服从参数为0.1的泊松分布,例2 敌坦克分队对我方阵地实施突袭
11、，其到达规律服从泊松分布，平均每分钟到达辆（1）模拟敌坦克在分钟内到达目标区的数量，以及在第、分钟内各到达几辆坦克（2）模拟在3分钟内每辆敌坦克的到达时刻。,（1）用poissrnd(4)进行模拟,To Matlab（poiss）,（2）坦克到达的间隔时间应服从参数为4的指数分布，用exprnd（1/4）模拟。,To Matlab（time）,2随机变量的模拟,利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟,(1)离散型随机变量的模拟设随机变量X的分布律为P (X=xi)=pi,i=1,2,令p(0)=0, p(n)= pi, n=1,2,将p(
12、n)作为分点，把区间(0 , 1)分为一系列小区间(p(n-1), p(n)对于均匀的随机变量RU(0,l)，则有,P( p(n-1)R p(n)= p(n)-p(n-1)= pn ,n=1,2,由此可知，事件(p(n-1)R p(n)和事件(X=xn) 有相同的发生的概率因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是：每产生一个(0,1)上均匀分布的随机数r，若p(n-1)r p(n)则理解为发生事件“X=xn”,例3 随机变量 x = 0，1，2表示每分钟到达超市收款台的人数，有分布列xk 0 1 2pk 0.4 0.3 0.3 模拟十分钟内顾客
13、到达收款台的状况.,To Matlab（paidui）,仿真分为静态仿真和动态仿真。动态仿真可分为连续系统仿真和离散系统仿真。离散系统是指状态变量只在某个离散时间点集合上发生变化的系统。例如：电梯系统服务，排队系统，通信网络服务的仿真等。连续系统是指状态变量随时间连续改变的系统。例如：传染病的检测和预报等。,为了仿真系统必须设置一个仿真时钟将时间从一个时刻向另一个时刻推进，并且可随时反映系统时间的当前值。模拟时间推进方式有两种：时间步长法和事件步长法。模拟离散系统常用事件步长法。连续系统常用时间步长法（也称固定增量推进法或步进式推进）。,三、时间步长法,时间步长法:,按照时间流逝的顺序，一步一
14、步地对系统的活动进行仿真。在整个仿真过程中，时间步长固定不变。,基本思想:,在进行系统仿真的过程中，可以把整个过程分成许多相等的时间间隔，时间步长的长度可以根据实际问题分别取作秒,分,小时，天等。程序中按照这个步长前进的时钟就是仿真的时钟。,基本步骤:,首先选取对象系统的一个初始起点作为仿真时钟的零点，然后根据实际问题的需要，选定一个时间步长。于是从仿真时钟的零点开始，每推进一个时间步长就对系统的活动和状态按照预定的规则和目的进行考察、分析、计算和记录，直到预定仿真结束时刻为止。,时间步长法,应用举例-池水含盐量问题,例6:池水含盐量问题某水池有2000m3水，其中含盐2kg，以每分钟6m3的
15、速率向水池内注入含盐率为0.5kg/m3的盐水，同时又以每分钟4m3的速率从水池流出搅拌均匀的盐水试用计算机仿真该水池内盐水的变化过程，并每隔10min计算水池中水的体积、含盐量和含盐率欲使池中盐水的含盐率达到0.2kg/m3，需经过多少时间？,池水含盐量随时间变化的示意图,分析:系统中，实体是水，属性是水的体积、含盐量和含盐率，活动是水的注入与流出，由于注入和流出活动的作用，池中水的体积与含盐量、含盐率均随时间变化，初始时刻含盐率为0.001kg/m3,以后每分钟注入含盐率为0.5 kg/m3的水6m3，流出混合均匀的盐水4m3,当池中水的含盐率达到0.2kg/m3时,结束仿真过程.,为了能
16、定量地考察系统在任一时刻的属性，引入下列记号：注水速度： WI6 m3/min排水速度： WO4 m3/min注入水的含盐率：SI0.5 kg/m3最终含盐率：SF0.2 kg/m3T时刻水的体积：VTm3T时刻水的含盐量：ST kgT时刻水的含盐率：SR=ST/VT kg/m3,对于这样一个连续系统仿真时，必须在一系列离散时间点t0 t1 t2 tn上来进行考察，这些离散时间点之间的间隔T= ti- ti-1(i=1,2,n)就是时间步长若取步长为1min时，就要每隔1min考察一次系统的状态特性，并作出相应的记录与分析在本题中每隔1min池水的动态变化过程是：每分钟水的体积增加 6-42
17、m3；每分钟向池内注入盐 60.5 3 kg；每分钟向池外流出盐 4SR kg；每分钟池内增加盐 3-4SR kg 取T0作为系统仿真的初始时刻，取T=0时系统的属性为系统的初始状态取池中盐水的含盐率达到0.2kg/m3的时刻作为仿真的终止时刻,池水含盐量仿真流程,初始化,按时间步长前进1min,计算池水体积VT,含盐量ST,含盐率SR,是否到达10min？,结束,输出结果,是,否,含盐率达到SF?,打印次数加1,输出时间T,水体积VT 含盐量ST,含盐率SR,下一个10min开始，计时单元清零,是,否,池水含盐量仿真结果,在这个问题中，池水的体积和含盐量、含盐率都是随时间连续变化的故这是一个
18、连续系统的仿真问题，我们是把这个连续过程离散化后用时间步长法进行仿真的我们还可以用微分方程建立这个问题的数学模型，并求出它的解析解，这个解析解就是问题的精确解有兴趣可以按照这一思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异也可进一步调整前面仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真精度的影响,四、事件步长法,事件步长法:,是以事件发生的时间为增量，按照事件发生的时间顺序，一步一步地对系统的行为进行仿真，直到预定的时间结束为止。,事件步长法,初始状态,事件步长加1,在当前步长内， 考察分析，计算和 记录系统的活动,仿真时间到否？,仿真结束,输出结果,是,否,事件步
19、长法与时间步长法的主要区别是：,(1)事件步长法与时间步长法都是以时间为增量来考察系统状态的变化但在时间步长法中，仿真时钟以等步长前进，而在事件步长法中，仿真时钟的步长取决于事件之间的间隔,(2)时间步长法在一个步长内，认为系统所处的状态相同，因而所选步长的大小将影响仿真的精度而在事件步长法中，每个事件的发生均有确切的时刻，不需要人为地选取步长，步长的大小对仿真精度影响较小,(3) 时间步长法每前进一个步长就要对整个系统进行一次全面考察，即使状态没有发生变化时也要扫描，而事件步长法只是在某一事件点上判断和比较事件是否出现因此，一般地讲，当判断比较的数目较大时，用时间步长法可以节省用机时间，而当
20、相继两个事件出现的平均间隔较长时，更适合于用事件步长法,例7 收款台前的排队过程的仿真假设：1、顾客的到达收款台是随机的，平均时间间隔为0.5分钟，即间隔时间服从=2的指数分布。2、对不同的顾客收款和装袋的时间服从正态分布N(1，1/9)。现要求模拟20位顾客到收款台前的排队情况， 主要关心问题：每个顾客的平均等待时间,队长,服务 员的工作效率.,收款排队系统主控程序图,顾客到达子程序图,服务结束子程序图,假设 t(i): 第i位顾客到达时刻，t2(i):第i位顾客受到的服务时间(随机变量)， T(i): 第i位顾客离去时刻. 将第i位顾客到达作为第i件事发生;t(i+1)- t(i)= r(
21、i) (随机变量) 平衡关系： 当 t(i+1)T(i) 时, T(i+1)=t(i+1)+t2(i+1); 否则, T(i+1)=T(i)+t2(i+1),Matlab(kefu)程序,五、Monte Carlo方法,基本思想：,Monte Carlo方法亦称为随机模拟（Random simulation）方法，有时也称作随机抽样技术或统计试验方法。此方法对研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数,当试验次数充分多时，某一事件出现的频率近似于该事件发生的概率.,例8：用蒙卡模拟求圆周率的估计值。,设二维随机变量(X,Y)在正方形服从均匀分布，,则(X,Y)
22、落在四分之一圆内的概率 为 PX2+Y2=1= /4,现产生n对二维随机点(xi,yi), xi,yi是区间0,1上均匀分布的随机数，若其中有k对满足 xi2+yi2=1，即相当于做n次投点试验，其中有k点落在1/4圆内。计算落入1/4圆内的频率为k/n，由大数定律，事件发生的频率依概率收敛于发生的概率，可得圆周率的估计值为4k/n。次数越多，精度越大。,Matlab(yuanzhoulv),例9：用Monte Carlo方法计算积分,分析：任取一列相互独立的、都具有a,b中均匀 分布的随机变量xi,则g(xi)也是一列相互独立的 随机变量,而且：,所以,只要求出,便能得到J 的数值。,为求E
23、g(xi),由大数定理知：,应用举例-报童的策略,例10 报童每天清晨从报社购进报纸零售,晚上将没有卖掉的报纸退回.每份报纸的购进价为1.3元,零售价为2元,退回价为0.2元.报童售出一份报纸赚0.7元 ,退回一份报纸赔1.1元.报童每天如果购进的报纸太少,不够卖时会少赚钱,如果购得太多卖不完时要赔钱. 试为报童筹划每天应如何确定购进的报纸数使得收益最大.报纸每捆10张，只能整捆购买，报纸可以分为3种类型的新闻日：好、一般、差，它们的概率分别为0.35,0.45和0.2,在这些新闻日中每天对报纸的需求分布的统计结果如下表：,试确定每天报童应该订购的报纸数量。,解：我们通过计算机仿真来解决此问题
24、。最优策略 应该是每天的利润最大。利润=销售收入-报纸成本+残值 这是一个随机现象的计算机仿真问题， 故先确定各种情况的随机数的对应关系。 新闻日和需求量对应的随机数分别如下面两个表格所示,计算机仿真的流程： 1）令每天的报纸订购数变化，40-100； 2）让时间从1开始变化（循环）到365； 3）产生新闻种类的随机数，确定当天的新闻类型； 4）产生需求量随机数，确定当天的报纸需求量； 5）计算当天的收入，计算累积利润， 6）比较得出最优定货量。,Matlab编程(baotong程序) x1=rand(365,1); x2=rand(365,1); for n=4:10paper=n*10;s
25、b(n)=0;for i=1:365if x1(i)0.35if x2(i)0.03 news=40;elseif x2(i)0.08news=50;elseif x2(i)0.23news=60;elseif x2(i)0.43news=70;,elseif x2(i)0.78news=80;elseif x2(i)0.93news=90;else news=100;end elseif x1(i)0.8if x2(i)0.10news=40;elseif x2(i)0.28news=50;elseif x2(i)0.68news=60;elseif x2(i)0.88news=70;,el
26、se news=80;end endif paper=newssale=news;remand=paper-news;elsesale=paper;remand=0;endsb(n)=sb(n)+2*sale- 1.3*paper+0.2*remand;end end,elseif x2(i)0.96news=80;else news=90;endelseif x2(i)0.44news=40;elseif x2(i)0.66news=50;elseif x2(i)0.82news=60elseif x2(i)0.94news=70;,optnews=40; optmoney=sb(4); 4
27、0,sb(4)/365 for n=5:10if sb(n)=optmoneyoptnews=n*10;optmoney=sb(n);endn,sb(n)/365 end optnews,optmoney,optmoney/365经过计算机仿真后得到最优购货量是每天60份， 平均每天利润34.4元。,应用举例-赶火车仿真,例11 赶火车过程仿真 一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站运行时间为均值30分钟、标准差为2分钟的正态随机变量.火车大约在下午1点离开A站。离开时刻的频率分布为,用计算机仿真火车开出、火车到达B站、这个人到 达B站情况，并给出能赶
28、上火车的仿真结果。,引入以下变量：T1 火车从A站开出的时刻；T2 火车从A站运行到B站所需要的时间；T3 此人到达B站的时刻；,这个人到达B站时的频率分布为：,开车时间的仿真测试s1=0; s2=0; s3=0;x=rand(10000,1);for i=1:10000if x(i)0.9s3=s3+1;endend s1/10000, 1-s1/10000-s3/10000,s3/10000,s1=0; s2=0; s3=0;s4=0;x=rand(10000,1);for i=1:10000if x(i)0.3s1=s1+1;elseif x(i)0.7(前面已运行过0.3的，系统默认会
29、踢出0.3的）s2=s2+1;elseif x(i)0.9s3=s3+1;else s4=s4+1;endend end s1/10000, s2/10000（默认踢出前一项）,s3/10000,s4/10000,人到达时刻仿真测试,火车运行时间的仿真测试 x=randn(10000,1); for i=1:10000y(i)=30+2*x(i); end,赶上火车的仿真结果 s=0; x1=rand(10000,1);（颠倒写也可） x2=rand(10000,1); x3=randn(10000,1); for i=1:10000 if x1(i)0.7 T1=0; elseif x1(i
30、)0.9 T1=5;elseT1=10;end T2=30+2*x3(i);,if x2(i)0.3T3=28; elseif x2(i)0.7T3=30;elseif x2(i)0.9T3=32;else T3=34;end if T3T1+T2s=s+1; end end s/10000,1、(射击命中率) 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人现模拟我方将要对敌人实施的20次打击结果，并确定有效射击的比率及毁伤敌方火炮的平均值。,实验作业,2、已知零件C由零件A和零件B连接而成，已知A、B的长度均为随机变量，具体数值如下。试抽取10个样本以计算C的平均长度。,
