离散事件系统仿真基础-课件

1、第四章 离散事件系统仿真基础主要内容基本概念随机变量模型的确定随机数的产生随机变量的产生第一节 基本概念离散事件系统状态仅在离散时间点上变化，且离散时间点一般不确定面向事件；反映系统各部分相互作用的一些事件，模型为反映事件状态的数集，仿真结果是产生处理这些事件的时间历程连续系统：时间常为均匀间隔计时；系统动力学模型由表征系统变量间关系的方程描写，结果常为变量随时间的变化历程例：单服务台排队系统系统工作时间长度固定顾客到达时间随机服务员服务时间随机要求通过仿真估计系统工作情况，以决定是否增加服务台显然，离散事件系统一般有固有的随机性（注意：连续系统也有随机性，如白噪声，但两者行为不同）研究的理论

2、基础：经典的概率及数理统计理论、随机过程理论简单系统可能有理论解析解，但对实际系统，只有靠计算机仿真计算才有可能提供较完整的结果常用概念实体永久实体：永久驻留在系统中，是系统处于活动的必要条件，如服务员临时实体：仅在系统中存在一段时间，按一定规律到达，如顾客关系：临时实体按一定规律不断产生，在永久实体作用下通过系统，最后离开系统事件引起系统状态发生变化的行为离散事件系统本质是由事件驱动的例：顾客到达事件使服务员状态由闲到忙，或使队列长度加1事件的发生一般与某一类实体相联系，放在事件表中管理，事件表通常记录事件类型、发生条件、时间及相关实体的有关属性活动导致系统状态变化的一个过程为活动活动表示两

3、个可区分事件之间的过程，标志着系统状态的转移如顾客到达事件与顾客开始接受服务事件之间为一活动，使服务员忙及队列长度减1进程相当于系统的子集或子系统，包含若干个事件及活动，并且描述了其所包含事件及活动间的逻辑关系和时序关系如某一顾客在系统中的全部活动为一进程事件、活动、进程的关系图仿真时钟离散事件动态系统的状态本来就只在离散时间点上发生变化，因而不需要进行离散化处理。离散事件系统一般不以时间推动，但事件间有时序关系，仿真中仍必须有控制时间的部件由于引起状态变化的事件发生时间的随机性，仿真钟的推进步长则完全是随机的两个相邻发生的事件之间系统状态不会发生任何变化，因而仿真钟可以跨过这些“不活动”周期

4、,仿真钟的推进呈现跳跃性，推进速度具有随机性。统计计数器因固有的随机性，某一次仿真运行得到的状态变化过程只不过是随机过程的一次取样，离散事件系统的仿真结果只有在统计意义下才有参考价值在仿真模型中, 需要有一个统计计数部件， 以便统计系统中的有关变量，如排队系统中的顾客等待时间、队列长度等离散事件系统仿真的一般步骤系统建模：一般用流程图描述，反映临时实体在系统内部历经的过程、永久实体对临时实体的作用及相互间逻辑关系关键：确定随机变量的模型确定仿真算法产生随机变量确定仿真建模策略事件调度法：面向事件建立仿真模型活动扫描法：面向活动建模进程交互法：面向进程建模三阶段法：结合活动扫描与事件调度图形仿真

5、方法：Petri网建立仿真模型定义状态变量、定义系统事件及有关属性、活动及进程、设计仿真钟的推进方法等仿真程序设计及运行仿真语言或高级语言长期运行或多次运行仿真结果分析统计结果、可信度分析等第二节 随机变量模型的确定无序中蕴含着有序，随机过程也有数学描述形式，可近似归纳总结为几种变量分布模式，使定量研究成为可能没有绝对的无序和有序，如混沌以单服务台排队系统中顾客到达时刻为例，总可以找到一种接近的随机变量分布通常需要从观测数据中寻找规律在寻找分布形式时，根据对随机变量（Random variable, r.v.）的特性了解程度，一般会遇到三种情况r.v.分布类型已知，需要由观测数据确定分布参数需

