模拟技术教学课件PPT.ppt

模拟技术 1模拟的基本知识2模拟确定性行为3随机数和随机变量的生成4随机模拟案例5模拟结果的统计分析 1 1模拟的概念及作用现实系统的数学或逻辑模型可能十分复杂 例如大多数具有随机因素的复杂系统 其中的一些随机性因素很难用准确的数学公式表述 从而也无法对整个系统采用解析法求解 模拟是处理这类实际问题的有力工具 1模拟的基本知识 monte carlo模拟 对于很难用解析方法加以处理的问题 模拟是一种有效的技术 对建模过程中的假设进行鉴定 对理论研究的结论加以检验 对不同的实现方案进行多次模拟 按照既定的目标函数对不同方案进行比较 从中选择最优方案 模拟通常借助于计算机进行 计算机模拟 在已经建立的数学 逻辑模型的基础之上 通过计算机试验 对一个系统按照一定的决策原则或作业规则 由一个状态变换为另一个状态的行为进行描述和分析 蒙特卡罗 monte carlo 方法 也称为随机模拟 randomsimulation 源于二战期间研制原子弹的 曼哈顿计划 s m ulam和j vonneumann用赌城montecarlo来命名这种方法 中子扩散 项目 基本思想 为了解决数学 物理 工程技术等方面的问题 首先建立一个概率模型或随机过程 使它的参数等于问题的解 然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征 最后给出所求解的近似值 运输系统模拟摩天大楼安全疏散系统模拟国民经济发展模拟金融衍生产品的定价及交易风险估算人口增长系统模拟供水系统模拟管理系统模拟雷达系统模拟战争系统模拟 应用领域 模拟的主要优点和缺点 通常 模拟时间是模拟的主要自变量 设计正确的模拟时间推进机理 模拟过程中应根据系统的特性正确推进模拟时间 使系统中各要素与发生的事件保持同步 1 2模拟的分类 下次事件法 将模拟时间由一个事件发生的时间点推进到紧接着的下一次事件发生的时间点 固定时间步长法 模拟时间每次均以相等的固定步长向前推进 每到达一个新的模拟时间点需检查相应时间段内是否发生了事件 需根据实际问题合理设置模拟时间发生改变的步长 推进模拟时间的基本方法 根据模拟过程中因变量的变化情况进行分类 1 离散型模拟 因变量在与事件时间有关的具体模拟时间点呈离散性变化 大多数系统 如排队服务系统 可采用离散型模拟 时间推进方法 一般采用下次事件法应当重点对系统状态可能发生改变的事件进行描述 并确定这些事件之间的逻辑关系 排队系统通常采用离散型模拟模型 其中 发生系统状态变化的事件有两个 一是有顾客到达 二是服务员完成服务 将最近发生上述两种事件之一的时刻设置为下次事件发生点 就可将服务过程描述为图2所示的模拟模型 初始化数据 确定顾客到达时间 有顾客到达 服务员空闲 确定服务时间 是 是 否 否 服务员空闲 服务员开始服务 顾客进入队列等候 推移时间至下次事件发生点 队列中有顾客 否 是 选择一位顾客 该顾客离开队列 图2排队系统中服务过程的模拟模型 2 连续型模拟 因变量随时间的改变呈连续性变化 在大多数计算机模拟过程中 按固定的步长推进模拟时间 通常需建立一系列的由系统状态变量组成的状态方程组 以描述状态变量与模拟时间的关系 3 混合型模拟 因变量随时间的推移而作连续性的变化并具有离散性的突变 如库存控制系统 终态模拟 在规定的时间t内进行模拟运行 时间达到t时 模拟终止 其性能指标明显取决于系统的初始状态 稳态模拟 随着模拟时间的推移 系统的性能逐渐趋于平稳 其目的是研究非终态系统长期运行条件下的稳态性能 模拟时间的长短取决于能否获得系统性能的优良估计 可由模拟输出的精度确定 1 3模拟的方式 明确问题 建立模型 正确描述待研究问题 明确规定模拟的目的和任务 确定衡量系统性能或模拟输出结果的目标函数 然后根据系统的结构及作业规则 分析系统各状态变量之间的关系 以此为基础建立所研究的系统模型 1 4模拟的一般步骤 收集和整理数据资料 模拟技术的正确运用 往往由大量的输入数据作依靠 在随机模拟中 应认真分析具体收集到的随机性数据资料 确定系统中随机性因素的概率分布特性 以此为依据产生模拟过程所必需的抽样数据 编制程序 模拟运行 分析模拟输出结果 模拟结果的统计特性 样本均值 方差 置信区间等 灵敏性分析 选择最优方案 1 频率法 例 估算定积分 构造一个矩形包含了曲边梯形 矩形区域高h满足 考虑的情形 由大数定律 事件发生频率依概率收敛于事件发生的概率 在矩形区域内产生n 足够大 个点 具体思路 落在曲边梯形内的点为m个 则 则所求定积分 其中矩形区域面积 求解程序 求定积分的monte carlo模拟结果 模型中的参数 如矩形区域的高h 对结果有影响吗 思考 这种方法对函数值可正可负的情形该如何应用 2 平均值法 由强大数定律知 以概率1成立 因此当n足够大时 可得近似公式 由此得到计算i的平均值法 3 计算作为i的估计值 1 产生 0 1 上均匀分布随机数 2 令 用样本平均值做总体数学期望的估计 求解程序 例 估算圆周率 模拟产生在正方形abcd中均匀分布的n个点 如果有m个在圆内 则圆的面积与正方形面积之比可近似为m n 即 模拟程序 估算圆周率 n 10000 xy 1 2 rand 2 n x xy 1 y xy 2 mypi 4 sum x 2 y 2 1 n 较简洁实现 模拟运行结果 思考 建立更精确的估算圆周率的数学模型 例 求解非线性规划 在估计的区域 内随机取若干试验点 从试验点中找出可行点 再从可行点中选择最小点 基本假设 试验点的第j个分量服从 内的均匀分布 例 matlab中编程需设计三个m文件 主程序 目标函数程序 约束条件测试程序 例 某湖泊中有机物新陈代谢系统模型可由下面的状态方程组表示 其中 时间t以年为单位 其它符号意义为 t时刻太阳已供能量 t时刻植物生成数量 t时刻吞食植物的虫类生成数量 t时刻食虫植物的生成数量 t时刻湖底的有机物沉淀量 t时刻已扩散于周围环境的总能量 以上各变量的初值为 3 1均匀分布随机数的生成 3 1 1平方取中法 middle squaremethod 1 从一个四位数开始 称为种子 2 将它平方得到一个8位数 必要时前面加0 3 取中间的4位数作为下一个随机数 按照以上方式 就能得到0 9999中的随机整数序列 因而可用来模拟 0 1 上均匀分布的随机数 3随机数和随机变量的生成 缺点 会退化为0或几个数的循环 所以现在一般不采用这种方法 思考 如何判断随机数生成算法的好坏 3 1 2线性同余法 linearcongruence 给定3个正整数a b c和种子 产生规则 特点 问题 序列最终陷入循环 循环周期至多c个数 序列各数之间缺乏统计独立性 c一般取得较大 如 思考 利用线性同余法实现以下功能 产生 0 1 上的随机数 设randint为线性同余法随机数生成函数 产生u a b 上均匀分布的随机数 产生0 1 n上的随机整数 产生从j到k j k 上的随机整数 randint c a b a randint c int n 1 randint c j int k j 1 randint c 3 2matlab随机模拟函数 常见分布随机变量的模拟 1 产生mxn阶 0 1 均匀分布的随机数矩阵 rand m n 2 产生mxn阶 a b 均匀分布u a b 的随机数矩阵 unifrnd a b m n 思考 在matlab中用rand和unifrnd分别产生10个均匀分布u 0 5 的随机数 3 产生mxn阶均值为 标准差为的正态分布的随机数矩阵 randn m n 产生标准正态分布的随机数 当研究现象视为大量相互独立的随机变量之和 且其中每一种变量对总和的影响都很小时 可以认为该现象服从正态分布 4 产生mxn阶期望为的指数分布的随机数矩阵 指数分布在排队论 可靠性分析中有广泛应用 排队服务系统中顾客到达率为常数时的到达间隔 故障率为常数时零件的寿命都服从指数分布 例 顾客到达某商店的间隔时间服从参数为0 1的指数分布 指数分布的均值为1 0 1 10 指两个顾客到达商店的平均间隔时间是10个单位时间 即平均10个单位时间到达1个顾客 顾客到达的间隔时间可用exprnd 10 模拟 5 产生mxn阶参数为的泊松分布的随机数矩阵 poisson分布在排队系统 产品检验 天文 物理等领域有广泛应用 指数分布与泊松分布的关系 如相继两个事件出现的间隔时间服从参数为的指数分布 则在单位时间间隔内事件出现的次数服从参数为的泊松分布 即单位时间内该事件出现k次的概率为 解 先明确随机变量及其分布 单位时间到达数量 poissrnd 4 到达时间间隔 exprnd 1 4 3 3其它随机变量的生成 1 逆变换法 令为 0 1 上均匀分布的随机变量 假设f x 为任一概率分布函数 且f x 存在反函数 令 由的性质知 即为具有分布函数f x 的随机变量 逆变换法的具体步骤 确定随机变量x的概率分布函数f x 产生 0 1 区间上服从均匀分布的随机数 令 求出的x即可作为随机变量x的一个值 例 a b 上的均匀分布 例 指数分布的分布函数为 先产生 0 1 上服从均匀分布的随机数 令 解得 可取 2 一般离散型分布问题 实际应用中 许多随机变量无法用某种规则的分布来近似 此时我们可以通过累积分布进行处理 例 设随机变量的概率分布为 则 的即可作为随机变量的抽样值 对 0 1 区间上均匀分布的随机变量 取满足 模拟方法 例 随机变量函数or累积分布 模拟方法 例 0 3 0 5 0 9 1 0 x f x 0 1 2 3 4 图 累积概率分布图 例 0 3 0 5 0 9 1 0 x f x 0 1 2 3 4 图 累积概率分布图 随机变量的建模 利用理论分布 基于对问题的实际的 合理的假设 选择适当的理论分布模拟随机变量优点是给出了各种理论结果出现的概率 便于进行数学分析和处理 但当问题越复杂 数学处理变得越困难 并且丢失了试验数据的信息 基于实际数据的频率做近似模拟优点在于完全与观察数据相符 并且随实际问题的复杂程度增大不会产生更大的困难 仅增大工作量而已 缺点是不便于进行数学分析 不得不依赖于模拟得到的统计结果 应用常将两种模拟方法结合起来使用 问题 如图 一列火车从a站开往b站 某人每天赶往b站上这趟火车 a b 火车运行方向 某人 他已了解到 1 火车从a站到b站的运行时间是均值为30分钟 标准差为2分钟的随机变量 思考 请研究他能否赶上这趟火车 4模拟模型案例 他到达b站的时刻的频率分布为 他能否及时赶上火车 2 火车在下午大约1点离开a站 离开时刻的频率分布如下 明确问题 他能及时赶上火车的概率是多少 i 分析法 用概率统计知识建立分析模型 求解析解 思考 ii 模拟法 用概率统计知识建立模型 通过模拟求近似解 即先建立模拟模型 然后通过计算机模拟得到问题的近似解 在同样条件下多次试验 计算他能及时赶上火车的频率 建模方向 思路 问题分析 能及时赶上火车的充要条件是 其中t1 火车从a站出发的时刻 t2 火车的运行时间 t3 他到达b站的时刻 是什么变量 如何模拟 ii 将午后1时记为t 0 设火车运行时间t2服从正态分布 t2 n 30 22 基本假设 i 假设t1 t2 t3都是相互独立的随机变量 火车出发时刻t1和人到达b站时刻t3的分布律分别为 建立模型 为了简化计算 将下午1点记为初始时刻 得到随机变量t1和t3的分布律如下 能及时赶上火车的概率p p t3 t1 t2 如果r为在 0 1 均匀分布的随机数 为了模拟随机变量t1和t3 可以通过如下方法 则t1和t3可分别用来模拟随机变量t1和t3 主要变量说明 n模拟次数k临时变量 存储当前累计模拟次数count存储赶上火车的次数 两种不同风格的算法描述 模拟算法设计 输入 赶火车次数 天数 输出 赶上火车的频率 第1步输入模拟次数n第2步k 1 count 0第3步当kt3 count count 1 end第10步执行第3步第11步输出赶上火车频率p count n 算法1 分步骤描述 i 初始化 输入模拟次数n count 0 ii 模拟n次fori 1ton 模拟随机变量t1 t2 t3 分别赋给t1 t2 t3 ift1 t2 t3 count count 1endifendforapp prb count n 算法2 伪代码描述 n input 输入模拟次数 count 0 fori 1 n rt1 rand 模拟随机变量t1 火车从a站出发的时刻 ifrt1 0 7endend forprob count n 模拟程序 模拟结果 每次模拟1000次或5000次 系统模拟注意事项 一次模拟结果毫无意义 模拟是试验性的 是思维结果的验证 必须进行足够多次的模拟 并对结果进行统计分析 系统模拟特点 系统模拟是研究系统 特别是动态系统的重要方法 对于 结构复杂的系统 很难用解析方法求出变量关系的系统 内部机理不明的 黑箱 系统 为验证用其他方法建立的模型及结果 应是较好的选择 练习 有30个电子部件它们的使用情况如下 一开始工作 若损坏则立即开始工作 若损坏则立即开始工作 依次类推 设的寿命服从参数为 小时 的指数分布且相互独立 令t为30个器件的总使用时间 使用m c 蒙特卡罗 模拟法求解和 例 某港口提供有足够的泊位供船舶停靠 但现在仅有一个供装卸的泊位 船舶先到则先进行装卸 假设船舶的到达时间间隔为15 145分钟 而船舶所需的装卸时间为45 90分钟 现要弄清该系统的性能 考察以下问题 1 船舶进入该港后等待装卸的滞留时间 2 等待队列的最大长度 3 卸位 即装卸用的泊位 的利用率 假设 1 假设船舶到达时间间隔和装卸时间服从各自区间上的均匀分布 2 时间取值为整数 在给出一般算法前 考虑五艘船的情形 图 船的等待时间与卸位空闲时间 表 港口系统模拟结果 5艘船情形 tb between 船舶到达港口的时间间隔 ta arrive 一艘船到达港口的时间 tf finish 前一艘船驶离港口的时间 tu unload 当前船舶装卸所需时间 tw waittime 已到港的船舶等待装卸时间总和 ti idletime 卸位已空闲时间总和 maxwait 船舶最大等待时间 maxwaitqueue 船舶等待队列最大长度 所用变量说明 以下考虑一般算法 初始化数据 ta tf tw ti 0 k 1 产生船舶到达间隔时间tb船舶到达时间ta ta tb k n 打印输出 产生装卸服务时间tu ta tf 统计卸位空闲时间ti ti ta tf 统计船舶等候时间tw tw tf ta 设置船舶离港时间tf ta tu 设置船舶离港时间tf tf tu 是 是 否 否 k k 1 表 港口对船舶服务的模拟结果 到达间隔时间 15 145 装卸时间 45 90 表 港口对船舶服务的模拟结果 到达间隔时间 15 145 装卸时间 35 75 到达间隔时间 10 120 装卸时间 35 75 表 港口对船舶服务的模拟结果 对进入该港口的100艘船进行实际调查 记录其活动情况 得到这100艘船到达港口的时间间隔和装卸时间的分布情况 表 100艘船到达港口的时间间隔频数表 表 100艘船装卸时间频数表 表 船舶到达港口的时间间隔累积分布表 表 船舶装卸时间累积分布表 为比较准确地反映系统的性能 我们应采用稳态模拟方式 可以考察系统运行360天的情况 假定港口每天24小时连续工作 模拟时间t设置为8640小时 考虑到船舶到达港口的事件是影响整个系统状态的主要因素 可以将模拟时间设置为该事件的发生时刻 初始化数据t tf tw ti 0 产生船舶到达间隔时间tb设置模拟时间t t tb t 8640 打印输出 产生装卸服务时间tu t tf 统计装卸位空闲时间ti ti t tf 统计船舶等候时间tw tw tf t 设置船舶离港时间tf t tu 设置船舶离港时间tf tf tu 是 是 否 否 到达间隔时间 10 120 装卸时间 35 75 表 港口对船舶服务的模拟结果 经过5次模拟计算 装卸泊位的平均利用率为95 4 但到达的船舶却因得不到及时的服务而造成平均每艘滞留9 30小时 每年到港船舶总滞留时间约6826小时 为此 可修改模拟模型 考虑增设一个装卸泊位 并重新模拟运行 将这些数据提供给决策者 以确定是否需要再新建一个装卸泊位 练习 某公共汽车站每隔30分钟到达一辆汽车 但可能有 0 3 分钟的误差 此误差大小与前一辆汽车的运行无关 汽车最多容纳50名旅客 到达该汽车站时车内旅客人数服从 20 50 的均匀分布 到站下车的旅客人数服从 3 7 的均匀分布 每名旅客下车的时间服从 1 7 秒的均匀分布 旅客按每30分钟到达12个人的泊松分布到达汽车站 单列排队等车 先到先上 如果某位旅客未能上车 他将不再等候 旅客上车时间服从 4 12 秒的均匀分布 上下车的规则是 先下后上 逐个上车 逐个下车 假设每天共发25辆汽车 现在要求模拟一天汽车的运行情况 了解一天中在站内等候汽车的总人数 能上车及不能上车的人数 旅客排队时间分布情况 不能上车人数的分布情况等 例 露天矿用电铲采掘矿石 然后用卡车运送到卸场 假设共有m辆卡车装运 有n台 n m 电铲同时采掘供n辆卡车同时装卸 卸场有s个卸位 可供s辆卡车同时卸料 装运过程以班为单位 每个班开始工作时 m辆中的n辆卡车由n台电铲装车 其余m n辆卡车排队等候 每辆卡车装满矿石后驶离装料处 排于队列首位的待装卡车驶至空闲电铲前并掉转车头接收装料 已装满的卡车驶至卸场时 若卸位空闲则掉头卸完矿石并重新驶回采掘场 排在待装卡车的队尾 1 2 s 卸位 1 2 n 电铲 空载运行 装车队 采场 卸场 重载运行 卸车队 图 露天矿矿石装运系统流程图 假定我们总共模拟20个班运行的情况 已知每个班连续工作6小时 并且每个班开始工作时的初始状态均一致 即全部卡车都在采掘场等待装料 为充分发挥每台电铲 每辆卡车的效率 提高班产量 就需要确定电铲 卡车及卸位之间的一个适当匹配数量 经过对该矿开采情况的收集与处理 得到以下有关数据 1 装车时间 n 1 32 0 27 0 27 2 重载运行时间近似一个常数4分钟 3 卡车掉头时间近似一个常数0 67分钟 4 卸车时间 e 1 0 74 5 卡车空载运行时间近似一个常数3分钟 6 每辆卡车装载量 n 22 5 0 83 0 83 记c i i 1 m 为第i辆卡车 每辆卡车在运行过程中经过装车 重载运行 卸车和空载运行四个阶段 因此 可为每辆卡车c i 定义一个状态c i s 代表卡车所处的各个不同阶段 假定每辆卡车开始时均在采掘场 第一辆卡车开始装料后就开始影响整个系统中各辆卡车的状态 随着时间的推进 各卡车进入各自的阶段 所以为每辆卡车c i 设置一个进入某一阶段的时间c i t 这样 就可将模拟时间推进到所有卡车进入各自阶段时刻的最小值 相应地安排该阶段的服务 服务结束后进入下一阶段并重新设置状态值 为简单计 仅考虑在只有一个电铲和一个卸位 即n s 1 的情况下 应该配备多少辆卡车m才能尽量提高电铲 卡车的效率以及班产量 所用变量说明 t 模拟时间 w 当前工作的班号 c i s 第i辆卡车所处的状态 1 等待装车 2 重载运行 3 等待卸车 4 空载返回 c i t 第i辆卡车切换至c i s状态时的时间 cn 满足c i t t的最小的卡车号 tsl 电铲有空可以开始装料的时刻 tbu 卸位有空可以开始卸料的时刻 tsf 第w班中电铲空闲时间总和 tbf 第w班中卸位空闲时间总和 tcwl 第w班中待装卡车等待时间总和 tcwu 第w班中待卸卡车等待时间总和 mine 第w班装运矿石总量 tsf total 20个班电铲空闲时间总和 tbf total 20个班卸位空闲时间总和 tcwl total 20个班待装卡车等待时间总和 tcwu total 20个班待卸卡车等待时间总和 mine total 20个班装运矿石总量 初始化数据w 0 电铲空闲 打印输出 否 初始化第w班数据 取t min c i t i 1 m c min i c i t t c c s的值为多少 是 统计电铲空闲时间 统计等待装车时间 设置开始装车时间 计算装完时刻c c t tsl c c s 2 1 计算到卸位时刻c c t c c s 3 2 3 4 统计该班产量 计算到采石场时刻c c t c c s 1 c c t 360 是 否 w w 1 统计前w个班的数据 w 20 是 否 卸位空闲 否 是 统计卸位空闲时间 统计等待卸车时间 设置开始卸车时间 计算卸完时刻c c t tbu c c s 4 表 不同卡车数的模拟结果 图 卡车数量与电铲和卡车利用率的关系 可以看出 电铲效率随卡车数m的增加而上升 而卡车效率却随m的增加而减少 若要求电铲的效率在85 以上 则应保证m 7 此时 为保证卡车有较高效率 75 可取m 8或9 从班产量的角度看 电铲效率越高 即卡车数m越大 班产量相应越高 因此应取m 7 从表中看出 将卡车数从9辆提高到10辆时 班产量只提高了约2 23 却付出了卡车效率降低6 07 的代价 所以最多只考虑9辆卡车的情况 练习 飞机起飞调度问题 1989年美国数学模型竞赛b题 5模拟结果的统计分析 在对一个具有随机因素的系统设计出计算机模拟模型后 不能只将单次模拟运行结果看作问题的真实答案 模拟的输出结果是分布特征未知的随机变量 每次运行的结果仅仅是对该随机变量所有观察值总体的一次抽样 对总体的代表性不够 虽然可以通过增加模拟运行的时间从而增加抽样次数 但这些数据总是由一个 种子 经过一定的算法而获得的伪随机序列 它们是自相关的 并不能构成统计上独立的随机样本 1 系统性能测度 考察系统的性能 就是根据随机样本的取值来估计系统真实参数的统计计量 一般可用点估计和区间估计 令ex dx分别为系统真实参数所对应的随机变量的数学期望及方差 设模拟输出的随机变量为则 分别是系统参数ex和dx的点估计 为了知道ex落在某个范围之内的可能性 并且希望这个范围越小越好 需对ex进行区间估计 当模拟结果构成独立同分布的随机样本时 点估计和分别是ex和dx的无偏估