计算机仿真原理及应用第三讲 离散事件系统仿真方法ppt课件

单位 物理电子学院主讲人 王刚Email buncan wang 计算机仿真原理与应用 3离散事件系统仿真方法 离散事件系统示例例 单人理发店系统 设上午9 00开门 下午5 00关门 顾客到达时间一般是随机的 为每个顾客服务的时间长度也是随机的 系统的状态 服务台的状态 忙或闲 顾客排队等待的队列长度 状态量的变化也只能在离散的随机时间点上发生 系统 一些具有特定功能 相互之间以一定的规律联系着的物体所组成的总体 系统边界 为了限制所研究问题涉及的范围 用系统边界把所研究的系统与影响系统的环境区分开来 实体 系统的对象 系统的组成元素都可以称为实体 是仿真系统中可单独识别和刻划的构成要素 属性 属性反映实体的某些性质 是实体特征的描述 用特征参数或变量表示 它可以是文字型 数字型或逻辑型 活动 实体在一段时间内持续进行的操作或过程 通常用于表示两个可以区分的事件之间的过程 标志着系统状态的转移 状态 系统的状态是指在某一时刻实体及其属性值的集合 事件 事件是引起离散事件系统状态发生变化的行为 进程 由若干个有序事件及若干有序活动组成 描述了它所包括的事件及活动间的相互逻辑关系及时序关系 专用术语 仿真策略 要将系统模型转换为计算机模型 首先要从总体上确定仿真模型的控制逻辑和仿真时钟推进机制 即确定仿真策略 仿真策略是仿真模型的核心 反映了仿真模型的本质 从根本上决定了仿真模型的结构 迄今为止 离散事件系统形成了三种基本仿真策略 事件调度法 eventschedule ES 活动扫描法 activityscanning AS 进程交互法 ProcessInteractive PI 事件调度法由兰德公司在1963年提出 在美国广泛采用 欧洲不很流行 基本思想 将事件例程作为仿真模型的基本模型单元 按照事件发生的先后顺序不断执行相应的事件例程 每一个有确定发生时间的事件 都有一个事件例程 用事件例程来处理事件发生后对实体状态所产生的影响 并安排后续事件 事件调度法用事件的观点分析真实系统 通过定义事件及每个事件的发生 引起系统状态的变化 按时间顺序 在每个事件发生时 确定并执行有关的逻辑关系 按这种策略建立模型时 所有的事件均放在事件表中 模型中设有一个时间控制器 从事件表中选择最早发生的事件并将仿真时钟改到该事件发生的时间 然后调用与该事件相应的事件处理模块 这样 事件的选择与处理不断地进行 直到仿真终止或程序事件产生为止 事件调度法的仿真策略 1 初始化置仿真的开始时间t0和结束时间tf 置实体的初始状态 置初始事件及其发生时间ts 2 仿真时钟TIME ts 3 确定在当前时钟TIME下发生的事件类型Ei i 1 2 n 并按规则排序 4 如果TIME tf 执行 caseEiofE1 执行E1的事件例程 产生后续事件类型及发生时间 En 执行En的事件例程 产生后续事件类型及发生时间 endcase 否则 转 6 事件调度法的仿真策略 5 将仿真时钟TIME推进到下一个最早事件发生时间 这一步体现了仿真时钟的推进机制 是将仿真时钟推进到下一个最早事件的发生时刻 它与连续系统仿真中的时间推进方法 固定时间增量法不同 反映了离散事件系统状态仅在离散时刻点上发生变化的特点 这种时间推进方法在离散事件系统仿真中普遍采用 称为下一事件增量法 简称事件增量法 6 仿真结束 1 系统建模由观测数据确定随机变量的分布和参数 一般可用流程图或网格图的方式描述 反映临时实体在系统内部历经的过程 永久实体对临时实体的作用以及它们相互之间的逻辑关系 2 确定仿真算法两个方面内容 如何产生所需求的随机变量 采用什么方法对离散事件系统仿真 确定仿真策略 3 建立仿真模型仿真时钟在各种算法中的推进方法 根据仿真算法建立仿真系统的计算机模型 变量定义及程序流程 4 设计仿真程序实现仿真模型 5 仿真结果分析每次仿真结果只是随机变量的一次取样 仿真结果的可信度 离散事件系统仿真步骤 已知的基本信息 等待区足够大 排队规则先进先出FIFO 到达间隔服从负指数分布 1 1 10 台 天 修理时间服从负指数分布 2 1 15 台 天 仿真时间长度为365天 编程序求解 机器的平均等待时间 机器的平均逗留时间 修理台利用率 实例 机器修理车间的仿真 这是一个典型的单服务员单队列的排队系统仿真模型 这类排队系统主要包括两个要素 顾客 即服务对象 和服务员 即服务设备 该系统由到达模式 服务模式 并行服务员数目 系统容量 排队规则来表示 由命题可知 被修理的机器为 顾客 而修理台为 服务员 该排队系统的到达模式用机器到达间隔时间的负指数分布表示 服务模式由修理时间的负指数分布表示 系统中并行服务员数目为1 系统容量足够大 排队规则采用先进先出FIFO方式 系统描述 系统建模 采用事件调度法 具体的仿真步骤如下 1 初始化 给出当前仿真时钟 系统状态量及统计量的初始值 2 扫描事件表 将当前仿真时钟增加到下一个最早发生事件的时间上 3 处理该事件 相应地改变系统状态 4 收集统计数据 5 若仿真时间未结束 则返回2 否则 执行下一步 6 分析收集的统计数据 产生报告 通过分析可知 该仿真模型只存在两类事件 第一类事件为 到达事件 第二类事件为 离开事件 那么下一事件的类型由变量EVTFLAG给出 仿真模型的总体结构图如下页所示 其中INIT为系统初始化子程序 TIMEDV为时间推进子程序 ARRIVE为到达事件处理子程序 DEPART为离开事件处理子程序 REPORT为报告生成子程序 仿真模型的总体结构图 建模变量说明 到达事件的处理流程 离开事件的处理流程 计算机仿真结果 由已知条件可知 到达间隔时间服从 1 1 10 台 天 的负指数分布 修理时间服从 2 1 15 台 天 的负指数分布 仿真时间长度为365天 故到达间隔时间均值EATI 1 1 10 天 修理时间均值ERT 1 2 15 天 仿真结束时间TIME 365 天 给定随机数发生器种子SEED 113 则仿真结果机器在系统中平均逗留时间为33天 机器在队列中平均等待时间为40天 修理台利用率为78 9 基本思想 针对待求问题 根据物理现象本身的统计规律 或人为构造一合适的依赖随机变量的概率模型 使某些随机变量的统计量为待求问题的解 进行大统计量 N 的统计实验方法或计算机随机模拟方法 理论依据 大数定理 均匀分布的算术平均收敛于真值中心极限定理 置信水平下的统计误差 两个例子 Buffen投针实验求 射击问题 打靶游戏 蒙特卡罗方法概述 基本思想 Buffon投针实验 1777年 求 N 大数定理 蒙特卡罗方法概述 基本思想 针与平行线垂直方向夹角为 则相交概率为 一些人进行了实验 其结果列于下表 1 设r表示射击运动员弹着点到靶心的距离 r 表示击中r处相应的得分数 环数 f r 为该运动员弹着点的分布密度函数 它们反映运动员的射击水平 该运动员的射击成绩为 用概率语言来说 是随机变量 r 的数学期望 即 假设该运动员进行了 次射击 每次射击的弹着点依次为r1 r2 rN 则 次得分g r1 g r2 g rN 的算术平均值 代表了该运动员的成绩 射击问题 打靶游戏 用 次试验所得成绩的算术平均值作为数学期望的估计值 例如 设射击运动员的弹着点分布为 用计算机作随机试验 射击 的方法为 选取一个随机数 按右边所列方法判断得到成绩 这样 就进行了一次随机试验 射击 得到了一次成绩 r 作 次试验后 得到该运动员射击成绩的近似值 收敛性 大数定理 作为所求解的近似值 由大数定律可知 如果X1 X2 XN独立同分布 且具有有限期望值 E X 则 即随机变量X的简单子样的算术平均值 当子样数 充分大时 以概率1收敛于它的期望值E X 由前面介绍可知 蒙特卡罗方法是由随机变量X的简单子样X1 X2 XN的算术平均值 f x 是X的分布密度函数 则当N充分大时 有如下的近似式 蒙特卡罗方法的近似值与真值的误差问题 概率论的中心极限定理给出了答案 该定理指出 如果随机变量序列X1 X2 XN独立同分布 且具有有限非零的方差 2 即 其中 称为置信度 1 称为置信水平 这表明 不等式 近似地以概率1 成立 且误差收敛速度的阶为 O N 1 2 蒙特卡洛法中的误差 中心极限定理 上式中 与置信度 是一一对应的 根据问题的要求确定出置信水平后 查标准正态分布表 就可以确定出 通常 蒙特卡罗方法的误差 定义为 给出几个常用的 与 的数值 两点说明 1 MC方法的误差为概率误差 这与其他数值计算方法是有区别的 2 误差中的均方差 是未知的 必须使用其估计值来代替 在计算所求量的同时 可计算出估计值 2 减小估计的均方差 比如降低一半 那误差就减小一半 这相当于N增大四倍的效果 减小方差的各种技巧 1 增大试验次数N 在 固定的情况下 要把精度提高一个数量级 试验次数N需增加两个数量级 因此 单纯增大N不是一个有效的办法 显然当给定置信度 后 误差 由 和N决定 要减小 效率 一般来说 降低方差的技巧 往往会使观察一个子样的时间增加 在固定时间内 使观察的样本数减少 所以 一种方法的优劣 需要由方差和观察一个子样的费用 使用计算机的时间 两者来衡量 这就是MC方法中效率的概念 它定义为 2c 其中c是观察一个子样的平均费用 显然 2c越小 方法越有效 1 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程 蒙特卡罗方法的优点 从这个意义上讲 蒙特卡罗方法可以部分代替物理实验 甚至可以得到物理实验难以得到的结果 用蒙特卡罗方法解决实际问题 可以直接从实际问题本身出发 而不从方程或数学表达式出发 它具有直观 形象的特点 2 受几何条件限制小 计算s维空间中的任一区域Ds上的积分 无论区域Ds的形状多么特殊 只要能给出描述Ds的几何特征的条件 就可以从Ds中均匀产生N个点 因此 在具有随机性质的问题中 如考虑的系统形状很复杂 难以用一般数值方法求解 而使用蒙特卡罗方法 不会有原则上的困难 其中Ds为区域Ds的体积 这是数值方法难以作到的 得到积分的近似值 3 收敛速度与问题的维数无关 由误差定义可知 在给定置信水平情况下 MC方法的误差为O N 1 2 与问题本身的维数无关 维数的变化 只引起抽样时间及估计量计算时间的变化 不影响误差 这一特点 决定了蒙特卡罗方法对多维问题的适应性 而一般数值方法 比如计算定积分时 计算时间随维数的幂次方而增加 而且由于分点数与维数的幂次方成正比 需占用相当数量的计算机内存 这些都是一般数值方法计算高维积分时难以克服的问题 4 具有同时计算多个方案与多个未知量的能力 2 使用蒙特卡罗方法还可以同时得到若干个所求量 1 对于那些需要计算多个方案的问题 使用蒙特卡罗方法有时不需要像常规方法那样逐个计算 而可以同时计算所有的方案 其全部计算量几乎与计算一个方案的计算量相当 例如对于屏蔽层为均匀介质的几何平板 要计算若干种厚度的穿透概率时 只需计算最厚的一种情况 其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到 例如在模拟粒子过程中 可以同时得到不同区域的通量 能谱 角分布等 而不像常规方法那样 需要逐一计算所求量 5 误差容易确定 根据蒙特卡罗方法的误差公式 可以在计算所求量的同时计算出误差 6 程序结构简单 易于实现 在计算机上进行蒙特卡罗方法计算时 程序结构简单 分块性强 易于实现 1 收敛速度慢 2 误差具有概率性 蒙特卡罗方法的收敛速度为O N 1 2 一般不容易得到精确度较高的近似结果 对于维数少 三维以下 的问题 不如其他方法好 由于蒙特卡罗方法的误差是在一定置信水平下估计的 所以它的误差具有概率性 而不是一般意义下的误差 1蒙特卡罗方法概述 MC缺点 3 计算结果与系统大小有关 例如在粒子输运问题中 经验表明 只有当系统的大小与粒子的平均自由程可以相比较时 一般在十个平均自由程左右 蒙特卡罗方法计算的结果较为满意 但对于大系统或小概率事件的计算问题 计算结果往往比真值偏低 在使用蒙特卡罗方法时 可以考虑把蒙特卡罗方法与解析 或数值 方法相结合 取长补短 既能解决解析 或数值 方法难以解决的问题 也可以解决单纯使用蒙特卡罗方法难以解决的问题 真随机数 不可重复 物理方法产生 随机数是实现由已知分布抽样的基本量 在由已知分布的抽样过程中 将随机数作为已知量 用适当的数学方法可以由它产生具有任意已知分布的简单子样 由具有已知分布的总体中抽取简单子样 在蒙特卡罗方法中占有非常重要的地位 总体和子样的关系 属于一般和个别的关系 或者说属于共性和个性的关系 随机数的产生和检验 蒙特卡罗模拟就是边产生随机数边进行随机模拟的方法 准随机数 不具随机性质 只要处理问题能得到正确结果 如放射性衰变 电子设备的热噪音 宇宙射线的触发时间等 伪随机数 可重复 数学方法产生 必须通过统计检验 区分 数列的随机性和随机数的分布 一个完美的随机数序列可能具有某种分布 如均匀分布 高斯分布等 但具有某种分布的数列却可能完全不是随机的 分布函数为 最简单 最基本 也是最重要的随机数是在单一区间 0 1 上的均匀分布的随机数 其分布密度函数为 由于随机数在蒙特卡罗方法中占有极其重要的位置 我们用专门的符号 表示 用 1 2 代表相互独立且具有相同单位均匀分布的随机数序列 独立性 均匀性是随机数必备的两个特点 如 掷筛子游戏 投掷硬币 用来作为随机数发生器的物理源主要有两种 一种是根据放射性物质的放射性 另一种是利用计算机的固有噪声 用物理方法产生的随机数序列无法重复实现 不能进行程序复算 给验证结果带来很大困难 而且 需要增加随机数发生器和电路联系等附加设备 费用昂贵 因此 该方法也不适合在计算机上使用 在计算机上用物理方法产生随机数的基本原理是 利用某些物理现象 在计算机上增加某些特殊设备 可以在计算机上直接产生随机数 这些特殊设备称为随机数发生器 1 冯 诺伊曼平方取中法 递推公式 以十进制数为例 平方取中法是把一个2S位的十进制自平方后 去头截尾只保留中间2S个数字 然后用102S来除 这样就可以得到在 0 1 上均匀分布的伪随机数序列 例如 设十进制数的2S 4 并取x1 6406 则有 相应的伪随机数序列是0 6406 0 3680 0 1354 0 8333 0 4388等 具有周期性 有些数甚至会紧接着重复出现 很少使用 由Lehmer在1951年提出来的 它的一般形式是 对于任一初始值x0 伪随机数序列由下面递推公式确定 2 Lehmer线性同余法 例如 乘同余法 x0称为种子 改变它的值就得到基本序列的不同区段随机数 a 乘子 c 增量 m 模 乘同余法具有在计算机上容易实现 快速等特点 所以乘同余法已被广泛采用 伪随机数的均匀性 伪随机数的独立性 均匀性是指在 0 1 区间内等长度子区间中随机数的数量是一样的 按先后顺序出现的随机数中 每个随机数的取值与其相距一定间隔的随机数取值之间无关 判断伪随机数序列是否满足均匀和相互独立的要求 要靠统计检验的方法实现 对于伪随机数的统计检验 一般包括两大类 均匀性检验和独立性检验 伪随机数的均匀性 将区间 0 1 分为K个子区间 统计随机数落在第k个子区间的实际频数nk 它应当趋近于理论频数mk 计算统计量 如果 2值很大 表示远远偏离理想值 因此要求 2值尽可能小 但如果它趋于0则有可能N已进入循环 通常求和中的每一项的大小约为1 因此 2的值约为K K的取值不能太大也不能太小 太大反映不出 小区间 的均匀性 太小反映不出 大区间 的均匀性 限制条件 1 概率论中的Pearson定理说明 上式的极限概率分布是 2分布 给出了 2 x的概率 整数 是系统自由度表示独立测量的次数 由于存在一个限制条件 故 K 1 给出了 2 x的概率 余函数 因此 当给定显著水平 或置信度1 后 由方程Q 2 或P 2 1 解出 值 或从已有的 2表中查得 值 如果由 1 式计算出来的 小于 则认为在此置信度下 伪随机数序列在 0 1 中是均匀分布的 1 顺序相关法 它用相邻两个随机数的自相关函数 或相关系数 来标识伪随机数序列的独立性情况 间距为l的自相关函数是 相关系数越小 独立性越好 2 多维频率检验 1 将伪随机数序列用任意一种办法进行组合 每S个随机数作为S维空间中的一个点的坐标值 于是可以构成一个点序列 2 把S维空间中的单位方体分成为K个子方体 方体边长 3 统计落在第k个子方体中的实际频数nk 它应当趋近于理论频数 随机数的独立性统计检验 随机变量的抽样 离散型 随机变量X x1 x2 x3 xN 例如 x可取3个值x1 x2 x3 它们出现的几率分别为2 8 5 8 1 8 则随机数小于2 8时实现x1 在区间2 8 7 8中时实现x2 大于7 8时实现x3 概率密度f p1 p2 p3 pN 如果从 0 1 区间中均匀抽样得到的随机数 满足下式时 则物理量x取值为xn 实际需要的大多数随机变量并不是 0 1 区间均匀分布的 而是有各种不同形式分布密度函数的随机变量 因此 对不均匀的随机变量抽样的关键问题是如何从均匀分布的伪随机变量样本中 抽取符合所要求的分布密度函数的简单子样 设连续型变量x在区间 a b 中取值 可视为将上述的离散情形取连续极限 要求变量x 可由上式解析反解出x 的函数表达式 即求反函数 这对一些简单的几率密度函数解析表达式是很容易做到的 则求反函数后得 随机变量的抽样 连续型 累积函数 注意 1 和 同样服从 0 1 的均匀分布 一维 变换抽样法的基本思想是将一个比较复杂的分布p x 的抽样 变换为已知的简单分布g y 的抽样 我们希望找到x y之间的对应关系 使得几率密度守恒 变换抽样法 例如 黑体辐射的谱密度按频率 表示时为 当希望将谱密度用波长 2 c 表示时 有 显然 当g y 取 0 1 均匀分布时 问题即化为 寻找y x 使其导数为p x 然后在 0 1 区间中对变量y抽样得到均匀分布的随机数 再由x y 关系得到对应几率密度函数p x 的随机抽样x 二维 有两个变量x和y的联合分布密度函数为p x y 欲变换至变量u和v 它们的联合分布密度函数为g u v