6、要由观测数据确定概率分布类型及参数难以由观测数据确定理论分布形式，需要定义实验分布一、分布参数的确定分布参数的类型 定义分布所采用的大多数参数，由物理或几何解释，可分为三个基本类型位置参数比例参数形状参数位置参数确定了一个分布函数取值范围的横坐标当改变时，分布函数仅平移而无其它变化，又称位移参数例均匀分布函数U(a,b)，密度函数其中a，b均可定义为位置参数比例参数决定分布函数在其取值范围内取值的比例尺的改变只压缩或扩张分布函数，不改变基本形状例：指数分布函数EXPO()，密度函数形状参数确定分布函数的形状，从而改变分布函数的性质例：韦伯分布Weibull(,)，密度函数：分布参数的估计常用方

7、法最大似然估计(maximum likelihood estimation)在已经得到试验结果的情况下，应该寻找使这个结果出现的可能性最大的那个参数值作为真值的估计最小二乘估计(least-square estimation)无偏估计(unbiased estimation)最大似然估计法估计分布参数讨论一个未知参数的情形，设观测数据为离散分布情形：可令 为该分布的概率质量函数，定义似然函数L()为： 的最大似然估计值 使L()取最大值连续分布情形：令 为概率密度函数，定义似然函数为连续分布例：指数分布密度函数为，被估计的参数由为求使取最大值的，取对数因严格递增，取最大值等价于最大求极值考虑到

8、当时，由于为正，显然为最大值离散分布例：泊松分布被估计参数质量函数由令故为最大似然估计值对多个参数的估计问题，上述方法原则上适用，只是计算更繁琐，有时不能解析计算，必须用数值计算方法常见分布的参数的最大似然估计均匀分布 密度函数 参数正态分布 密度函数 位置参数 比例参数二项分布 质量函数t为正整数， ，t已知时二、分布类型的假设对观测数据进行预处理后，才能对分布类型进行假设连续分布数据的预处理方法点统计法、直方图法、概率图法离散分布数据的预处理方法点统计法、线图法定义实验分布连续分布类型的假设（1）点统计法：基于连续分布的偏差系数特征进行分布类型假设偏差系数 （方差、均值）常见分布的偏差系数

9、取值范围指数分布均匀分布 ， 取值 ，除0以外正态分布 ， ，取值同上预处理 ，则显然不能唯一确定分布类型，且计算不一定无偏，但对确定理论分布具有指导作用连续分布类型的假设（2）直方图法将 的取值范围等分为k个区间 ，令 表示第j个区间上观测数据个数 占总数的比例，即定义函数作出h(x)的直方图，并与理论分布的密度函数图形比较，选择接近的分布困难：区间长度b的取值太大：丢失很多信息太小：不能抑制观测噪声一般应多选几个b，取使直方图包络线较光滑的区间长度连续分布类型的假设（3）概率图法通过比较分布函数确定理论分布设 有m个取值，记为x(1)x(m)定义实验分布 其中 表示小于或等于x(i)的观测

10、数据的个数，显然为避免由有限个观测数据得到的实验分布函数值1，修正为概率图不能直接比较（基本都是S型），而是采用“分位点”比较法定义：设0g0时称混合乘同余法，c=0称乘同余法素数取模乘同余法（PMMLCG）m是小于2b的最大素数，而a的选择满足al-1被m整除的最小整数l=m-1，也就是说能被m整除的al-1的最小整数为am-1-1，那么得到的zj的周期为m-1，且在每个周期内，1,2,m-1这些整数严格地只出现一次。两个经过检验性能较好的PMMLCG (mod 235 31) (mod 231 1) 二、组合发生器为提高性能，可用一个发生器控制另一个发生器产生随机数，两种方式发生器1产生

11、，发生器2产生1,k上均匀分布的随机整数I，取出并重新产生UI发生器1，2产生 和 ， 在二进制下循环移位 次得 ，再次与 “异或”后得思路：减少递推公式的自相关性，提高独立性加长发生器周期，提高随机数密度和均匀性三、随机数发生器的测试1.均匀性检验常用频率检验：将随机数发生器取值范围0,1分为k个等长区间，由该发生器产生N个随机数，则落在每个子区间上随机数个数理论值为n=N/k，称为理论频率，实际第j个区间上的数据个数nj 总会有偏差，由此可得 ，其大小反映均匀性程度2.独立性检验计算该随机数序列相邻一定间隔的随机数之间的相关系数，判断相关程度先由随机数发生器产生N个随机数Ui，则前后相隔为

