计算机模拟版本3.ppt

数学建模培训 主讲 何仁斌 概述 计算机模拟是用计算机对系统的结构和行为进行动态演示 以评价或预测系统的行为效果 为决策者提供信息的一种方法 它是解决较复杂的实际问题的一条有效途径 计算机模拟也可以说是用计算机程序直接建立真实系统的模型 并通过计算了解系统随时间变化的行为或特征 应用领域 航空 机电 冶金 社会经济 交通运输 生态系统等 计算机模拟分为连续系统模拟和离散系统模拟 状态随时间连续变化的系统称为连续系统 通常该系统的模型一般可以用微分方程的形式表达 通过一些物理机理推导出来 模拟结果往往是近似的 例如 追逐问题 浓度问题 一 连续系统 1 追逐问题 A BB CC DD A 状态变量A x t y t 假设 v 1m s 追逐的方向是对准追逐目标 A 0 10 B 10 10 D 0 0 C 0 10 试确定四人追逐的运动轨迹 A x t t y t t 模拟方法 1 建立平面直角坐标系 2 以时间间隔 t进行采样 并且计算各个时刻下的状态 A x1 t y1 t x1 t t y1 t t B x2 t y2 t x2 t t y2 t t x1 t t y1 t t 由几何原理 x1 t v tcos y1 t v tsin 在t时刻下 3 选取足够小的 t 模拟到任意两人的距离小于v t为止 4 连接四人在各时刻下的位置 就得到所求的运动轨迹 Matlab simu2 m v 1 dt 0 05 d 20 x x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 000101010100 x 9 x 1 x 10 x 2 holdonaxis equal axis 010010 while d 0 1 fori 1 2 7d sqrt x i x i 2 2 x i 1 x i 3 2 x i x i v dt x i 2 x i d x i 1 x i 1 v dt x i 3 x i 1 d plot x i x i 1 endx 9 x 1 x 10 x 2 endholdon A B C D 某军的一导弹基地发现正北方向120km处海面上有敌艇一艘以90km h的速度向正东方向行驶 该基地立即发射导弹跟踪追击敌艇 导弹速度为450km h 自动导航系统使导弹在任一时刻都能对准敌艇 试问导弹在何时何处击中敌艇 2 导弹跟踪问题 P x t y t O x M vet H y A 0 H 正东方向 正北方向 H 120 千米 km h km h 敌导 xk t yk t simu1 mstest1 mhuatu m 微分方程建模 计算机模拟 消去t 微分方程建模 1 2 x dsolve D2x 120 y sqrt Dx 2 1 90 450 Dx 0 0 x 0 0 y 微分方程求解 simplify x ans 45 120 1 5 y 120 4 5 1800 1 4 5 1 3 5 202500 1 5 y 120 6 5 72 1 3 5 1 2 45 120 1 5 y 120 4 5 1800 1 4 5 1 3 5 202500 1 5 y 120 6 5 72 1 3 5 1 2 y 120 eval x ans 25 0000 0 0000i 25 0000 0 0000i 如何决定导弹位置P2 x2 y2 计算机模拟 1 当t 0时 导弹位置O 0 0 敌艇位置A 0 120 2 当t 时 导弹位置P1 x1 y1 敌艇位置M1 90 120 3 当t 2 时 导弹位置P2 x2 y2 敌艇位置M2 90 2 120 P x t y t O x M vet H y A 0 H 正东方向 正北方向 H 120 千米 km h km h 敌导 xk t yk t function num y j L T simu1 x y t eps k 0 whilek 1000 simu1 m p 90 k t x q 120 y d1 p p 2 q 2 0 5 d2 q p 2 q 2 0 5 x x 450 t d1 y y 450 t d2 if abs q eps break endk k 1 endnum k L x y j y T L 90 stest1 m t 0 1 x 0 y 0 num y j L T simu1 x y t eps 计算结果 Num 2 y j 131 2155 L 22 675 km T 0 2519 h 导弹的位置 x t y t 能否用图形描述它们的追逐运动过程 敌艇的位置 水平方向 90 k t 120 x 0 0 huatu m axis equal axis 0 140 0 140 gridonholdont 0 001 步长while abs x 2 120 0 1 终止条件 fork 1 280p 90 k t x 1 q 120 x 2 d1 p p 2 q 2 0 5 d2 q p 2 q 2 0 5 x 1 x 1 450 t d1 x 2 x 2 450 t d2 x1 1 90 k t x1 2 120 h1 line color 0 0 2 0 4 linewidth 2 h2 line color 0 0 6 0 9 linewidth 3 set h1 xdata x1 1 ydata x1 2 set h2 xdata x 1 ydata x 2 endendholdon 2 如果当基地发射导弹的同时 敌艇立即由仪器发觉 假定敌艇为高速快艇 它即刻以135千米 小时的速度与导弹方向成垂直的方向逃逸 问导弹何时何地击中敌艇 不停地改变逃跑策略 运动轨迹如何 修正敌艇速度 90 135千米 小时 思考 二 离散系统 离散系统是指系统状态只在有限的时间点或可数的时间点上发生变化的系统 并假设离散系统状态的变化是在一个时间点上瞬间完成 模型一般用流程图或网络来表示 其间可能涉及到随机事件等 关键 模拟步骤 数据收集 模拟时钟 在我方某前沿防守地域 敌人以一个炮排 含两门火炮 为单位对我方进行干扰和破坏 为躲避我方打击 敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现 我方指挥所对敌方目标的指示有50 是准确的 而我方火力单位 在指示正确时 有1 3的射击效果能毁伤敌人一门火炮 有1 6的射击效果能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来 确定有效射击的比率及毁伤敌方火炮的平均值 分析 这是一个概率问题 可以通过理论计算得到相应的概率和期望值 但这样只能给出作战行动的最终静态结果 而显示不出作战行动的动态过程 为了能显示我方20次射击的过程 现采用模拟的方式 射击问题 需要模拟出以下两件事 1 问题分析 2 当指示正确时 我方火力单位的射击结果情况 1 观察所对目标的指示正确与否 模拟试验有两种结果 每一种结果出现的概率都是1 2 因此 可用投掷一枚硬币的方式予以确定 当硬币出现正面时为指示正确 反之为不正确 模拟试验有三种结果 毁伤一门火炮的可能性为1 3 即2 6 毁伤两门的可能性为1 6 没能毁伤敌火炮的可能性为1 2 即3 6 这时可用投掷骰子的方法来确定 如果出现的是 三个点 则认为没能击中敌人 如果出现的是 点 则认为毁伤敌人一门火炮 若出现的是 点 则认为毁伤敌人两门火炮 2 符号假设 i 要模拟的打击次数 k1 没击中敌人火炮的射击总数 k2 击中敌人一门火炮的射击总数 k3 击中敌人两门火炮的射击总数 E 有效射击比率 E1 20次射击平均每次毁伤敌人的火炮数 3 模拟框图 4 模拟结果 5 理论计算 6 结果比较 虽然模拟结果与理论计算不完全一致 但它却能更加真实地表达实际战斗动态过程 蒙特卡罗方法 MonteCarlomethod 蒙特卡罗方法 MonteCarlomethod 或称计算机随机模拟方法 是一种基于 随机数 的计算方法 源于美国在第二次世界大战研制原子弹的 曼哈顿计划 该计划的主持人之一数学家冯 诺伊曼用驰名世界的赌城 摩纳哥的MonteCarlo 来命名这种方法 为它蒙上了一层神秘色彩 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用 早在17世纪 人们就知道用事件发生的 频率 来决定事件的 概率 19世纪人们用蒲丰投针的方法来计算圆周率 上世纪40年代电子计算机的出现 特别是近年来高速电子计算机的出现 使得用数学方法在计算机上大量 快速地模拟这样的试验成为可能 蒲丰投针实验近似计算圆周率 蒲丰投针实验 法国科学家蒲丰 Buffon 在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子 蒲丰投针实验的重要性并非是为了求得比其它方法更精确的 值 而是它开创了使用随机数处理确定性数学问题的先河 是用偶然性方法去解决确定性计算的前导 由此可以领略到从 概率土壤 上开出的一朵瑰丽的鲜花 蒙特卡罗方法 MC 蒲丰投针实验可归结为下面的数学问题 平面上画有距离为a的一些平行线 向平面上任意投一根长为l l a 的针 假设针落在任意位置的可能性相同 试求针与平行线相交的概率P 从而求 蒲丰投针实验近似计算圆周率 蒲丰投针实验 如右图所示 以M表示针落下后的中点 以x表示M到最近一条平行线的距离 以 表示针与此线的交角 针落地的所有可能结果满足 其样本空间视作矩形区域 面积是 针与平行线相交的条件 它是样本空间 子集A 面积是 symslphi int l 2 sin phi phi 0 pi ans l因此 针与平行线相交的概率为 从而有 特别当时 蒲丰投针实验近似计算圆周率 蒲丰投针实验近似计算圆周率 蒲丰投针实验的计算机模拟 formatlong 设置15位显示精度a 1 l 0 6 两平行线间的宽度和针长figure axis 0 pi 0 a 2 初始化绘图板set gca nextplot add 初始化绘图方式为叠加counter 0 n 2010 初始化计数器和设定投针次数x unifrnd 0 a 2 1 n phi unifrnd 0 pi 1 n 样本空间 fori 1 nifx i l sin phi i 2 满足此条件表示针与线的相交plot phi i x i r frame i getframe 描点并取帧counter counter 1 统计针与线相交的次数endendfren counter n pihat 2 l a fren 用频率近似计算 蒲丰投针实验近似计算圆周率 蒲丰投针实验的计算机模拟 意大利数学家拉泽里尼得到了准确到6位小数的 值 不过他的实验因为太准确而受到了质疑 蒲丰投针实验计算圆周率 蒙特卡罗投点法是蒲丰投针实验的推广 在一个边长为a的正方形内随机投点 该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值 即n 10000 a 2 m 0 fori 1 nx rand 1 a y rand 1 a if x a 2 2 y a 2 2 a 2 2 m m 1 endenddisp 近似为 num2str 4 m n 以下为向量化程序x rand 1 10000 y rand 1 10000 num2str length find x 1 2 y 1 2 1 4 10000 一些MonteCarlo思想的例子 例1 相遇问题 甲 乙两船在24小时内独立地随机到达码头 设两船到达码头时刻都服从 0 24 上的均匀分布 如果甲船到达码头后停留2小时 乙船到达码头后停留1小时 问两船相遇的概率有多大 可用几何概率 此处略 分析 如果知道甲 乙两船到达的时刻x和y 两船能相遇的条件就是 两船到达的时刻x和y可以随机生成 生成一组到达时刻 可以确定是否能相遇 如果重复很多次 统计相遇的比例就可近似为相遇的概率 这也是MonteCarlo思想 两船到达码头时刻服从 0 24 上的均匀分布 甲船停留2小时 乙船停留1小时 相遇概率 方法一 利用软件产生一组 0 24 上的随机数xi 10 45 30 39 020 16 87 8 5 4分别代表甲船到达码头时间 再产生一组 0 24 上的随机数yi 6 03 48 117 50 813 4 14 0分别代表乙船到达码头时间 把 xi yi i 1 2 当作一天中甲船乙船到达的时间 统计出能相遇的频数 计算出相遇的频率做为相遇的概率 可以使用软件 excel C等实现 两船到达码头时刻服从 0 24 上的均匀分布 甲船停留2小时 乙船停留1小时 相遇概率 方法二 clc N 1000 c 0 模拟次数N 相遇次数c清零fori 1 N 重复N次到达时间x unifrnd 0 24 1 1 甲船到达时间x 随机数 y unifrnd 0 24 1 1 乙船到达时间y 随机数 if x y 如果能相遇 则计数器加1endendP c N 显示相遇的概率近似值 注意 每次运行的结果一般都不一样 计算机模拟程序 报童诀窍的简化版 报童每天清晨从报社购进报纸零售 晚上退回没有卖掉的报纸 若每份报纸的购进价为b 0 75元 售出价为a 1元 退回价为c 0 6元 每天需求量X是离散型随机变量 其分布为 例2 问 如果报童每天购进报纸为n 510份 每天的平均利润是多少 方法二 如果我们知道每天的需求量 可直接计算利润 而每天需求量可以按分布生成 随机模拟思想 方法一 概率方法 略 售出价a 1 购进价b 0 75 退回价c 0 6 购进数量n 510份 按照需求量的分布规律 随机生成N 20个数据 510520500510520500500500500510500500500520510520510510代表20天的需求量 计算出报童在这20天的总利润和平均利润 用平均利润来近似报童的平均收入 这也是MonteCarlo方法 问题 如何按分布规律产生随机数据 随机数据很多时 如何编程 报童诀窍的简化版 问题1 如何产生以下分布规律的随机数据 注 rand m n 可以生成 0 1 上均匀分布随机数 把 0 1 分成长度为0 34 0 36 0 30的三个区间 0 0 34 0 34 0 70 0 70 1 用rand 1 1 产生1个 0 1 上均匀分布随机数 如该数在 0 0 34 0 34 0 70 或 0 70 1 内 相当于该天的需求量相应为500 510和520重复多次就可以若干天的需求量了 生成N 20天的需求量的matlab代码可以为 报童诀窍的简化版 生成N 20天的需求量的matlab代码 定义函数 functiony randfun1 N y zeros 1 N fori 1 Nt rand 1 1 ift 0 34X 500 elseift 0 70X 510 elseX 520 endy i X end 产生1行N列的全0矩阵 目的分配好矩阵大小 可以省略 文件名为randfun1 m 报童诀窍的简化版 问题2 如何从需求量计算利润 售出价a 1 购进价b 0 75 退回价c 0 6 购进数n 510 函数文件名fun2 mfunctiony fun2 x a 1 售出价b 0 75 购进价c 0 6 退回价n 510 购进数量if x n y n a b elsey x a b n x b c end 模拟程序代码N 1000 x randfun1 N y 0 fori 1 Ny y fun2 x i endy N 报童诀窍的简化版 a 1 b 0 75 c 0 6 n 510 n 510 N 1000 y 0 fori 1 Nt rand 1 1 if tn y y n a b elsey y x a b n x b c endEndfprintf 1 平均利润 3f y N 也可完整模拟程序 根据随机数t 计算需求量x值 根据需求量x 计算利润并累加到y中 显示平均利润 售出价a购进价b退回价c购进数n总利润y模拟天数N 报童诀窍的简化版 粒子分离器的某关键参数 记作y 由7个零件的参数 记作x1 x2 x7 决定 均值 标定值标准差 容差等级 标定值 3 关键参数y的目标值是1 50 当偏离为 0 1但未达到 0 3时 产品为次品 损失为1000元 当偏离达到 0 3时 产品为废品 损失为9000元 由于工艺原因 7个零件参数可以看着是正态随机变量 在后面的标定值及容差等级情况下 求产品的平均损失 容差等级A 1 B 5 C 15 零件参数设计 1997A 例3 均值 标定值标准差 容差等级 标定值 3 容差等级A 1 B 5 C 15 零件损失Q是一个随机变量 平均损失就是期望E Q 由于y的表达式很复杂 要想计算y的分布和上述概率很困难 我们必须寻找较为有效的近似方法 模拟 通过产生指定分布的随机数 来代表7个零件的参数值 计算y值 确定损失大小 多做几次可得均值 零件参数 编写计算y值的函数 y funY x1 x2 x3 x4 x5 x6 x7 编写计算损失费函数 Q funQ y D1 0 1 0 3 0 1 0 1 1 5 16 0 75 标定值D2 5 5 5 15 15 5 5 100 容差等级a D1b D1 D2 3x1 normrnd a 1 b 1 1 1 x7 normrnd a 7 b 7 1 1 y funY x1 x2 x3 x4 x5 x6 x7 Q funQ y 产生一组零件参数 相当制作一个产品 重复N 1000次 即生产N个产品 求出损失费用的平均 代码略 按指定分布产生随机数 作为零件参数 计算该产品的关键参数y值和损失的大小 零件期望a零件标准差b 零件参数 用蒙特卡洛法解非线性规划问题 基本假设 试验点的第j个分量xj服从 aj bj 内的均匀分布 符号假设 求解过程 先产生一个随机数作为初始试验点 以后则将上一个试验点的第j个分量随机产生 其它分量不变而产生一新的试验点 这样 每产生一个新试验点只需一个新的随机数分量 当K MAXK或P MAXP时停止迭代 在Matlab软件包中编程 共需三个 文件 randlp m mylp m lpconst m 主程序为randlp m mylp mfunctionz mylp x 目标函数z 2 x 1 2 x 2 2 x 1 x 2 8 x 1 3 x 2 转化为求最小值问题 randlp mfunction sol r1 r2 randlp a b n 随机模拟解非线性规划debug 1 a 0 试验点下界b 10 试验点上界n 1000 试验点个数r1 unifrnd a b n 1 n 1阶的 a b 均匀分布随机数矩阵r2 unifrnd a b n 1 sol r1 1 r2 1 z0 inf fori 1 nx1 r1 i x2 r2 i lpc lpconst x1x2 iflpc 1z mylp x1x2 ifz z0z0 z sol x1x2 endendend ToMatlab randlp 返回 随机数的产生 1 均匀随机数 均匀分布U 0 1 rand 2 产生其他随机数的方法逆变换法 舍选法 近似抽样法等 设概率分布函数F x 是严格单调增的 F的反函数记作F 1 先产生U U 0 1 再取X F 1即为所求 如指数分布 分布函数F x 1 exp 可得 3 非 常见分布随机数如何产生 离散 经验分布函数法 设一组实际数据 将它们分组整理形成频率图或表格 x 0 由均匀随机数R 0 1 决定X的抽样值 当0 R 0 3 抽样值xi 0 3 5 xi 3 5当0 3 R 0 8 抽样值xi 3 5 5 5 xi 5 5当R 0 8 抽样值xi 5 5 8 xi 8 r rand if0 rendy 产生一个随机数 Matlab程序 r rand 10 lisanrnd1 my fori 1 10if0 r i endendy 产生随机变量数组 计算结果y Columns1through75 55 53 55 53 55 55 5Columns8through108 05 55 5 MATLAB的实现 常见分布的随机变量获取抽样值 1 U 0 1 R rand m n 2 U a b R unifrnd a b m n 3 E R exprnd 1 m n 4 N 0 1 R randn m n 5 N mu sigma R normrnd mu sigma m n 6 B N p R binornd N p m n 7 P R poissrnd m n 以上语句都是产生m n的随机矩阵 hist x 产生模拟随机数的计算机命令 2 产生m n阶 均匀分布的随机数矩阵 rand m n 产生一个 均匀分布的随机数 rand 1 产生m n阶 a b 均匀分布U a b 的随机数矩阵 unifrnd a b m n 产生一个 a b 均匀分布的随机数 unifrnd a b 当只知道一个随机变量取值在 a b 内 但不知道 也没理由假设 它在何处取值的概率大 在何处取值的概率小 就只好用U a b 来模拟它 当研究对象视为大量相互独立的随机变量之和 且其中每一种变量对总和的影响都很小时 可以认为该对象服从正态分布 机械加工得到的零件尺寸的偏差 射击命中点与目标的偏差 各种测量误差 人的身高 体重等 都可近似看成服从正态分布 若连续型随机变量X的概率密度函数为其中 0为常数 则称X服从参数为的指数分布 指数分布的期望值为 排队服务系统中顾客到达率为常数时的到达间隔 故障率为常数时零件的寿命都服从指数分布 指数分布在排队论 可靠性分析中有广泛应用 注意 Matlab中 产生参数为的指数分布的命令为exprnd 例顾客到达某商店的间隔时间服从参数为0 1的指数分布 指数分布的均值为1 0 1 10 指两个顾客到达商店的平均间隔时间是10个单位时间 即平均10个单位时间到达1个顾客 顾客到达的间隔时间可用exprnd 10 模拟 设离散型随机变量X的所有可能取值为0 1 2 且取各个值的概率为其中 0为常数 则称X服从参数为的帕松分布 帕松分布在排队系统 产品检验 天文 物理等领域有广泛应用 帕松分布的期望值为 指数分布与帕松分布的关系 1 指两个顾客到达商店的平均间隔时间是10个单位时间 即平均10个单位时间到达1个顾客 2 指一个单位时间内平均到达0 1个顾客 例 1 顾客到达某商店的间隔时间服从参数为0 1的指数分布 2 该商店在单位时间内到达的顾客数服从参数为0 1的帕松分布 排队问题 机械故障等候修理 飞机跑道 日常生活中经常遇到的排队问题 自选商场收款台 医院里病人等候就诊 输入情况 顾客到达时间和服务时间 系统状态 排队等候的顾客数目 队长 L t 服务员是否在工作或服务效率等 简图 关系 系统在什么条件下处于空闲状态 yi xi 1 排队系统中 顾客到达时刻数据如何收集 对每个顾客的服务时间如何 X x1 x2 xn 计算机遵循某种规则进行随机抽样 S s1 x2 xn 排队服务系统的模拟 1 模拟目的2 系统假设与输入3 系统的状态4 初始条件和终止条件5 模拟过程的运行6 系统的性能指标7 与分析模型的比较8 其它排队系统 1 模拟目的 如果知道了顾客到达的时间和服务时间的统计规律 大量实际数据或概率分布 怎样通过模拟得到衡量系统的性能指标 并与分析模型的结果进行比较 2 系统假设与输入 首行作如下基本假设 1 顾客源是无穷的 2 排队的长度没有限制 3 到达系统的顾客按先后顺序进入服务 假设由概率分布产生了如下的数据 其中ik表示第k位顾客与第k 1位顾客到达的时间间隔 sk是第k位顾客的服务时间 下面就以这组数据作为系统的输入 3 系统的状态 队长 排队等候的顾客数目 记作L t 服务员的状态用S t 表示 当服务员工作时 令S t 1 空闲时令S t 0 事件 引起系统状态L t 和S t 改变的行为 顾客到达 和 顾客离开 设第k位顾客到达和离开时刻分别为ak和dk 则有 4 初始条件和终止条件 不妨假设t 0是第1位顾客到达的时刻 此时服务员处于空闲状态 模拟的终止时刻是T 给定值 或模拟直到第n个 给定值 顾客进入服务的时间 5 模拟过程 进行模拟时 可以事先产生一批数据ik sk 再计算ak dk 然后让时钟t按ak dk从小到大的顺序跳动 也可在时钟t跳到某一事件发生时 才产生一个所需要的i或s 这里采用后一种方法 设当前时钟为t 队长L t 记作WL 服务员状态S t 记作SS t以后下一个 顾客到达 事件的发生时刻记作AT t以后下一个 顾客离开 事件的发生时刻记作DT AT DT WL 0 SS 0 AT 0 DT 999 置t AT SS 0 置t DT WL 0 WL WL 1 令SS 1 产生s 置DT t s 产生i 置AT t i t T 令SS 0 WL WL 1 产生s 置DT t s 停止 置DT 999 是 否 是 否 是 否 否 是 st 6 it 5 WL 0 SS 0 AT 0 DT 99999 T 10 t 0 ttt 0 www 0 sss 0 while t0WL WL 1 s exprnd 1 st 1 1 DT t s elseSS 0 DT 99999 endendttt ttt t www www WL sss sss SS endttt 2 length ttt www 2 length ttt sss 2 length ttt 6 系统的性能指标 平均队长平均待时间服务利用率 7 与分析模型的比较 对于一般的排队服务系统 很难建立分析模型并求出结果 即使是单服务员系统 也只有在顾客到达间隔和服务时间服从某种特殊的概率分布时 才可以在稳态情况下由分析模型得到上述三个指标的解析形式的结果 分析模型基本假设 1 到达间隔服从参数 的指数分布 2 服务时间服从参数为 的指数分布 3 排队规则是 先到先服务 顾客数目和排队长度都无限制 排队论中将这样的单服务员系统记作M M 1 服务强度 稳态 1 8 其它的排队服务系统 有两个 或多个 服务员串联服务系统 每个顾客接受多次服务从排队顾客中选服务时间最短的最先服务队的长度有限制带有优化目标的排队服务系统 库存太多造成浪费或资金积压 库存太少不能满足需求也造成损失 需要进行决策 何时进货 进多少 使得平均费用最少 而收益最大 库存系统 方案甲 按前一天的销售量作为当天的生产量 方案乙 按前二天的平均销售量作为当天的生产量 假定市场对该产品的每天需要量是一个随机变量 但从以往的统计分析得知它服从正态分布N 135 22 42 例1某企业生产易变质的产品 如蛋糕 当天生产的产品必须当天售出 否则就会变质 该产品单位成本为2 5元 单位产品售价为5元 企业为避免存货过多而造成损失 拟从以下两种货存方案中选出一个较优的方案 如何决定当天的生产数量 模拟的基本思路 需要获得市场对该产品的需要量的样本值 产生抽样值 寻找货存量与需要量之间的关系 按照两种不同方案 货存量 需要量 计算出经T天后企业的利润值 累计 比较大小 从中选出一个较优的方案 MATLAB编制的程序 kucun m 运行程序 kutest m初始状态 T 100 S1 136 S21 126 S22 148 TL1 TL2 kucun T S1 S21 S22 观察计算结果如果多数情况有 TL1 TL2 则认为方案乙较优 思考 1 一般模拟结果波动性较大 如何减少这种波动 2 修改上述程序 可靠性问题 一设备上有三个相同的轴承 每个轴承正常工作寿命为随机变量 概率分布如下 有轴承损坏 设备停止工作 检修工准备开始更换部件 称为一个延迟时间 它也是随机变量 分布如下 主要费用 1 设备停工损失费 5元 分钟 2 检修工人的工时费 12元 小时 3 轴承的成本费 16元 个 更换轴承所需要的时间 一个两个三个203040 min 问题 现在有两种方案 方案一 损坏一个轴承只更换一个轴承 方案二 一旦有轴承损坏就全部更换 试通过计算机模拟对以上两种方案做出评价 随机数怎样产生 模拟时选用时间步长法还是事件步长法 关于随机数的产生见Lrnd m 零件寿命 Yrnd m 延迟时间 方案一的数学模型 kekao1 m目标函数minc Ui T 其中 ti表示延迟时间 设事件发生时刻Ti 它是由轴承工作寿命Li 延迟时间和修理时间迭加产生 总的运行时间为T 10000 例 手工模拟 方案一的初始事件表 第一次刷新后的事件表 c1 5 c2 12 c3 16 kekao1 m g1 20 60 g2 30 60 g3 40 60 Lrnd RandomfunctionLrndgetsarandomlife three yrnd Randomfunctionyrndgetsarandomde