SIMULINK仿真基础之离散时间系统分析.ppt

《SIMULINK仿真基础之离散时间系统分析.ppt》由会员分享，可在线阅读，更多相关《SIMULINK仿真基础之离散时间系统分析.ppt（70页珍藏版）》请在三一办公上搜索。

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为均匀分布的随机数正态分布的随机数的产生:y=3.2+0.8v,v为标准正态分布的随机数。假设：开始时服务台前无顾客，顾客在服务台前不离去；要求：对100个顾客到收款台

10、缴款排队过程进行仿真,关心的问题：每个顾客的平均等待时间atime,最大队长maxl,服务员的工作效率work,；,事件步长法-例子,实体：服务员（忙1，闲0），顾客(到达时刻ca，服务时间cs，等待 时间ct)队列（长度lq）；事件event：顾客到达事件1（交款，排队）服务结束事件2（闲忙）事件表：序号 事件类型（1、2）发生时刻t,事件步长法-例子,顾客到达子程序图,事件步长法-例子,服务结束子程序图,事件步长法-例子,例:机器修理系统仿真有m台机器，由c个工人共同负责修理，并设：各台机器质量相同，机器的连续运转时间相互独立且服从同一负指数分布，平均寿命为1/v（v0）。每个工人技术相同

11、,且修理时间相互独立并服从同一负指数分布,平均修理时间为1/u(u0).修复后的机器其寿命分布不变。机器停止运转每单位时间的损失费为c1元，工人单位时间的产值为c2元。若机器的等待时间为E，工人总的空闲时间为F，则系统总的损失费为 S=c1E+c2F试求当机器数m固定时，为使系统的总损失费最小，应配备多少工人为最优？,事件步长法-例子,解:采用最短时间事件步长法进行仿真，首先要确定一个初始状态。不妨假定开始时所有的机器都正常运转，工人处于空闲状态。设T为仿真终止时间。依次仿真计算修理工人数C分别取1、2、3、时的系统损失费。最后根据系统损失费S的极小值来确定最优工人数。根据以上思想并采用最短时

12、间事件步长法仿真。,事件步长法-例子,事件步长法-例子,框图中各标识符号含义如下：T:仿真终止时间m:机器数c:工人数B(i):第i台机器所处状态A(i):第i台机器现在状态的当前时刻D(j):第j个工人修复机器的释放时刻c1:一个工人单位时间的产值c2:一台机器单位故障时间的损失费,B(i)=,1 机器正常运转,0 机器发生故障,事件步长法-例子,U:一个工人单位时间平均能修复的机器数V:一台机器在单位运转时间内发生故障的平均次数S:总损失费用E:机器故障时间累计F:工人空闲时间累计K:最短时间机器的当前时刻H:最短时间机器对应的机器号L:工人修完机器的释放时刻x:最短释放时间工人的当前时刻

13、y:最短释放时间对应的工人序号,事件步长法-例子,下表列出了当 m=86，1/V=500小时，1/U=34小时，c1=3.46元/小时，c23.2元/小时时的仿真结果，其中每次预定仿真时间为一万个小时，连续仿真五次取其平均值作为仿真结果。由仿真结果可知当工人数c10时，最优工人数为7,应用举例-可靠性问题,例：一设备上三个相同的轴承，每个轴承正常工作寿命为随机变量，其概率分布如表所示,任何一个轴承损坏都可以使设备停止工作，从有轴承损坏，设备停止工作，到检修工到达开始更换部件为止，称为一个延迟时间延迟时间也是随机变量，其概率分布如下表所示,应用举例-可靠性问题,设备停工时每分钟损失5元,检修工每

14、小时工时费12元，轴承每个成本 16元更换一个轴承需要 20 min，同时更换两个轴承需要30min，同时更换三个轴承需要40min,现在有两种方案：方案一是损坏一个更换一个；方案二是一旦有轴承损坏就全部更换试通过计算机仿真对这两种方案做出评价,在这一问题中，轴承寿命在1000到1900h之间，而延迟在5到 15min之间，故若用时间步长法时，步长选取有些困难步长小浪费很大，步长大又不精确，所以采用事件步长法在事件发生时再考虑系统状态的变化情况，这就比较合理,应用举例-可靠性问题,为了进行仿真，首先对轴承寿命和延迟时间与随机数对应，对应规则分别如下两个表,轴承寿命随机数,延迟时间随机数,应用举

