1、Chap 8 离散事件系统仿真离散事件系统仿真前面讨论的系统，其状态变量的取值是连续变化的（时间上可以连续也可以离散），这类系统的仿真称为连续系统仿真。现开始讨论另一类性质完全不同的系统，其状态只是在离散时间点上发生变化，且这些离散时间点一般是不确定的，称为离散事件系统离散事件系统仿真仿真。例如单人理发馆系统，设上午9点开门，晚上11点关门，顾客的到达时间一般是随机的，为每个顾客服务的时间长度也是随机的。描述该系统的状态是服务台的状态（忙或闲）、顾客排队等待的队长。显然这些状态变量的变化也只能在离散的随机时间点上发生。类似的如订票系统、库存系统、加工制造系统、交通控制系统、计算机系统等。由于离
2、散事件系统固有的随机性，对这类系统的研究往往十分困难，经典的概率及数理统计理论、随机过程理论虽然为研究这类系统提供了理论基础，并能对一些简单系统提供解析解，但对工程实际中的大量系统，惟有依靠计算机仿真才能提供较为完整的结果。8.1 基本概念基本概念1、实体实体 实体是描述系统的三个基本要素之一，在离散事件系统中的实体可分为两大类：临时实体及永久实体。在系统中只存在一段时间的实体叫临时实体。永久驻留在系统中的实体称为永久实体。临时实体按一定规律不断地到达（产生），在永久实体作用下通过系统，最后离开系统，整个系统呈现出动态过程。2、事件事件 事件是引起系统状态发生变化的行为。从某种意义上讲，离散系
3、统是由事件来驱动的。如，理发馆系统中，可以定义“顾客到达”为一类事件，由于顾客的到达，系统的状态将发生变化服务员可能从闲变忙（如果无人排队），或排队的队长会增加。类似的，可以定义服务开始事件、服务结束事件。在一个系统中，往往有许多类事件，而事件的发生一般与某一类实体相联系，有些事件的发生还可能引起别的事件的发生，或是另一类事件发生的条件。为了实现对系统中事件进行管理，仿真模型中必须建立事件表，表中记录每一发生了的或将要发生的事件的类型和发生时间，以及与该事件相联的实体的有关属性等等。3、活动活动 离散事件系统中的活动，通常用于表示两个可以区分的事件之间的过程，它标志着系统状态的转移。如前例中，
4、顾客的到达事件与该顾客开始接受服务事件之间可以称为一个活动（排队活动），该活动使系统的状态（队长）发生变化。4、进程进程 进程是由若干个有序的事件及若干有序活动组成，一个进程描述了它所包括的事件及活动间的相互逻辑关系及时序关系。如顾客到达、经过排队、接受服务、服务完毕后离去可称为一个进程。进程排队活动服务活动顾客到达事件服务开始事件服务结束事件5、仿真钟仿真钟 仿真钟用于表示仿真时间的变化。在连续系统仿真中，将连续模型进行离散化而成为仿真模型时，仿真时间的变化基于仿真步长的确定，可以是定步长也可以是变步长，称为时间步长法。对于离散事件系统而言，其状态本来就只在离散时间点上发生变化，因而不需要进
5、行离散化处理。但是由于引起状态变化的事件发生时间的随机性，仿真钟的推进步长则完全是随机的，所以说仿真模型中时间控制部件必不可少，应按一定规律来控制仿真钟的推进。8、统计计数器统计计数器 连续系统仿真的目的是要得到状态变量的动态变化过程并由此分析系统的性能。离散事件系统的状态随事件的不断发生也呈现出动态变化的过程，但仿真的目的主要不是要得到这些状态是如何变化的。由于状态的变化是随机的，某一次仿真运行得到的状态变化过程只不过是随机过程的一次取样。如果进行另一次独立的仿真运行所得到的状态变化过程可能完全是另一种情况。所以它们只有在统计意义下才有参考价值。在前例中，由于顾客到达的时间间隔具有随机性，服
6、务员为每个顾客服务的时间长度也是随机的，因而在某一时刻，顾客排队的队长或服务台的忙闲情况完全是不确定的，在分析该系统时，感兴趣的是系统的平均队长、顾客的平均等待时间或服务员的利用率等。所以在仿真模型中，需要有一个统计计数部件，以便统计系统中的有关变量。8.2 仿真钟的推进仿真钟的推进 离散事件系统仿真的仿真钟推进方法有两种：一种是按下一最早发生事件的发生时间推进，称为事件调度法，亦称为事件步长法；另一种是固定增量推进法。事件步长法事件步长法事件步长法是以事件发生的时间为增量，按照事件发生的时间顺序，一步一步地对系统的行为进行仿真，直到预定的时间结束为止。在多数随机系统中，可以有多种性质的事件发
7、生，通常按照发生时间的先后顺序逐个处理，换句话说，首先处理发生时刻距初始时刻最短的事件，这种处理方法称为是最短时间的事件步长法。事件步长法事件步长法初始状态事件步长加1在当前步长内，考察分析，计算和记录系统的活动仿真时间到否？仿真结束输出结果是否事件步长法事件步长法事件步长法与时间步长法的主要区别是：（1）事件步长法与时间步长法都是以时间为增量来考察系统状态的变化但在时间步长法中，仿真时钟以等步长前进，而在事件步长法中，仿真时钟的步长取决于事件之间的间隔（2）时间步长法在一个步长内，认为系统所处的状态相同，因而所选步长的大小将影响仿真的精度而在事件步长法中，每个事件的发生均有确切的时刻，不需要
8、人为地选取步长，步长的大小对仿真精度影响较小事件步长法事件步长法（3）时间步长法每前进一个步长就要对整个系统进行一次全面考察，即使状态没有发生变化时也要扫描，而事件步长法只是在某一事件点上判断和比较事件是否出现因此，一般地讲，当判断比较的数目较大时，用时间步长法可以节省用机时间，而当相继两个事件出现的平均间隔较长时，更适合于用事件步长法事件步长法事件步长法事件表法的基本思路事件表法的主要思路是将系统的仿真过程看成一个事事件点序列件点序列，根据事件出现的时序，用一个称之为事件表的表格来调度事件执行的顺序对于那些当前需处理的事件，列入事件表中，从中取出最接近的事件进行处理，处理完毕后自动退出事件表
9、在处理当前事件的过程中，往往又会产生一个后继事件，因此，必须预测出这一后继事件的出现时间，并将它列入事件表中事件步长法事件步长法-例子例子例例1收款台前的排队过程的仿真。收款台前的排队过程的仿真。考虑一个收款台的排队系统。某个杂货店只有一个收款台,顾客的到达时间时服从均值为4.5的负指数分布，每一个顾客的服务时间服从均值为3.2，方差为0.8的正态分布。这里时间的单位是分钟,且服务的时间不取负值.试对收款台前的排队过程进行仿真。事件步长法事件步长法-例子例子收款排队系统主控程序图事件步长法事件步长法-例子例子负指数分布的随机数的产生：负指数分布的随机数的产生：x=-4.5ln(u),u为均匀分
10、布的随机数正态分布的随机数的产生正态分布的随机数的产生:y=3.2+0.8v,v为标准正态分布的随机数。假设假设：开始时服务台前无顾客，顾客在服务台前不离去；要求要求：对100个顾客到收款台缴款排队过程进行仿真关心的问题关心的问题：每个顾客的平均等待时间atime,最大队长maxl,服务员的工作效率work,；事件步长法事件步长法-例子例子实体实体：服务员（忙1，闲0），顾客(到达时刻ca，服务时间cs，等待 时间ct)队列（长度lq）；事件事件event：顾客到达事件1（交款，排队）服务结束事件2（闲忙）事件表事件表：序号 事件类型（1、2）发生时刻t事件步长法事件步长法-例子例子顾客到达子
11、程序图事件步长法事件步长法-例子例子服务结束子程序图事件步长法事件步长法-例子例子例:机器修理系统仿真有m台机器，由c个工人共同负责修理，并设：各台机器质量相同，机器的连续运转时间相互独立且服从同一负指数分布，平均寿命为1/v（v0）。每个工人技术相同,且修理时间相互独立并服从同一负指数分布,平均修理时间为1/u(u0).修复后的机器其寿命分布不变。机器停止运转每单位时间的损失费为c1元，工人单位时间的产值为c2元。若机器的等待时间为E E，工人总的空闲时间为F F，则系统总的损失费为 S S=c1E E+c2 2F F试求当机器数m固定时，为使系统的总损失费最小，应配备多少工人为最优？事件步
12、长法事件步长法-例子例子解:采用最短时间事件步长法进行仿真，首先要确定一个初始状态。不妨假定开始时所有的机器都正常运转，工人处于空闲状态。设T为仿真终止时间。依次仿真计算修理工人数C分别取1、2、3、时的系统损失费。最后根据系统损失费S的极小值来确定最优工人数。根据以上思想并采用最短时间事件步长法仿真。事事件件步步长长法法-例例子子输入原始数据c1,c2,U,V,m,T给出m台机器的指数寿命B(H)=1?结束F=F+x-KP=T3/T输出C,SC=1E=0 F=0初始状态:B(i)=1,i=1,2,m;D(j)=0,j=1,2,.,c选取寿命最短的机器i:H=i,K=A(i)选取最早释放时刻的
13、工人j:y=j,x=D(j)xK?F=F+x-K改变第y个人工的释放时刻:D(y)=A(H)机器由故障释放改为正常:B(H)=1机器由正常改为故障:B(H)=0C=C+1KT?选取寿命最短的机器i:H=i,K=A(i)给出机器修复时间仿真时钟前进一个步长A(H)=K-1/vLOG(RND(1)给出机器正常运行时间仿真时钟前进一个步长A(H)=K-1/vLOG(RND(1)是否否是是否事件步长法事件步长法-例子例子框图中各标识符号含义如下：T:仿真终止时间m:机器数c:工人数B(i):第i台机器所处状态A(i):第i台机器现在状态的当前时刻D(j):第j个工人修复机器的释放时刻c1:一个工人单位
14、时间的产值c2:一台机器单位故障时间的损失费B(i)=1 机器正常运转0 机器发生故障事件步长法事件步长法-例子例子U:一个工人单位时间平均能修复的机器数V:一台机器在单位运转时间内发生故障的平均次数S:总损失费用E:机器故障时间累计F:工人空闲时间累计K:最短时间机器的当前时刻H:最短时间机器对应的机器号L:工人修完机器的释放时刻x:最短释放时间工人的当前时刻y:最短释放时间对应的工人序号事件步长法事件步长法-例子例子下表列出了当 m=86，1/V=500小时，1/U=34小时，c1=3.46元/小时，c23.2元/小时时的仿真结果，其中每次预定仿真时间为一万个小时，连续仿真五次取其平均值作
15、为仿真结果。由仿真结果可知当工人数c10时，最优工人数为7工人数C2345878910每小时平均损失E 183.0 128.9 76.0 37.624.59.09.410.2 15.2应用举例应用举例-可靠性问题例：一设备上三个相同的轴承，每个轴承正常工作寿命为随机变量，其概率分布如表所示寿命h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05任何一个轴承损坏都可以使设备停止工作，从有轴承损坏，设备停止工作，到检修工到达开始更换部件为止，称为一个延迟时间延迟时间也是随机变量，其概率
16、分布如下表所示延迟时间min51015概率0.60.30.1应用举例应用举例-可靠性问题设备停工时每分钟损失5元,检修工每小时工时费12元，轴承每个成本 16元更换一个轴承需要 20 min，同时更换两个轴承需要30min，同时更换三个轴承需要40min现在有两种方案：方案一是损坏一个更换一个；方案二是一旦有轴承损坏就全部更换试通过计算机仿真对这两种方案做出评价在这一问题中，轴承寿命在1000到1900h之间，而延迟在5到 15min之间，故若用时间步长法时，步长选取有些困难步长小浪费很大，步长大又不精确，所以采用事件步长法在事件发生时再考虑系统状态的变化情况，这就比较合理应用举例应用举例-可
17、靠性问题为了进行仿真，首先对轴承寿命和延迟时间与随机数对应，对应规则分别如下两个表轴承寿命h频率随机数区间10000.10(0,0.10)11000.130.10,0.23)12000.250.23,0.48)13000.130.48,0.61)14000.090.61,0.70)15000.120.70,0.82)16000.020.82,0.84)17000.060.84,0.90)18000.050.90,0.95)19000.050.95,1.00)延迟时间频率随机数区间50.60(0,0.6)100.300.6,0.9)150.100.9,1.0)轴承寿命随机数延迟时间随机数应用举例
18、应用举例-可靠性问题由于在这一问题中各个轴承的寿命完全决定了系统的运行状态，也即决定了两个方案的费用大小，故我们选择轴承发生故障作为事件，这三个轴承发生故障的事件分别记为A、B、C（1）方案一的仿真 产生初始事件表事件类型发生时刻延迟时间A1400h5minB1500h15minC1500h15min应用举例应用举例-可靠性问题仿真时钟步进，计算费用，产生下一个事件由表上看出，最早发生的事件是A，所以t 1400hcost(5+20)54116145元，下一个A事件发生的时刻为第2 400小时25分钟（随机产生轴承寿命为1000h），刷新事件表，即删去老的A事件，产生新的A事件刷新后的事件表如
19、下表事件类型发生时刻延迟时间A2400h25min5minB1500h15minC1500h15min应用举例应用举例-可靠性问题寻找事件表中的最早事件进行处理由上表看出，B、C事件同时发生在第1500 h，故同时处理时钟步进为 t 1500 h，再根据费用的计算方法得：cost=145(15+30)5 6216408元 最后利用随机数产生新事件B和C，刷新事件表，得到的新事件表:事件类型发生时刻延迟时间A2400h25min5minB2700h45min10minC2900h45min5min应用举例应用举例-可靠性问题重复，t=2400 h 25 min，cost=408(5+20)5 4
20、 116=553元事件类型发生时刻延迟时间A3700h50min5minB2700h45min10minC2900h45min5min重复这一过程，一直需要的时间结束即可得到方案的费用应用举例应用举例-可靠性问题（2）方案二的仿真方案二与方案一的区别就是一旦故障发生，就更换3个轴承设初始事件表仍为方案一初始表表中最早的事件是A，处理事件A时要考虑延时，更换3个轴承的时间和费用事件类型发生时刻延迟时间A1400h5minB1500h15minC1500h15min cost=(5+40)58316281元应用举例应用举例-可靠性问题根据下一次三个轴承发生故障的时刻刷新后的事件表:事件类型发生时刻
21、延迟时间A2400h45min5minB2600h45min10minC2800h45min5min再重复上述过程，累加费用，即可得到方案二的总费用最后比较两种方案的费用大小即可确定选取那一种应用举例应用举例-可靠性问题程序运行结果：T=100000h方案一：cost=32705元方案二：cost=24429元从而得出方案二较方案一优。应应用用举举例例-可靠性问题/可靠性问题#include#include#include#define TIMES 100000 struct table float begin;int delay;int make;struct table tab3;floa
22、t cost=0;int sort3;int make_event1(void)int first=0,sum;sort0=sort1=sort2=-1;float b=tab0.begin;if(btab1.begin)first=1,b=tab1.begin;if(btab2.begin)first=2,b=tab2.begin;int i=0;应应用用举举例例-可靠性问题 if(first=0)sorti+=0;if(tab1.begin=tab0.begin)sorti+=1;if(tab2.begin=tab0.begin)sorti+=2;else if(first=1)sorti
23、+=1;if(tab2.begin=tab1.begin)sorti+=2;else if(first=2)sorti+=2;for(sum=0,i=0;i=0)sum+;switch(sum)case 1:tabsort0.make=20;cost+=(20+tabsort0.delay)*5+4+16;break;case 2:tabsort0.make=30;cost+=(30+tabsort0.delay)*5+6+32;break;case 3:tabsort0.make=40;cost+=(40+tabsort0.delay)*5+8+48;break;return sum;应应用
24、用举举例例-可靠性问题int make_event2(void)int first=0,sum;sort0=sort1=sort2=-1;float b=tab0.begin;if(btab1.begin)first=1,b=tab1.begin;if(btab2.begin)first=2,b=tab2.begin;int i=0;if(first=0)sorti+=0;if(tab1.begin=tab0.begin)sorti+=1;if(tab2.begin=tab0.begin)sorti+=2;else if(first=1)sorti+=1;if(tab2.begin=tab1.
25、begin)sorti+=2;else if(first=2)sorti+=2;tabsort0.make=40;cost+=(40+tabsort0.delay)*5+8+48;return sum;应应用用举举例例-可靠性问题 void create_event1(struct table*t)float x,y;struct table*t1;x=random(1000);x/=1000;t-begin+=(t-delay+t-make)/80.0;if(xbegin+=1000;else if(xbegin+=1200;else if(xbegin+=1300;else if(xbeg
26、in+=1400;else if(xbegin+=1500;else if(xbegin+=1800;else if(xbegin+=1700;else if(xbegin+=1800;else if(xbegin+=1900;for(t1=&tab0;t1begin=t-begin)t-delay=t1-delay;return;y=random(1000);y/=1000;if(ydelay=5;else if(ydelay=10;else if(ydelay=15;return;应应用用举举例例-可靠性问题 void create_event2(struct table*t,struct
27、 table*t0)float x,y;struct table*t1;x=random(1000);x/=1000;if(t0!=0)t-begin=t0-begin+(t0-delay+t0-make)/80.0;else t-begin+=(t-delay+t-make)/80.0;if(xbegin+=1000;else if(xbegin+=1200;else if(xbegin+=1300;else if(xbegin+=1400;else if(xbegin+=1500;else if(xbegin+=1800;else if(xbegin+=1700;else if(xbegi
28、n+=1800;else if(xbegin+=1900;for(t1=&tab0;t1begin=t-begin)t-delay=t1-delay;return;y=random(1000);y/=1000;if(ydelay=5;else if(ydelay=10;else if(ydelay=15;return;应应用用举举例例-可靠性问题 void printf_result(void)for(int i=0;i3;i+)printf(%dt%ft%dn,i+1,tabi.begin,tabi.delay);printf(cost=%fn,cost);return;void main1
29、(void)int i;float T=0;memset(tab,0,sizeof(tab);cost=0;randomize();create_event1(&tab0);create_event1(&tab1);create_event1(&tab2);printf_result();domake_event1();T=tabsort0.begin;printf_result();for(i=0;iTIMES);printf(1:tCost=%f T=%fn,cost,T);应应用用举举例例-可靠性问题void main2(void)int i;float T=0;memset(tab,0
30、,sizeof(tab);cost=0;randomize();create_event2(&tab0,0);create_event2(&tab1,0);create_event2(&tab2,0);printf_result();domake_event2();T=tabsort0.begin;printf_result();for(i=0;iTIMES);printf(2:tCost=%f T=%fn,cost,T);void main(void)main1();main2();固定增量法固定增量法 选择适当的时间单位T作为仿真钟推进进的增量，每推进一步进行如下处理：1）该步内若无事件发
31、生，则仿真钟再推进一个时间单位；2）若在该步内有若干个事件发生，则认为这些事件均发生在该步的结束时刻。缺点是：仿真钟每推进一步都要检查事件表以确定是否有事件发生，增加了执行时间；将发生的事件均视为发生在这一步的结束时刻，如果T选得过大，会引入较大的误差；且要求确定各类事件处理的顺序，增加了建模的复杂性。主要用于系统事件发生时间具有较强周期性的模型，如定期订货的库存系统，以年、月为单位的经济计划系统等。应用举例应用举例-报童的策略例例 报童每天清晨从报社购进报纸零售,晚上将没有卖掉的报纸退回.每份报纸的购进价为1.3元,零售价为2元,退回价为0.2元.报童售出一份报纸赚0.7元,退回一份报纸赔1
32、.1元.报童每天如果购进的报纸太少,不够卖时会少赚钱,如果购得太多卖不完时要赔钱.试为报童筹划每天应如何确定购进的报纸数使得收益最大.报纸每捆10张，只能整捆购买，报纸可以分为3种类型的新闻日：好、一般、差，它们的概率分别为0.35,0.45和0.2,在这些新闻日中每天对报纸的需求分布的统计结果下图：应用举例应用举例-报童的策略需求量好新闻的需求概率一般新闻的需求概率 差新闻的需求概率 40 0.03 0.10 0.44 50 0.05 0.18 0.22 60 0.15 0.40 0.16 70 0.20 0.20 0.12 80 0.35 0.08 0.06 90 0.15 0.04 0.
33、00 100 0.07 0.00 0.00试确定每天报童应该订购的报纸数量应用举例应用举例-报童的策略解：我们通过计算机仿真来解决此问题。最优策略应该是每天的利润最大。利润=销售收入-报纸成本-损失+残值这是一个随机现象的计算机仿真问题，故先确定各种情况的随机数的对应关系。新闻日和需求量对应的随机数分别如下面两个表格所示新闻种类 出现概率对应的随机数区间好新闻 0.35 (0.00,0.35)一般新闻 0.45 0.35,0.80)差新闻 0.20 0.80,1.00)应用举例应用举例-报童的策略计算机仿真的流程：1）令每天的报纸订购数变化，40-100；2）让时间从1开始变化（循环）到365
34、；3）产生新闻种类的随机数，确定当天的新闻类型；4）产生需求量随机数，确定当天的报纸需求量；5）计算当天的收入，计算累积利润，8）比较得出最优定货量。需求量好新闻的随机数区间一般新闻的随机数区间差新闻的随机数区间 40 （0.00,0.03 （0.00,0.10 （0.00,0.44 50 0.03,0.08)0.10,0.28)0.44,0.66)60 0.08,0.23)0.28,0.68)0.66,0.82)70 0.23,0.43)0.68,0.88)0.82,0.94)80 0.43,0.78)0.880.96)0.94,1.00)90 0.78,0.93)0.96,1.00)100
35、 0.93,1.00应用举例应用举例-报童的策略具体的计算由具体的计算由Matlab编程计算实现。编程计算实现。x1=rand(365,1);x2=rand(365,1);for n=4:10 paper=n*10;购买的报纸量 lr(n)=0;利润 for i=1:365 if x1(i)0.35 if x2(i)0.03 news=40;elseif x2(i)0.08 news=50;elseif x2(i)0.23 news=60;elseif x2(i)0.43 news=70;elseif x2(i)0.78 news=80;elseif x2(i)0.93 news=90;els
36、e news=100;endelseif x1(i)0.8 if x2(i)0.10 news=40;elseif x2(i)0.28 news=50;elseif x2(i)0.88 news=60;elseif x2(i)0.88 news=70;elseif x2(i)0.98 news=80;应用举例应用举例-报童的策略 else news=90;end else if x2(i)0.44 news=40;elseif x2(i)0.88 news=50;elseif x2(i)0.82 news=60;elseif x2(i)=news sale=news;实际销售量 remand=
37、paper-news;剩余量 应用举例应用举例-报童的策略 else sale=paper;remand=0;end lr(n)=lr(n)+2*sale-1.3*paper+0.2*remand;endendoptnews=40;最佳量optmoney=lr(4);最佳利润 40,lr(4)/365for n=5:10 if lr(n)=optmoney optnews=n*10;optmoney=sb(n);end n,lr(n)/365endoptnews,optmoney,optmoney/365 Matlab程序 经过计算机仿真后得到最优购货量是每天60份，平均每天利润34.4元。8
38、.3 随机数的产生随机数的产生 现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响要使我们的数学模型能够较真实地刻画实际对象，必须面对这个现实概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具但是它还难以用来处理复杂系统中的随机性这里介绍使用计算机来模拟随机现象的方法，它经常应用于复杂系统的动态仿真的研究当中是处理复杂系统中随机性的计算机模型也是使用计算机研究和解决实际问题的一条重要途径随机现象的模拟随机现象的模拟1均匀分布的随机数及其产生对随机现象进行模拟，实质上是要给出随机变量的模拟也就是说利用计算机随机地产生一系列数值，它们的出现服从一定的概率分布,则称这些数
39、值为随机数。最常用的是在(0,1)区间内均匀分布的随机数，也就是我们得到的这组数值可以看作是(0,l)区间内均匀分布的随机变量的一组独立的样本值以后我们将指出其它分布的随机数可利用均匀分布的随机数产生随机现象的模拟随机现象的模拟记 XU(0,1),0,1)区间上均匀分布的随机变量X的概率密度f(x)和概率分布函数F(x)分别为:111000)(0101)(xxxxxFxxf其它数学期望:E(X)=1/2;方 差:=1/12 随机现象的模拟随机现象的模拟下图为均匀分布的密度函数和分布函数示意图：随机现象的模拟随机现象的模拟均匀随机数是产生其它随机数的基础。例如，抛硬币、抽签、统计经验分布都可以由
40、它产生。产生随机数的方法：（1）硬件设备硬件设备：从真实物理现象的随机因素中产生随机数，放射性粒子的放射源，电子晶体管的固有噪音等，单位时间内放射出的粒子数是随机的。优点：真正的随机数；缺点：外部设备，无法重复随机现象的模拟随机现象的模拟（2）数学公式数学公式：产生伪随机数 用数学公式或位移寄存器的移位操作来产生的随机数,称为伪随机数.因为真实的随机数,只能从客观真实的随机现象本身产生出来,所以产生理想的伪随机数列不是一件容易的事.一般对于产生伪随机数的方法,有如下几点要求:1)要求伪随机数列有较理想的随机性和均匀性,就是对其随机性和均匀性进行统计检验时,有合乎要求的精度;2)产生伪随机数的程
41、序应当简短、运算速度快,占用计算机的内存单元少;3)伪随机数列的循环周期应当尽可能地大,以满足模拟的需要4)伪随机数列中,前后之间和各子列之间,要求相互是独立无关的。随机现象的模拟随机现象的模拟当我们需要一系列均匀分布的随机数时,可以按照一定的算法由计算机计算产生的伪随机数产生随机数的方法很多，其中以乘同余法使用较广用以产生均匀分布随机数的乘同余法的递推公式为：xn+1=xn(mod M),rn=xn/M其中，是乘因子，M是模数前面式子的右端称为以M为模数（modulus）的同余式。给定了一个初值x0(称为种子)后，计算出的rn即为(0,l)上均匀分布的随机数随机现象的模拟随机现象的模拟无论用
42、哪一种方法产生的随机数都存在这样的问题，必须对它进行统计检验看看它们是否具有较好的独立性和均匀性一般在计算机(或计算器)及其使用的算法语言中都有随机数生成的命令，它们所生成的随机数都是经过检验并且可用的.这里就不再详细介绍检验的方法了随机现象的模拟随机现象的模拟C语言中与随机数有关的函数：randomize(void);初始化x=random(int M);产生0M 之间的随机数x=rand(void);产生0215-1之间的随机数Matlab中与随机数有关的函数：Rand(n)n*n个0,1之间均匀分布随机数Rand(m,n)m*n 个0,1之间均匀分布随机数randn(n)n*n个N(0,
43、1)标准正态分布随机数randn(m,n)m*n个N(0,1)标准正态分布随机数随机现象的模拟随机现象的模拟2随机变量的模拟利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟（l）离散型随机变量的模拟设随机变量X的分布律为Pr(X=xi)=pi,i=1,2,令p(0)=0,p(n)=pi,n=1,2,将p(n)作为分点，把区间(0,1)分为一系列小区间(p(n-1),p(n)对于均匀的随机变量RU(0,l)，则有Pr(p(n-1)R p(n)=p(n)-p(n-1)=pi,n=1,2,ni=1随机现象的模拟随机现象的模拟由此可知，事件(p(n-1)R
44、 p(n)和事件(X=xn)有相同的发生的概率因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是：每产生一个(0,1)上均匀分布的随机数（以后简称随机数）r，若p(n-1)r p(n)则理解为发生事件“X=xn”于是就可以模拟随机变量的取值情况随机现象的模拟随机现象的模拟（2）连续型随机变量的模拟一般说来，具有给定分布的连续型随机变量可以利用在区间（0，1）上均匀分布的随机数来模拟最常用的方法是反函数法由概率论的理论可以证明，若随机变量Y有连续的分布函数F(y)，而X是区间（0，l）上均匀分布的随机变量，令ZF-1(X)，则Z与Y有相同的分布由此，若已知
45、Y的概率密度为f(y)，由Y=F-1(X)可得XF(Y)=f(y)dyY-随机现象的模拟随机现象的模拟是区间（0，l）上均匀分布的随机变量如果给定区间（0，1）上均匀分布的随机数ri，则具有给定分布的随机数yi可由方程中解出ri F(Y)=f(y)dyyi0例当需要模拟服从参数为的指数分布时,由ri e-xdy=1-e-yiyi0可得yi=-(1/)ln(1-ri).因为(1-ri)和ri 同为(0，1)区间上的均匀分布的随机数，故上式可以简化为yi=-(lnri)/.随机现象的模拟随机现象的模拟反函数法是一种普通的方法，但当反函数不存在或难以求出时，就不宜于使用了舍选法是另一种方法，其实质是
46、从许多均匀随机数中选出一部分，使之成为具有给定分布的随机数步骤如下：设随机变量X有密度f(x)，又存在实数ab，使 Pr(aXb)=1.1).选取常数，使f(x)l，x(a,b);2).产生均匀的随机数 r1和 r2，令 y=a(b-a)r1；随机现象的模拟随机现象的模拟3).若r2 f(y)，则令x=y，否则剔除r1和r2，重返步骤2，如此重复循环，产生的随机数x1,x1,.,xn,服从的概率分布由密度函数f(x)确定若不存在上述的有限区间，可以选择一个有限区间(a1,b1)使得 Pr(a1X1-其中 是充分小的正数重复上面的步骤，所产生的随机数仅会出现较小的误差随机现象的模拟随机现象的模拟
47、（3）正态随机数的模拟,除了可用反函数和舍选法模拟正态随机变量外,还有两种常用的方法坐标变换法:设r1,r2是相互独立的均匀随机数,令 x1(-2lnr1)1/2cos(2r2)x2=(-2lnr1)1/2sin(2r2)则x1,x2是相互独立的标准正态的随机数随机现象的模拟随机现象的模拟另一个方法是利用中心极限定理设xi,i=1,2,n 是n个相互独立的(0,1)上的均匀的随机变量，有E(xi)=1/2,D(xi)=1/12，由中心极限定理知y=(xi-)将渐近服从正态分布N(0,l)因此，取n个均匀的随机数ri,则有 y=(ri-)n/12li=1nn2n/12li=1nn2随机现象的模拟
48、随机现象的模拟将可看成是标准正态分布N(0,1)的随机数这里的n应取得足够大，一般取n=10即可若取n=12,则上式简化为 z=ri-8再由y=x+得到正态分布N(,2)的随机数y.i=1n计算机仿真举例计算机仿真举例出发时刻（T）1：00 1：051：10 频率 0.7 0.2 0.1例 赶火车过程仿真一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站运行时间为均值30分钟、标准差为2分钟的正态随机变量.火车大约在下午1点离开A站。离开时刻的频率分布为计算机仿真举例计算机仿真举例 0.1 0.2 0.4 0.3 频率1：34 1:32 1:30 1:28到达时刻
49、（T）用计算机仿真火车开出、火车到达B站、这个人到达B站情况，并给出能赶上火车的仿真结果。引入以下变量：T1 火车从A站开出的时刻；T2 火车从A站运行到B站所需要的时间；T3 此人到达B站的时刻；0.1 0.2 0.7 P（频率）10 5 0 T1(分)这个人到达B站时的频率分布为：计算机仿真举例计算机仿真举例 T1,T2,T3 是随机变量，其概率分布为 x1=0.7,x2=0.9,y1=0.3,y2=0.7,y3=0.9,开车时间的仿真测试 s1=0;s2=0;s3=0;求概率0.7 0.2 0.1 x=rand(10000,1);for i=1:10000 if x(i)0.9 s3=s
50、3+1;end ends1/10000,1-s1/10000-s3/10000,s3/10000计算机仿真举例计算机仿真举例 T2（分）28 30 32 34 P(频率)0.3 0.4 0.2 0.1 s1=0;s2=0;s3=0;s4=0;x=rand(10000,1);for i=1:10000 if x(i)0.3 s1=s1+1;elseif x(i)0.7 s2=s2+1;else if x(i)0.9 s3=s3+1;else s4=s4+1;end endends1/10000,s2/10000,s3/10000,s4/10000人到达时刻仿真测试人到达时刻仿真测试计算机仿真举例
