《蒙特卡罗模拟应用》由会员分享，可在线阅读，更多相关《蒙特卡罗模拟应用（50页珍藏版）》请在金锄头文库上搜索。
1、 一. 随机现象的模拟 1. 随机变量及其分布 随机事件：在一定条件下有可能发生的事件。 概率：随机事件发生的可能性的度量 P(A)， 0 P(A) 1. 随机变量：在一定的范围内随机取值的变量, “X=ak” (k=1,2,n), 或 随机发生. 随机变量的分布： 若已知 则称 ak a1 a2 an P(“X=ak”) p1 p2 pn为随机变量 X 的分布列, 简称 X 的分布 若已知则称 p(x) 为随机变量 X 的分布密度, 简称X的分布 称 为两点分布单次贝努里试验的结果称 为二项分布n 重贝努里试验的结果称 为离散的均匀分布以相同的概率取所有可能的数值称 为泊松分布发生率较低次数
2、无限增大时贝努里试验的极限称 为均匀分布密度在a,b的任何相等的子区间上取值的概率相同称 为正态分布密度许多偶然因素作用结果的总和。称 为指数分布密度 质点于随机时间陆续到达的时间间隔， 2. 随机数和随机变量的模拟 10. 随机数（RND）: 计算机随机数发生器产生的数串, 它在（0, 1) 中的分布是均匀的。一般称之为伪随机数。 20. 具特定分布随机变量的模拟 变量 X 有分布列 令 则有 以 p(k)为分点，将0，1分为 n 个小区间取随机数 R , 则容易证明 P( “ p(k-1) R p(k) ” ) = pk = P ( “ X = ak” )随机事件 “ p(k-1) R p
3、(k) ” 与 “ X=ak” 有相同的概率分布。可以使用随机数在各小区间出现的情况来模拟随机事件 “ X=ak” 发生的状况。例 随机变量 x = 0，1，2表示每分钟到达超市收款台的人数，有分布列 xk 0 1 2 pk 0.4 0.3 0.3模拟十分钟内顾客到达收款台的状况用用MATLAB模拟随机事件模拟随机事件 rand(20, 1) %生成生成20个均匀随机数向量。个均匀随机数向量。 randn(20,1) %生成生成20个正态随机数向量。个正态随机数向量。r=rand(1,10);for i=1:10; if r(i)0.4 n(i)=0; elseif 0.4=r(i)&r(i)
4、 rr =0.5678 0.7942 0.0592 0.6029 0.0503 0.4565 0.0185 0.8214 0.4447 0.6154 nn = 1 2 0 1 0 1 0 2 1 1rr =0.2311 0.6068 0.4860 0.8913 0.7621 0.7919 0.9218 0.7382 0.1763 0.4057nn = 0 1 1 2 2 2 2 2 0 1 二. 系统仿真（Simulation） 1. 系统仿真：使用计算机对一个系统的结构和行为进行动态模拟。为决策提供必要的参考信息。 特点：对象真实、复杂，进行模仿。 2. 仿真模型：由计算机程序控制运行从数值
5、上模仿实际系统的动态行为。 3. 仿真过程 1. 现实系统的分析：了解背景，明确目的，提出总体方案。 2. 组建模型：确定变量, 明确关系, 设计流程,编制程序 3. 运行检验： 确定初始状态，参量数值， 运行程序，检验结果，改进模型。 4. 输出结果 三. 动态系统的仿真 1. 时间步长法：把整个仿真过程分为许多相等的时间间隔每个间隔为一个时间单位时间步长。在每个时间步长内模拟系统的动态。 仿真时钟：用以控制时间步进的过程 （每一次步进一个步长） 例 3.15 池水含盐池中有水 2000 m3，含盐 2 kg，以 6m3 / 分 的速率向池中注入浓度为 0.5 kg / m3 的盐水，又以
6、4 m3 / 分的速率从池中流出混合后的盐水问欲使池中盐水浓度达到 0.2 kg / m3，需要多长时间？ 系统分析：池中有盐水，匀速注入浓盐水，匀速流出混合后的盐水，池中盐水的浓度变化。 目的：仿真池中盐水浓度的变化，给出达到给定浓度的时间。 变量、参量 时间 t，体积 V(t), 盐量 S(t), 浓度 p(t); 流入流速 rI, 流入浓度 pI, 流出流速 rO, 流出浓度 p(t), 给定浓度 p* 时间步长 t , 打印步长 T.关系: 在 t, t+ t 内有 动态系统模拟的伪代码运算 池水含盐动态系统模拟变量 V(n)=时刻 n 池中盐水体积 p(n)=时刻 n 池中盐水浓度
7、S(n)=时刻 n 池中盐水含盐量 t = 时间单位 N = 模拟时间长度输入 t，V(0), p(0), S(0), N 过程过程 Begin for n=1 to N do 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,
8、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.模型clft=1; v=2000;s=2;p=1/1000;%初态初态V=v(end); S=s(end); P=p(end); x=0;%终终态态while p(end)=0.2 %最终阈值最终阈值 T=0; while T0.2 T=20;%转输出 end; end; x=x t-1;%时间记录 V=V v(end); S=S s(end); P=P p(end);end; V1
9、=103.*V; a=x,V1,S,P %输出调节 MATLAB实现 x=x/2+1， 0 t 10, x(0)=0建立M文件 fun . M function y=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 function fyy=rhf(t,x) fyy=y(1).*(1-y(2).2)-y(2);y(1);运行结果t0=
10、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 . M function y=fun(t,x) y=(3-6*x)/(2000+2*t);t0=0; tf=200; x0=2; t, x=ode23(fun,t0 tf,x0);Plot(t, x); 问题1. 在池水含盐的问题中令 rO=rI=6m3/分 10. 池中盐水的浓度如何变化？ 20. 若当p(t) = 0.3kg/m3 时令pI=0,
11、需要多少时间达到 p*= 0.2kg / m3?2. 若池中盐水的初始浓度为 p*， 对于不同的初始体积V0，当pI=0时， 计算池中盐水浓度降低一半所用的时间 例 3.16 市场服务 超市有两个出口的收款台，两项服务：收款、装袋。两名职工在出口处工作。有两种安排方案：开一个出口，一人收款、一人装袋；开两个出口，每个人既收款又装袋。问商店经理应选择哪一种收款台的服务方案。 假设： 1. 顾客的到达是随机的。 2. 收款装袋的时间是相同的。 3. 第一种方案中，收款与装袋同时进行。参量、变量： n(t) 到达顾客人数，L(t) 队列长, 服务时间 T1 排队时间累加, T2 服务时间累加模型n(
12、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=ze
13、ros(1,31); %T2服务时间的累加服务时间的累加,L1=zeros(1,31); %L1到达顾客人数累加。到达顾客人数累加。t=1;tau=1;x=0:30;r=rand(1,30);%随机数随机数%随机模拟随机模拟for i=1:30; t=t+1; if 0=r(i) & r(i)0.4 n=0; elseif 0.4=r(i) & r(i)eL)L3=sum(L(L2)%未被服务的顾客等待时间总和未被服务的顾客等待时间总和g1=(T1(end)-L3)/eL %平均等待时间平均等待时间g2=g1+tau %平均逗留时间平均逗留时间g3=eL/30 %平均每分钟服务的顾客人数平均每
14、分钟服务的顾客人数t RND n L T1 T21 .81 2 1 1 12 .02 0 0 0 13 .60 1 0 0 14 .04 0 0 0 05 .46 1 0 0 16 .31 0 0 0 07 .67 1 0 0 18 .25 0 0 0 09 .24 0 0 0 010 .10 0 0 0 011 .40 1 0 0 112 .02 0 0 0 0 t RND n L T1 T213 .39 0 0 0 0 14 .68 1 0 0 115 .08 0 0 0 016 .59 1 0 0 117 .66 1 0 0 118 .90 2 1 1 119 .12 0 0 0 120
15、 .64 1 0 0 121 .79 2 1 1 122 .31 0 0 0 123 .86 2 1 1 124 .68 1 1 1 1 2. 事件表法（面向事件法）：每处理一个事件就前进一步（每步的时间可能不同），以事件为中心安排。对系统中的一系列不同类型的事件按发生的前后顺序逐个进行分析、处理。 事件表：记录今后将要发生的事件及其属性（发生的时间，事件的类型等），调度事件执行的顺序。它随着程序的执行过程不断地更新和补充以保证仿真过程有序地进行。系统仿真初始事件表找出最靠前事件事件属性顾客到达服务结束仿真完毕？输出结果下一顾客到来时刻，记事件表出纳忙？为顾客服务产生服务结束时刻，记事件表统计
16、数据顾客排队有顾客排队？为顾客服务收款台闲置产生服务结束时刻，记事件表统计数据结束服务YNNY例 市场服务假设：1.顾客的到达收款台是随机的，平均时间间隔为0.5分钟，即间隔时间服从=2的(负)指数分布。2.对不同的顾客收款和装袋的时间服从正态分布 N(1，1/3)。参量，变量： t(i): 第i位顾客到达时间， t2(i):第i位顾客受到的服务时间(随机变量)，T(i): 第i位顾客离去时间.将第i位顾客到达作为第i件事发生; t(i+1)- t(i)= t1(i) (随机变量)平衡关系：当 t(i+1)T(i) 时, T(i+1)=t(i+1)+t2(i);否则, T(i+1)=T(i)+
17、t2(i)模拟20位顾客到收款台前的排队情况。clft=zeros(1,21);%每位顾客到达时间T=zeros(1,21);%每位顾客离去时间w=zeros(1,21);%顾客等待时间累加ww=zeros(1,21);%收款台空闲时间累加r=exprnd(2,1,21);%服从指数分布的随机数 t2=normrnd(1,1/3,1,21);%服从正态分布的随机数for i=1:1:20 t(i+1)=t(i)+t1(i); if t(i+1)=T(i); w(i+1)=w(i); T(i+1)=t(i+1)+t2(i); ww(i+1)=t(i+1)-T(i)+ww(i); else w(i
18、+1)=T(i)-t(i+1)+w(i); T(i+1)=T(i)+t2(i); ww(i+1)=ww(i); end;end;b=t,T,w,ww;%求队伍长brow,bcol=size(b);b=b,zeros(brow,1);b(1,bcol+1)=0;for j=2:brow l=0; if j-l-10 while b(j,1)=b(j-l-1,2) l=l+1; end; b(j,bcol+1)=l;end;end; bg1=w(end)/20 %平均等待时间g2=sum(T-t)/20 %平均逗留时间g3=T(end)/20 %平均每分钟服务的顾客人数 3. 系统仿真软件 MAT
19、LAB-SIMULINK Simulink 是一个用来对动态系统进行建模、仿真和分析的软件包。它支持线性和非线性系统，连续和离散时间模型等。 Smulink 提供了一个图形化的用户界面，可以用鼠标点击和拖拉模块的图标建模。通过图形界面，可以象用铅笔在纸上画图一样画模型图。 Simulink 包括一个复杂的由接受器、信号源、线性和非线性组件及连接件组成的模块库，每个组件是包含若干模块的模块集。当然也可以定制或者创建用户自己的模块。Simulink 可以对模型进行仿真，使用显示模块可以在运行仿真时观察到仿真的结果。还可以改变参数并且立即就可以看到它的变化。仿真结果放在工作空间（workspace）
20、中以待进一步的处理或可视化。 Simulink 使用不同窗口分别显示模块库 simulink library browser, 模型 model space, 和仿真图形的输出work space . 使用命令“simulink” 或工具条中的快捷钮进入模块库浏览器。 在“File”下点击“Model”就可以进入模型空间组建模型。 直接从模块库中将模块拖拉过来就可以拼接成所需要的模型。 它可以以 .M 文件的形式保存下来。 打开 Scope 模块以观察仿真的结果。 选择菜单“Simulation”的“Parameters”以打开“Simulation parameters”窗口设置参数。 选择“simulation | start”菜单开始仿真，并观察Scope模块的输出。习题。2. 进一步建模分析市场服务的问题 10. 假设对不同的顾客收款和装袋的时间服从正态分布 N(1，0.1)。对同一位顾客两者相等。 20. 探讨顾客的到来对出口设置的影响50 以上有不当之处，请大家给与批评指正，以上有不当之处，请大家给与批评指正，谢谢大家！谢谢大家！