15、例-可靠性问题,由于在这一问题中各个轴承的寿命完全决定了系统的运行状态，也即决定了两个方案的费用大小，故我们选择轴承发生故障作为事件，这三个轴承发生故障的事件分别记为A、B、C,（1）方案一的仿真,产生初始事件表,应用举例-可靠性问题,仿真时钟步进，计算费用，产生下一个事件由表上看出，最早发生的事件是A，所以,t 1400hcost(5+20)54116145元，,下一个A事件发生的时刻为第2 400小时25分钟（随机产生轴承寿命为1000h），刷新事件表，即删去老的A事件，产生新的A事件刷新后的事件表如下表,应用举例-可靠性问题,寻找事件表中的最早事件进行处理,由上表看出，B、C事件同时发生

16、在第1500 h，故同时处理时钟步进为 t 1500 h，再根据费用的计算方法得：,cost=145(15+30)5 6216408元 最后利用随机数产生新事件B和C，刷新事件表，得到的新事件表:,应用举例-可靠性问题,重复，t=2400 h 25 min，cost=408(5+20)5 4 116=553元,重复这一过程，一直需要的时间结束即可得到方案的费用,应用举例-可靠性问题,（2）方案二的仿真,方案二与方案一的区别就是一旦故障发生，就更换3个轴承设初始事件表仍为方案一初始表表中最早的事件是A，处理事件A时要考虑延时，更换3个轴承的时间和费用,cost=(5+40)58316281元,应

17、用举例-可靠性问题,根据下一次三个轴承发生故障的时刻刷新后的事件表:,再重复上述过程，累加费用，即可得到方案二的总费用,最后比较两种方案的费用大小即可确定选取那一种,应用举例-可靠性问题,程序运行结果：T=100000h方案一：cost=32705元方案二：cost=24429元从而得出方案二较方案一优。,应用举例-可靠性问题,/可靠性问题#include#include#include#define TIMES 100000 struct table float begin;int delay;int make;struct table tab3;float cost=0;int sort3

