蒙特卡罗模拟应用

•一.随机现象的模拟•1.随机变量及其分布•随机事件：在一定条件下有可能发生的事件。•概率：随机事件发生的可能性的度量P(A)，•0≤P(A)≤1.•随机变量：在一定的范围内随机取值的变量,•“X=ak”(k=1,2,…,n),或随机发生.•badttpbaXP)()],[(],[baX•随机变量的分布：•若已知则称•aka1a2…an•P(“X=ak”)p1p2…pn•为随机变量X的分布列,简称X的分布•若已知•则称p(x)为随机变量X的分布密度,简称X的分布1,)(1nkkkkppaXPbadttpbaXP)()],[(1)(dttp•称为两点分布•单次贝努里试验的结果•称为二项分布•n重贝努里试验的结果•称为离散的均匀分布•以相同的概率取所有可能的数值•称为泊松分布•发生率较低次数无限增大时贝努里试验的极限））（（222exp21)(xxP)1(,)1()(pppCaXPknkknknknaXPk,,2,1,/1)(,2,1,!)(kekkXPk011)(xpxpxP•称为均匀分布密度•在[a,b]的任何相等的子区间上取值的概率相同•称为正态分布密度•许多偶然因素作用结果的总和。•称为指数分布密度质点于随机时间陆续到达的时间间隔，其它,0],[,1)(baxabxP））（（222exp21)(xxP0,00,)(xxexpx•2.随机数和随机变量的模拟•10.随机数（RND）:•计算机随机数发生器产生的数串,•它在（0,1)中的分布是均匀的。•一般称之为伪随机数。•20.具特定分布随机变量的模拟•变量X有分布列•令•则有1,)(1nkkkkppaXP1,)()1()(nkkppp0,)0(1)(pppkiik1,)()1()(nkkppp•以p(k)为分点，将[0，1]分为n个小区间•取随机数R,则容易证明P(“p(k-1)Rp(k)”)=pk=P(“X=ak”)•随机事件“p(k-1)Rp(k)”与“X=ak”有相同的概率分布。•可以使用随机数在各小区间出现的情况•来模拟随机事件“X=ak”发生的状况。•例随机变量x={0，1，2}表示每分钟到达超市收款台的人数，有分布列•xk012•pk0.40.30.3•模拟十分钟内顾客到达收款台的状况•用MATLAB模拟随机事件•rand(20,1)%生成20个均匀随机数向量。•randn(20,1)%生成20个正态随机数向量。•r=rand(1,10);•fori=1:10;•ifr(i)0.4•n(i)=0;•elseif0.4=r(i)&r(i)0.7•n(i)=1;•elsen(i)=2;•end;•end•r•r=0.56780.79420.05920.60290.0503•0.45650.01850.82140.44470.6154•n•n=1201010211•r•r=0.23110.60680.48600.89130.7621•0.79190.92180.73820.17630.4057•n•n=0112222201•二.系统仿真（Simulation）•1.系统仿真：使用计算机对一个系统的结构和行为进行动态模拟。•为决策提供必要的参考信息。•特点：对象真实、复杂，进行模仿。•2.仿真模型：由计算机程序控制运行•从数值上模仿实际系统的动态行为。•3.仿真过程•1.现实系统的分析：•了解背景，明确目的，提出总体方案。•2.组建模型：•确定变量,明确关系,设计流程,编制程序•3.运行检验：•确定初始状态，参量数值，•运行程序，检验结果，改进模型。•4.输出结果•三.动态系统的仿真•1.时间步长法：•把整个仿真过程分为许多相等的时间间隔•每个间隔为一个时间单位—时间步长。•在每个时间步长内模拟系统的动态。•仿真时钟：用以控制时间步进的过程•（每一次步进一个步长）•例3.15池水含盐•池中有水2000m3，含盐2kg，•以6m3/分的速率向池中注入浓度为0.5kg/m3的盐水，•又以4m3/分的速率从池中流出混合后的盐水•问欲使池中盐水浓度达到0.2kg/m3，需要多长时间？••系统分析：•池中有盐水，•匀速注入浓盐水，•匀速流出混合后的盐水，•池中盐水的浓度变化。•目的：仿真池中盐水浓度的变化，给出达到给定浓度的时间。•变量、参量•时间t，体积V(t),盐量S(t),浓度p(t);•流入流速rI,流入浓度pI,•流出流速rO,流出浓度p(t),给定浓度p*•时间步长Δt,打印步长T.•关系:在[t,t+Δt]内有trrtVttVOI)()()(ttprprtSttSOII)]([)()()(/)()(ttVttSttp•动态系统模拟的伪代码•运算池水含盐动态系统模拟•变量V(n)=时刻n池中盐水体积•p(n)=时刻n池中盐水浓度•S(n)=时刻n池中盐水含盐量•Δt=时间单位•N=模拟时间长度•输入Δt，V(0),p(0),S(0),N•过程Begin•forn=1toNdo•Begin•V(n+1)←V(n)+(rI-r0)Δt•S(n+1)←S(n)+[ripi-r0p(n)]Δt•p(n+1)←S(n+1)/V(n+1)•End•End•输出V(1),V(2),…,V(n)•S(1),S(2),…,S(n)•p(1),p(2),…,p(n)系统仿真流程图•N•Y初始化V(0),S(0)仿真时钟t=0打印时钟T=0计算V(t+Δt),S(t+Δt),p(t+Δt)时钟步进t=t+Δt,T=T+1p(t)p*Tm打印输出•参数rI=6，rO=4，pI=0.5，p*=0.2•初始V(0)=2000,S(0)=2•步长Δt=1,m=10•结果表3.15，t=185分,盐水浓度为0.2.•模型IIOprtStVrdtdS)()()0()(,VrttVrrrdtdVOIIIOprtStVrdtdS)()(clft=1;v=[2000];s=[2];p=[1/1000];%初态V=[v(end)];S=[s(end)];P=[p(end)];x=[0];%终态whilep(end)=0.2%最终阈值T=0;whileT10%打印阈值T=T+1;t=t+1;%时钟步进v=[v1];s=[s1];p=[p0];%变量步进v(t)=v(t-1)+2;s(t)=s(t-1)+3-4*p(t-1);%模型p(end)=s(end)/v(end);ifp(end)0.2T=20;%转输出end;end;x=[xt-1];%时间记录V=[Vv(end)];S=[Ss(end)];P=[Pp(end)];end;V1=10^3.*V;a=[x‘,V1’,S‘,P’]%输出调节•MATLAB实现•x’=x/2+1，0≤t≤10,x(0)=0•建立M文件fun.M•functiony=fun(t,x)•y=x/2+1;•运行结果•t0=0;tf=10;•x0=0;•[t,x]=ode23(‘fun’,t0,tf,x0);•plot(t,x)•MATLAB实现•y’’+(y2-1)y’-y=0,y(0)=0.25,y’(0)=0•化为一阶方程组•y1’=y1(1-y22)-y2,y2’=y1,•建立M文件rhf.M•functionfyy=rhf(t,x)•fyy=[y(1).*(1-y(2).^2)-y(2);y(1)];•运行结果•t0=0;tfinal=15;Y0=[0.25,0]’;•[T,Y]=ode23(‘rhf’,t0,tfinal,Y0);•plot(T,Y)求表达式（符号运算）S=dsolve(‘Dx=(3-6*x)/(2000+2*t)’);求数值解建立M文件fun.Mfunctiony=fun(t,x)y=(3-6*x)/(2000+2*t);t0=0;tf=200;x0=2;[t,x]=ode23(‘fun’,[t0tf],x0);Plot(t,x);drdrVtVpptptptrdttdptVOtIII)]()([)()0()]()()[()()(000••问题•1.在池水含盐的问题中令•rO=rI=6m3/分•10.池中盐水的浓度如何变化？•20.若当p(t)=0.3kg/m3时令pI=0,•需要多少时间达到p*=0.2kg/m3?•2.若池中盐水的初始浓度为p*，•对于不同的初始体积V0，当pI=0时，•计算池中盐水浓度降低一半所用的时间•例3.16市场服务•超市有两个出口的收款台，•两项服务：收款、装袋。•两名职工在出口处工作。•有两种安排方案：•开一个出口，一人收款、一人装袋；•开两个出口，每个人既收款又装袋。•问商店经理应选择哪一种收款台的服务方案。••假设：•1.顾客的到达是随机的。•2.收款装袋的时间是相同的。•3.第一种方案中，收款与装袋同时进行。•参量、变量：•n(t)到达顾客人数，L(t)队列长,τ服务时间•T1排队时间累加,T2服务时间累加•模型•n(t)=0,L(t)=0则L(t+Δt)=0,T2(t+Δt)=T2(t)•L(t)0则L(t+Δt)=L(t)+n-1,T2(t+Δt)=•n(t)=1,L(t)=0则L(t+Δt)=0,T2(t+Δt)=T2(t)+τ•L(t)0则L(t+Δt)=L(t)+n-1,T2(t+Δt)=•n(t)=2,L(t)=0则L(t+Δt)=n-1,T2(t+Δt)=•L(t)0则L(t+Δt)=L(t)+n-1,T2(t+Δt)=clfL=zeros(1,31);%L等待的顾客人数，T1=zeros(1,31);%T1等待时间的累加，T2=zeros(1,31);%T2服务时间的累加,L1=zeros(1,31);%L1到达顾客人数累加。t=1;tau=1;x=0:30;r=rand(1,30);%随机数%随机模拟fori=1:30;t=t+1;if0=r(i)&r(i)0.4n=0;elseif0.4=r(i)&r(i)0.7n=1;elsen=2;end;%排队分析ifL(t-1)==0&n==0L(t)=L(t-1);T1(t)=T1(t-1);%模型T2(t)=T2(t-1);L1(t)=L1(t-1);elseL(t)=L(t-1)+n-1;T1(t)=T1(t-1)+L(t);T2(t)=T2(t-1)+tau;L1(t)=L1(t-1)+n;end;end;r=[0r];a=[x',r',L',L1',T1',T2']eL=T2(end)/tau%已被服务的人数L2=(find(L1eL))L3=sum(L(L2))%未被服务的顾客等待时间总和g1=(T1(end)-L3)/eL%平均等待时间g2=g1+tau%平均逗留时间g3=eL/30%平均每分钟服务的顾客人数•tRNDnLT1T2•1.812111•2.020001•3.601001•4.040000•5.461001•6.310000•7.671001•8.250000•9.240000•10.100000•11.401001•12.020000tRNDnLT1T213.39000014.68100115.08000016.59100117.66100118.90211119.12000120.64100121.79211122.31000123.86211124.681111•2.事件表法（面向事件法）：每处理一个事件就前进一步（每步的时间可能不同），以事件为中心安排。•对系统中的一系列不同类型的事件按发生的前后顺序逐个进行分析、处理。•事件表：记录今后将要发生的事件及其属性（发生的时间，事件的类型等），调度事件执行的顺序。•它随着程序的执行过程不断地更新和补充以保证仿真过程有序地进行。系