12、j个数的相关系数的均值为 其中S2为随机数方差的估计值： 再根据统计理论处理C语言中与随机数有关的函数randomize(void); 初始化x=random(int M); 产生0M 之间的随机数x=rand(void); 产生0215-1之间的随机数Matlab中与随机数有关的函数Rand(n)：n个0,1之间均匀分布随机数Rand(m,n)：m*n 个0,1之间均匀分布随机数randn(n):n个N(0,1)标准正态分布随机数randn(m,n)：m*n个N(0,1)标准正态分布随机数例 射击命中率 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打

13、击，敌方对其阵地进行了伪装并经常变换射击地点经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值说明这是一个概率问题，只涉及到简单的随机数的产生和处理，可以通过理论计算得到相应的概率和期望值。但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程为了能显示我方20次射击的过程，考虑采用计算机模拟1.问题分析需要模拟出以下两件事：(1)观察所对目标的指示正确与否模拟试验有

14、两种结果，每一种结果出现的概率都是1/2因此，可以用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确(2)当指示正确时，我方火力的射击结果模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)这时可用投掷骰子的方法来确定：出现、点：则认为没击中敌人；出现、点： 则认为击毁敌人一门火炮；出现点： 则认为击毁敌人两门火炮2. 符号假设i：要模拟的打击次数；k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击比率； E1：20次射击平均每次毁伤敌

15、人的火炮数3. 程序流图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止硬币正面?YNNY1，2，34，564. 模拟结果从以上模拟结果可以计算出：有效射击 E=7/20=0.35平均值5. 理论计算则由全概率公式：6. 结果比较模拟结果与理论计算虽然结果不完全一致，但仍反映了事件发生的随机性。只要多次实验求平均值，模拟值就会很接近理论值例：赶火车的仿真一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站运

16、行时间为均值30分钟、标准差为2分钟的正态随机变量。火车大约在下午1点离开A站，离开时刻的频率分布为出发时刻（T） 1：00 1：051：10 频率 0.7 0.2 0.1此人到达B站时的频率分布为： 模拟火车开出、火车到达B站、此人到达B站的情况，并给出能赶上火车的仿真结果。引入以下变量：T1 火车从A站开出的时刻；T2 火车从A站运行到B站所需要的时间；T3 此人到达B站的时刻； 0.1 0.2 0.4 0.3 频率1：34 1:32 1:30 1:28到达时刻（T）T1, T2, T3 是随机变量，其概率分布为 x=0.7, x2=0.9, y1=0.3, y2=0.7, y3=0.9

17、0.1 0.2 0.7 P（频率） 10 5 0 T1(分) T3（分） 28 30 32 34 P( 频率) 0.3 0.4 0.2 0.1开车时间的仿真测试 s1=0; s2=0; s3=0; 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 i

18、f x(i)0.9 s3=s3+1; else s4=s4+1; end endends1/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:10000if x1(i)0.7 T1=0; elseif x1(i)0.9 T1=5; else T1=10;endT2=30+2*x3(i);if x2(

19、i)0.3 T3=28;elseif x2(i)0.7 T3=30; else if x2(i)0.9 T3=32; else T3=34; endendif T3T1+T2 s=s+1;endends/10000第四节 随机变量的产生原理对产生随机变量的方法的性能要求准确性：即由这种方法产生的随机变量应准确地具有所要求的分布 快速性要求，在离散事件仿真中，一次运行往往需要产生几万甚至几十万个随机变量，产生随机变量的速度将极大地影响着仿真执行的效率。 产生随机变量的常用方法反变换法、组合法、卷积法及取舍法 一、反变换法常用且直观，以概率积分变换定理为基础设随机变量x的分布函数为F(x)，则x的

