随机模拟(仿真)-simulation

随机模拟（仿真）-simulation仿真（也称为模拟）：就是用计算机程序在计算机上模仿各种实际系统的运行过程，并通过计算了解系统随时间变化的行为或特性。计算机仿真：是在已经建立的数学、逻辑模型之上，通过计算机实验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另外一个状态的行为进行描述和分析。实际问题数学、逻辑模型计算机模型数学、计算机解实际解为什么要进行仿真•实际系统建立之前，要对系统的行为或结果进行分析研究；•有些真实系统做实验会影响系统运行，例如，在生产中任意改变工艺系数可能导致废品，在经济活动中随意将一个决策付诸行动会导致经济混乱；•在系统上做多次试验，很难保证每次试验的操作条件相同，因而对实验结果好坏很难作出正确的判断；•当人是系统的一部分时，他的行为往往实验结果有所影响，这时，最好进行模拟研究；•实验时间太长，费用太大，或者有危险，使得试验不容易进行；•有些系统一旦建立起来后无法复原，例如，建立大型企业，要分析社会和经济效益，不能用建立起来试试看的办法。那些问题适合计算机仿真解决•难以用数学公式表示的系统，或者没有求解数学模型的有效方法；•虽然可以用解析的方法解决问题，但是问题的分析与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法；•希望在较短时间内观察到系统的发展全过程，以估计某参数对系统行为的影响；•难以在实际环境中进行试验和观察，计算机仿真是唯一的方法；•需要对系统或过程进行长期运行的比较，从大量方案中寻找最优。仿真的分类•模拟是系统状态随时间而变化的动态写照，因此，通常时间是模拟的主要自变量，其它的变量为因变量。•（1）按照模拟过程中因变量的变化情况，可以将模拟分为离散、连续、混合3种类型；•（2）如果采用模拟计算机、采用数字计算机以及联合使用则分为模拟仿真、数字仿真以及混合仿真；•（3）根据仿真变量的特征分为随机模拟仿真和模糊模拟仿真。模拟的方法•设计正确的模拟时间推进机理是进行模拟的一个非常重要的问题，模拟过程应该根据系统的特征正确推进模拟时间，使系统中各项要素与发生的事件保持同步，推进时间模拟的基本方法有：•（1）下次事件法：是将模拟时间由一个事件发生时间点推进到紧接着下一个事件发生的时间点。既时间变化幅度由事件变化确定；•（2）固定时间步长法：此种方法模拟时间每次均以相等的固定步长向前推进，每到达一个新的模拟时间点需要检查相应的时间段内是否发生了事件；（有时候需要动态调整步长）模拟的一般步骤•明确问题，建立模型：正确描述研究的问题，明确规定模拟的目标和任务，确定衡量系统性能或模拟输出结果的目标函数，然后根据系统的结构及作业规则，分析系统各状态变量之间的关系，以次为基础建立所研究的系统模型；•收集和整理数据资料：模拟的实现往往离不开大量数据的输入，且需要确定随机因素的概率分布特性，并以此为抽样的根据；•编制程序：模拟运行，选择适当的计算机语言，按照系统数学、逻辑模型编写计算机程序。•分析模拟输出结果：一般包括如下几个方面•（1）模拟结果的统计特性：均值、方差以及置信区间；•（2）灵敏度分析；•（3）根据确定的目标函数，在众多的实现方案中选取最优方案。仿真的基础•概率论的大数定理；•微积分基础；•各种分布：在matlab模拟中，常用的随机分布为unif(均匀分布),exp(指数分布),norm(正态分布),chi2(χ2分布),t(t分布),f(F分布),bino(二项分布),poiss(泊松分布),unid(整数均匀分布)；•有关分布的计算功能：pdf(概率密度),cdf(分布函数),inv(逆概率分布),stat(均值与方差),rnd(随机数生成)；•将分布与计算功能结合在一起，就是实现某种及算功能，例如，expinv就是计算指数分布的逆分布。用随机模拟计算积分例1如下图所示，在正方形内有1/4单位圆。向正方形内投小石头，假设每次都能够投进正方形内且可以落在正方形内任何一点。问，小石头落在1/4单位圆内(包含边界)的概率多大？01x1y分析：假设头入正方形内的石头有n块，有k块落入了1/4单位圆内。P为小石头落入1/4单位圆内的概率。那么根据Bernoull(伯努利)大数定理，有1}|pnk{|plimn即，当实验次数n充分大时，频率和概率之差小于任意数ε的概率趋于1。而另外一方面由几何概型有4Ap样本空间的测度的测度事件在实际操作中，实验次数n不可能趋于无穷大，所以有nk4p44p(数学模型)对于估计nk4只有不断重复做实验，这种试验可以具体去操作，(均匀投石块，然后数数，这样需要较高成本)。也可以让计算机去重复试验，但是需要将数学模型转化为计算机模拟模型(让计算机完成均匀投石块，自动计数，也需要成本)。用计算机模拟投石块过程和步骤如下：1、自动生成随机点[0,1]x[0,1]，模拟石块在正方形内的任意位置，用(xi,yi)表示，共n个点；2、判别(xi,yi)是否满足xi2+yi2≤1，即判别石块是否落在1/4单位圆内，共k个点满足；3、整理、统计模拟结果，用4k/n估计π。利用matlab编程计算过程1、编写M文件：function[pai]=fangzhenpai(n)%生成均匀分布;X=unifrnd(0,1,n,2);%判别是否落在1/4圆内;k=0;fork1=1:ny=X(k1,1)^2+X(k1,2)^2;ify=1k=k+1;elsek=k;endend%整理、统计、估计结果;pai=4*k/n;2、计算过程n=[2,20,200,2000,5000,10000,50000,100000];fork=1:length(n)pai(k)=fangzhenpai(n(k));endpaipai=2.00003.20003.34003.17403.14643.14603.13933.1423值得注意的是，每次模拟即或程序相同，不同人不同计算时间，结果一般不同(因为随机)。例2计算定积分dxx1322分析：若随机变量X的概率分布密度是P(x)(a≤x≤b)，则随机变量Y=f(X)的数学期望为badx)x(p)x(f)]X(f[E)Y(E当X是[a,b]上的均匀分布时，有bxa,ab1)x(p将p(x)代入上面的期望算子，有)]X(f[E)ab(dx)x(fdxab)x(f)]X(f[Ebaba而另外一方面，根据大数定理，设y1,y2,…,yn为来自总体Y的一组样本，有1}|)y(Ey{|limn其中)]X(f[E)Y(E)y(E,yn1yn1ii所以有niin1iiba)x(fnabny)ab(y)ab()]X(f[E)ab()Y(E)ab(dx)x(f其中,x1,x2,…,xn为来自服从[a,b]上均匀分布总体X的一组样本。且yi=f(xi)，i=1,2,…,n。dxx1322那么计算的步骤如下：1、在[2，3]上抽样x1,x2,…,xn；2、计算n,...,2,1i,x1y2ii3、整理估计定积分；n1i2i322x1n23dxx1将上述过程编写M文件如下functionI=jifen1(a,b,n)X=unifrnd(2,3,n,1);L=(b-a)/n;fork=1:nY(k)=(1+X(k)^2)^(1/2);endI=L*sum(Y);a=2;b=3;n=30000;I=jifen1(a,b,n)I=2.6912计算结果如下：而该积分的准确值为：symsxint((1+x^2)^(1/2),2,3)ans=3/2*10^(1/2)-1/2*log(-3+10^(1/2))-5^(1/2)-1/2*log(2+5^(1/2))double(ans)ans=2.6948例3估计事件发生的概率}9XX,3XXPr{243221其中X1服从均匀分布U[2,5]，X2服从指数分布exp(3)，X3和X4分别服从正态分布N(3,2)和N(1,1)。1、在总体X=(X1,X2,X3,X4)中抽样，得到样本(x1k,x2k,x3k,x4k)，k=1,2,…,N；2、判断随机样本是否满足不等式组，假设满足不等式组的有n个；3、用n/N估计θ。步骤：编写M文件：functionp=fangzhenguji(N)x1=unifrnd(2,5,N,1);x2=exprnd(3,N,1);x3=normrnd(3,2,N,1);x4=normrnd(1,1,N,1);n=0;fork=1:Ny1=x1(k)+x2(k)^2;y2=x3(k)+x4(k)^2;ify1=3&y2=9n=n+1;elsen=n;endEndp=n/N;计算结果p=fangzhenguji(10000)p=0.8415例4求使下式成立的最大f值8.0}fXXX{pr33221其中，X1服从均匀分布U[1,3]，X2服从指数分布exp(1)，X3服从正态分布N(2,1)。步骤：1、产生N个独立随机序列(X1k,X2k,X3k),k=1,2,…,N；2、计算fk=f(X1k,X2k,X3k),k=1,2,…,N;3、从序列f1,f2,…,fN中取[0.8N]作为f的估计值。编写文件：function[fbar]=jieji(m,a)x1=unifrnd(1,3,m,1);x2=exprnd(1,m,1);x3=normrnd(2,1,m,1);fork=1:mf(k)=x1(k)+x2(k)^2+x3(k)^3;endn=floor(a*m);f1=[];fork=1:mb=max(f);bb=find(fb-0.00001);rr=bb(1);f1=[f1;f(rr)];f(rr)=[];endbar=f1(n);fbar=jieji(10000,0.8)fbar=4.9698计算结果：例5一个排队系统的仿真某店只有一个收款台，顾客到达收款台的时间间隔服从均值为4.5的负指数分布，每个顾客的服务时间服从均值为3.2，标准差为0.6的正态分布，对100位顾客去收款台缴款的排队过程进行仿真。分析：排队系统涉及的参考变量为：顾客的等待时间Wq，顾客的逗留时间W,排队长Lq,队长L，服务台的空闲时间I，服务台的繁忙时间B，排队系统的状态概率P。最后需要统计出有关指标的平均值均方差。合理假设1、进入系统的顾客没有人因为不愿意等而离开系统；2、先到先服务；仿真步骤：1、产生初始事件表，并从初始事件表中找出最靠前的事件；2、启动仿真时钟；3、判别当前事件是哪一个事件，如果是顾客到达事件，则启动顾客到达子程序；如果是服务事件，则启动服务结束子程序；4、重复2-3，直到仿真完毕；5、整理、分析并输出结果。顾客到达的子程序1、产生下一个顾客到来的时刻，记入事件表；2、判别服务台是否忙s=1?，若是，则顾客则排队长L=L+1;否则，服务台空闲,L=0;3、产生一个服务结束时刻，记入事件表；4、统计所需的数据；服务结束子程序1、结束服务，产生结束时间，记入事件表；2、判别是否有顾客排队，L=0?，若L=0，则服务台空闲；否则，排队长L=L-1;3、产生服务结束时刻，记入事件表；4、统计所需的数据；产生初始事件表从事件表中找出最靠前的事件仿真时钟步进是哪一类事件顾客到达服务结束产生下一个顾客到达纪录，记入事件表结束服务，产生结束时间，记入事件表出纳是否忙？出纳空闲，s=0顾客排队，Lq=Lq+1否是产生服务结束时刻，记入事件表是否顾客排队队长减1，Lq-1出纳空闲，s=0否产生服务结束时刻，记入事件表是统计所需数据仿真是否完毕是结束程序否M文件程序:paiduisimu.mfunction[ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu]=paiduisimu(n,a,b,c)%L时钟进步步长；n需要处理的人数；a表示负指数分布参数；b,c表示正态分布参数;t=0;%起步时钟;LL=[];%