18、;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+=1;if(tab2.begin=tab

19、1.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;,应用举例-可靠性问题,int make_eve

20、nt2(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.begin)sorti+=2;else if(

21、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(xbegin+=1400;else if(xbegin+=1

22、500;else if(xbegin+=1800;else if(xbegin+=1700;else if(xbegin+=1800;else if(xbegin+=1900;for(t1=,y=random(1000);y/=1000;if(ydelay=5;else if(ydelay=10;else if(ydelay=15;return;,应用举例-可靠性问题,void create_event2(struct table*t,struct table*t0)float x,y;struct table*t1;x=random(1000);x/=1000;if(t0!=0)t-begi

23、n=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(xbegin+=1800;else if(xbegin+=1900;for(t1=,y=random(1000);y/=1000;if(ydelay=5;el

24、se if(ydelay=10;else if(ydelay=15;return;,应用举例-可靠性问题,void printf_result(void)for(int i=0;iTIMES);printf(1:tCost=%f T=%fn,cost,T);,应用举例-可靠性问题,void main2(void)int i;float T=0;memset(tab,0,sizeof(tab);cost=0;randomize();create_event2(,固定增量法,选择适当的时间单位T作为仿真钟推进进的增量，每推进一步进行如下处理：1）该步内若无事件发生，则仿真钟再推进一个时间单位；2）

25、若在该步内有若干个事件发生，则认为这些事件均发生在该步的结束时刻。缺点是：仿真钟每推进一步都要检查事件表以确定是否有事件发生，增加了执行时间；将发生的事件均视为发生在这一步的结束时刻，如果T选得过大，会引入较大的误差；且要求确定各类事件处理的顺序，增加了建模的复杂性。主要用于系统事件发生时间具有较强周期性的模型，如定期订货的库存系统，以年、月为单位的经济计划系统等。,应用举例-报童的策略,例 报童每天清晨从报社购进报纸零售,晚上将没有卖掉的报纸退回.每份报纸的购进价为1.3元,零售价为2元,退回价为0.2元.报童售出一份报纸赚0.7元,退回一份报纸赔1.1元.报童每天如果购进的报纸太少,不够卖

26、时会少赚钱,如果购得太多卖不完时要赔钱.试为报童筹划每天应如何确定购进的报纸数使得收益最大.报纸每捆10张，只能整捆购买，报纸可以分为3种类型的新闻日：好、一般、差，它们的概率分别为0.35,0.45和0.2,在这些新闻日中每天对报纸的需求分布的统计结果下图：,应用举例-报童的策略,试确定每天报童应该订购的报纸数量,应用举例-报童的策略,解：我们通过计算机仿真来解决此问题。最优策略应该是每天的利润最大。利润=销售收入-报纸成本-损失+残值这是一个随机现象的计算机仿真问题，故先确定各种情况的随机数的对应关系。新闻日和需求量对应的随机数分别如下面两个表格所示,应用举例-报童的策略,计算机仿真的流程

27、：1）令每天的报纸订购数变化，40-100；2）让时间从1开始变化（循环）到365；3）产生新闻种类的随机数，确定当天的新闻类型；4）产生需求量随机数，确定当天的报纸需求量；5）计算当天的收入，计算累积利润，8）比较得出最优定货量。,应用举例-报童的策略,具体的计算由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)

28、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;else 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)=news sale=news

29、;实际销售量 remand=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份，平均每

30、天利润34.4元。,8.3 随机数的产生,现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响要使我们的数学模型能够较真实地刻画实际对象，必须面对这个现实概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具但是它还难以用来处理复杂系统中的随机性这里介绍使用计算机来模拟随机现象的方法，它经常应用于复杂系统的动态仿真的研究当中是处理复杂系统中随机性的计算机模型也是使用计算机研究和解决实际问题的一条重要途径,随机现象的模拟,1均匀分布的随机数及其产生,对随机现象进行模拟，实质上是要给出随机变量的模拟也就是说利用计算机随机地产生一系列数值，它们的出现服从一定的概率分布,则称这些

31、数值为随机数。最常用的是在(0,1)区间内均匀分布的随机数，也就是我们得到的这组数值可以看作是(0,l)区间内均匀分布的随机变量的一组独立的样本值以后我们将指出其它分布的随机数可利用均匀分布的随机数产生,随机现象的模拟,记 XU(0,1),0,1)区间上均匀分布的随机变量X的概率密度f(x)和概率分布函数F(x)分别为:,数学期望:E(X)=1/2;方 差:=1/12,随机现象的模拟,下图为均匀分布的密度函数和分布函数示意图：,随机现象的模拟,均匀随机数是产生其它随机数的基础。例如，抛硬币、抽签、统计经验分布都可以由它产生。,产生随机数的方法：,（1）硬件设备：从真实物理现象的随机因素中产生随

32、机数，放射性粒子的放射源，电子晶体管的固有噪音等，单位时间内放射出的粒子数是随机的。优点：真正的随机数；缺点：外部设备，无法重复,随机现象的模拟,（2）数学公式：产生伪随机数 用数学公式或位移寄存器的移位操作来产生的随机数,称为伪随机数.因为真实的随机数,只能从客观真实的随机现象本身产生出来,所以产生理想的伪随机数列不是一件容易的事.,一般对于产生伪随机数的方法,有如下几点要求:1)要求伪随机数列有较理想的随机性和均匀性,就是对其随机性和均匀性进行统计检验时,有合乎要求的精度;2)产生伪随机数的程序应当简短、运算速度快,占用计算机的内存单元少;3)伪随机数列的循环周期应当尽可能地大,以满足模拟

33、的需要4)伪随机数列中,前后之间和各子列之间,要求相互是独立无关的。,随机现象的模拟,当我们需要一系列均匀分布的随机数时,可以按照一定的算法由计算机计算产生的伪随机数产生随机数的方法很多，其中以乘同余法使用较广用以产生均匀分布随机数的乘同余法的递推公式为：,xn+1=xn(mod M),rn=xn/M,其中，是乘因子，M是模数前面式子的右端称为以M为模数（modulus）的同余式。给定了一个初值x0(称为种子)后，计算出的rn即为(0,l)上均匀分布的随机数,随机现象的模拟,无论用哪一种方法产生的随机数都存在这样的问题，必须对它进行统计检验看看它们是否具有较好的独立性和均匀性一般在计算机(或计

34、算器)及其使用的算法语言中都有随机数生成的命令，它们所生成的随机数都是经过检验并且可用的.这里就不再详细介绍检验的方法了,随机现象的模拟,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,1)标准正态分布随机数randn(m,n)m*n个N(0,1)标准正态分布随机数,随机现象的模拟,2随机变量的模拟,利用均匀

35、分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟,（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,n,i=1,随机现象的模拟,由此可知，事件(p(n-1)R p(n)和事件(X=xn)有相同的发生的概率因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是

36、：每产生一个(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有相同的分布由此，若已知Y的概率密度为f(y)，由Y=F-1(X)可得,XF(Y)=f(y)dy,Y,-,随机现象的模拟,是区间（0，l）上均匀分布的随机变量如果给定

37、区间（0，1）上均匀分布的随机数ri，则具有给定分布的随机数yi可由方程中解出,ri F(Y)=f(y)dy,yi,0,例当需要模拟服从参数为的指数分布时,由,ri e-xdy=1-e-yi,yi,0,可得yi=-(1/)ln(1-ri).因为(1-ri)和ri 同为(0，1)区间上的均匀分布的随机数，故上式可以简化为yi=-(lnri)/.,随机现象的模拟,反函数法是一种普通的方法，但当反函数不存在或难以求出时，就不宜于使用了舍选法是另一种方法，其实质是从许多均匀随机数中选出一部分，使之成为具有给定分布的随机数步骤如下：,设随机变量X有密度f(x)，又存在实数ab，使 Pr(aXb)=1.,

38、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(a11-其中 是充分小的正数重复上面的步骤，所产生的随机数仅会出现较小的误差,随机现象的模拟,（3）正态随机数的模拟,除了可用反函数和舍选法模拟正态随机变量外,还有两种常用的方法,坐标变换法:设r1,r2是相互独立的均匀随机数,令 x1(-2lnr

39、1)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/12,l,i=1,n,n,2,n/12,l,i=1,n,n,2,随机现象的模拟,将可看成是标准正态分布N(0,1)的随机数这里的n应取得足够大，一般取n=10即可若取n=12,则上式简化为 z=ri-8再由y=x+

40、得到正态分布N(,2)的随机数y.,i=1,n,计算机仿真举例,例 赶火车过程仿真一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站运行时间为均值30分钟、标准差为2分钟的正态随机变量.火车大约在下午1点离开A站。离开时刻的频率分布为,计算机仿真举例,用计算机仿真火车开出、火车到达B站、这个人到达B站情况，并给出能赶上火车的仿真结果。,引入以下变量：T1 火车从A站开出的时刻；T2 火车从A站运行到B站所需要的时间；T3 此人到达B站的时刻；,这个人到达B站时的频率分布为：,计算机仿真举例,T1,T2,T3 是随机变量，其概率分布为 x1=0.7,x2=0.9,y

41、1=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=s3+1;end ends1/10000,1-s1/10000-s3/10000,s3/10000,计算机仿真举例,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 enden

42、ds1/10000,s2/10000,s3/10000,s4/10000,人到达时刻仿真测试,计算机仿真举例,火车运行时间的仿真测试x=randn(10000,1);for i=1:10000 y(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)0.9 T1=5;else T1=10;endT2=30+2*x3(i);,if x2(i)0.3 T3=28;elseif x2(i)0.7 T3=30;elseif x2(i)0.9 T3=32;else T3=34;end if T3T1+T2 s=s+1;endends/10000,