20、抽样值可按如下方法求： 产生在0,1上的均匀分布的独立随机变量u，由反分布函数 得到的值即为所需要的随机变量x该方法对分布函数进行反变换，故此得名例：连续随机变量对指数分布 产生x抽样值 解：用随机数发生器产生令则因同分布则有由上面的例子可以看到，用反变换法产生随机变量时首先必须用随机数发生器产生在0,1上均匀分布的独立的u，以此为基础得到的随机变量x才能保证分布的正确性，可见，选择一个均匀性及独立性较好的随机数发生器在产生随机变量中是十分重要的当x是离散随机变量时，其反变换法的形式略有不同，原因在于离散随机变量的分布函数也是离散的，因而不能直接利用反函数来获得随机变量的抽样值对有离散随机变量

21、 x，分别以概率 取值 ， 其中 ，由此可划分n个子区间： 用随机数发生器产生 ，若 则取x=xi即可显然，区间搜索策略直接影响抽样值产生速度例：离散随机变量质量函数及累积分布函数如下 用反变换法产生随机变量x 先由随机数发生器产生0,1区间上均匀分布的随机变量u，设u=0.72，按反变换法，逐区间搜索比较，得x=x3=3二、组合法若分布函数可表示成若干其它分布函数之和，且这些分布函数更易取样时，选择组合法设随机变量x的分布函数 或密度函数 其中 ，且 ， 易于取样 则组合法步骤：(1) 产生一个随机整数 J，满足 即决定采用哪一个分布函数进行取样(2) 产生具有分布函数Fi(x)或密度函数f

22、i(x) 的随机变量，则x=xi第一步可用离散反变换法实现 第二步可用反变换法或以下介绍的其它方法显然，两步都需要产生相互独立的例：设密度函数为产生服从该分布的随机变量 x 直接用反变换法产生随机变量 x很困难，分析f(x)的特点，可知它以纵轴为对称轴分为两部分f(x)可写成如下形式 其中用组合法产生随机变量x的方法如下：由U(0,1)产生随机变量u1及u2若u10.5，则由f1(x)的分布函数产生，可由反变换法得到若 ，则由f2(x)的分布函数产生，即令从而三、卷积法若 ， 为独立同分布的随机变量，则 x的分布函数与 的分布函数相同，称 x的分布为Yi分布的m重卷积为产生 x，可先独立地产生

23、 ，并简单相加即得 x，称为卷积法例：产生均值为的m维Erlang分布随机变量xm维Erlang分布En()的概率密度为 难以直接产生xErlang(m,)可表示为m个均值为 的独立的指数随机变量之和，即 以及解: (1)独立产生m个U(0,1)随机数Ui (2)用反变换法分别产生 (3)令 即可改进算法实际仿真中，为提高算法速度，考虑到对数运算速度较慢，可将(2)、(3)步改进为 (2)计算 (3)则四、取舍法反变换法、组合法以及卷积法有一个共同的特点，即直接面向分布函数，因而称为直接法，它们以反变换法为基础当反变换法难于使用时(例如随机变量的分布函数不存在封闭形式)，取舍法是主要方法之一，

24、其基本思想是根据概率密度函数(质量函数)的几何特征，拟合一个简单的分布形式，以后者产生抽样值并按照某种判据进行取舍基本思想设随机变量 x密度函数为f (x)，f (x)最大值为C，x取值范围0,1若独立产生两个0,1区间均匀分布的随机变量u1，u2，则Cu1是在0,C内均匀分布的随机变量，且满足 的概率为若上式成立，则认为 x=u2，否则拒绝u2显然，抽样成功的概率为f (x)下的面积除以总面积C(见图)，直观上看成功抽样的点符合所需要的分布一般情形，是根据f (x)的特征规定一个函数t(x)，其曲线处于f (x)上方且形状接近f (x) ，再按照判据取舍对t(x)的要求是 易于从t(x)进行反变换如果令 ，则 从而可将r(x)看作是一个密度函数，并用r(x)代替f (x)取样，以得到所需要的随机变量由于r(x)不是要求的f (x)，故此产生选取与舍弃问题，取舍算法如下产生由r(x)独立地产生随机变量u2检验如下不等式 若不等式成立，则令x=u2，否则返回第一步r(x)通常可以选择均匀、指数、正态分布等五、典型随机变量的产生正态分布正态分布 密度函数因分布函数在极坐标下有封闭形式，可采用反变换法设x1，x2是两个独立的N(0,1)随机变量，其联合密度函数转换成极坐标形式则其中|J|