1、随机模拟方法在用传统方法难以解决的问题中，某些问题含有不确定的随机因素，分析起来通常比确定性的问题困难。有的模型难做定量分析，得不到解析的结果或者是有解析结果，但计算代价太大以至不能使用，在这种情况下，可以考虑随机模拟的方法即 Monte Carlo 方法。该方法是一类以概率统计理论为指导的非常重要的数值计算方法，也是一种用于解决数值问题的基于计算机的统计抽样方法。目前，随机模拟方法已广泛应用于诸如生物信息学、统计物理学、计算机科学、材料科学、金融学和经济学等领域。基本知识基本思想为了求解物理、数学、工程技术以及生产管理等方面的问题，首先建立一个概率或者随机过程，使它的参数等于问题的解；然后通
2、过对模型或过程的观察或者抽样实验来计算所求参数的统计特征，最后给出所求解的近似值。而解的精确度可用估计值的标准误差来表示。随机模拟方法是一种独具风格的数值计算方法，其优点大致有如下三方面：（A）方法的程序结构简单；（ B）算法的概率性和问题的维数无关；（C ）方法的适应强。随机数和伪随机数用 Monte Carlo 方法模拟某过程的时候，需要产生各种概率分布的随机变量。最基本、最简单、最重要的随机变量是在 上均匀分布的随机变量。为了方便，通常0,1把 上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都可以0,1借助于随机数来实现，因此，随机数是随机抽样的基本工具。在计算机上用数学的
3、方法产生随机数是目前广泛使用的方法，它的特点是占用内存少、产生速度快、又便于重复产生，比如说平方取中法、移位指令加法、同余法等等。然而这种随机数是根据确定的递推公式求得的，存在着周期现象，初值确定后所有随机的数便被唯一确定下来，不满足真正随机数的要求，所以通常称数学方法产生的随机数为伪随机数。在实际应用中，只要这些伪随机数序列通过一系列的统计检验，还是可以把它当称“真正”的随机数来使用。产生随机数的命令在 Matlab 软件中，可以直接产生满足各种分布的随机数，相关命令如下：(1)产生 阶 均匀分布 的随机数矩阵：unifrnd (a,b,m, n)；mn,ab(,)Uab产生一个 均匀分布的
4、随机数：unifrnd (a,b)；,ab(2)产生 阶 均匀分布的随机数矩阵：rand (m, n)；mn01产生一个 均匀分布的随机数：rand；,(3) 产生 阶均值为 ，方差为 的正态分布的随机数矩阵： normrnd ( ,m, 2,n)；(4) 产生 阶期望值为 的指数分布的随机数矩阵：exprnd( ,m,n)mn若连续型随机变量 的概率密度函数为 ,X0()xef其中 为常数，则称 服从参数为 的指数分布。排队服务系统中顾客到0达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。注：参数为 的指数分布的期望值为 ，故在 MA
5、TLAB 中产生参数为 的指数分1布的命令为：exprnd( )1/（5）产生 阶参数为 的泊松分布的随机数矩阵：poissrnd( ,m,n)mn注：参数为 的泊松分布的期望值为 。应用示例例一：排队论中的模拟排队论主要研究随机服务系统的工作过程。在排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对每个个别的随机服务现象的统计研究，找出反映这些随机现象平均特性的规律，从而为设计新的服务系统和改进现有服务系统的工作提供依据。对于排队服务系统, 顾客常常注意排队的人是否太多, 等候的时间是否长, 而服务员则关心他空闲的时间是否太短。 于是人们常用排队的长度、等待的时间及服务利用率等
6、指标来衡量系统的性能。单服务员的排队模型：在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客。当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店。设：(1)顾客到来间隔时间服从参数为 0.1 的指数分布；(2)对顾客的服务时间服从4,15 上的均匀分布；(3)排队按先到先服务规则，队长无限制。假定一个工作日为 8 小时，时间以分钟为单位。要求：1 模拟一个工作日内完成服务的个数及顾客平均等待时间 。t2 模拟 100 个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间。假设：（1）顾客源是无穷的；（2）排队的长度没有限制；（3）到达系统的顾客按先后顺序依次进入
7、服务， 即“先到先服务” 。符号说明 ：总等待时间； ：第 个顾客的到达时刻；wic：第 个顾客开始服务时刻； ：第 个顾客服务结束时刻；ibie：第 个顾客与第 个顾客之间到达的间隔时间；ix1i：对第 个顾客的服务时间。iyi显然，如下关系式成立：； ；1iicxiieby1max(,)iice2e12345c2b11 5b4e3e模拟框图：4b3b初始化：令 1,0,iew产生间隔时间随机数 参数为 0.1 的指数分布ix:,iiicxb产生服务时间随机数 4, 15的均匀分布iyiiey累计等待时间： iwbc准备下一次服务： 1i产生间隔时间随机数 参数为 0.1 的指数分布ix:1
8、iicx确定开始服务时间： 1max(,)iibce480?ib1;/itwi输出结果：完成服务个数：mi平均等待时间： w停 止Y N图 8-1 程序框图求解经过编写程序计算，模拟出一个工作日内完成服务的顾客数 人，平均34m等待时间约 分钟。注意每次运行的结果会不一样，甚至结果的偏差较大。42t类似地，模拟出 100 个工作日内完成服务的顾客数 人，平均等待时间约43分钟。注意每次运行的结果会不一样，但是结果的偏差较小，总的来说，数27t值结果比较稳定。例二：零件的参数设计一件产品由若干零件组装而成，标志产品性能的某个参数取决于这些零件的参数。零件参数包括标定值和容差两部分。进行成批生产时
9、，标定值表示一批零件该参数的平均值，容差则给出了参数偏离标定值的容许范围。若将零件参数视为随机变量，则标定值代表期望值，在生产部门无特殊要求时，容差通常规定为均方差的 3 倍。进行零件参数设计，就是要确定其标定值和容差。这时要考虑两个方面因素：当各零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大。零件容差的大小决定了其制造成本，容差设计得越小，成本越高。试通过如下的具体问题给出一般的零件参数设计方法。粒子分离器某参数（记作 ）由 7 个零件的参数（记作 ）决定，经验公y127,x式为：7616.24356.02485.012351 .1.24.7 xxy
10、 的目标值（记作 ）为 1.50。当 偏离 0.1 时，产品为次品，质量损失为0y01000 元；当当 偏离 0.3 时，产品为废品，质量损失为 9000 元。y零件参数的标定值有一定的容许变化范围；容差分为 A、B 、C 三个等级，用与标定值的相对值表示，A 等为 1%，B 等为 5%，C 等为 10%。7 个零件参数标定值的容许范围及不同容差等级零件的成本（元）如下表（符号 / 表示无此等级零件）：表 8-1 零件参数的相关数据标定值容许范围 C 等 B 等 A 等1x0.075, 0.125 / 25 /20.225, 0.375 20 50 /3x0.075, 0.125 20 50
11、20040.075, 0.125 50 100 5005x1.125, 1.875 50 / /612, 20 10 25 1007x0.5625, 0.935 / 25 100现成批生产，每批产量 1000 个。在原设计中，7 个零件参数的标定值为： = 1x0.1， = 0.3， = 0.1， = 0.1， = 1.5， = 16， = 0.75；容差均取最便宜的2x34x56x7等级。请你综合考虑 偏离 造成的损失和零件成本，重新设计零件参数（包括标定y0值和容差） ，并与原设计比较，总费用降低了多少。（1） 关于零件参数的假设由于零件在加工制造过程中存在许多随机因素，因此，由中心极限定
12、理知道零件的参数可以看成是服从正态分布的随机变量。设七个零件的加工是独立的，则七个零件参数可以视作是相互独立的正态随机变量。即，20(,)1,.7iixNi:其中， 是第 个零件参数 的标定值，也就是随机变量 的均值。0ixi ix另一方面，零件参数的容差的选择均取最便宜的等级（C 等： 10%） ，又根r据概率统计中的 原则，即容差通常是均方差的 3 倍，因此，3容差= ，0iixr所以， 。0/,12,.7iixr（2） 产品参数 的分布y为求产品质量的损失费用，必须求出产品的次品率和废品率，因而必须求出产品参数 的分布，设 由七个零件的参数表示的经验公式为：，127(,.)yfx从理论上
13、，利用 的分布函数以及函数 的形式，可推导出 的近似分布。此处，ix y采用随机模拟的方法来假设、验证 的分布。步骤一：利用计算机命令产生 组相互独立的正态分布随机数n， 。127(,.)jjjxx1,2.步骤二：由 得到 的一组样本值，画出频率直方图（如图） ，7()jjjyfxy用分布拟合的 检验法检验 服从正态分布的假设。2y1.4 1.6 1.8 2 2.2 2.4 2.6050100150200250300350图 8-2 产品参数 的直方图y（3） 模型的建立显然，目标函数应该是产品的总费用最小，而产品总费用=零件总成本+质量损失费用，因为问题带有随机因素，所以费用应该选择期望值。
14、引入：第 个参数取第 个容差等级时所需的成本， ， 。ijcj 1,2.7i1,23j：0-1 变量，如果第 个参数取第 个容差等级则取 1，否则取 0。ijdij根据（2）的分析结果， ，概率密度函数 ，利2(,)yN: 2()()yfye用一般正态分布和标准正态分布的关系，所以正品的概率,1.50()pfyd类似地，次品的概率 1.41.826()()pfydf,.21.8.6()()yyyy表示标准正态分布的分布函数。():废品的概率 1.231.8()()pfydfy.1.8.2()()yy因此目标函数为.732310(090)ijicdp问题就是要使目标函数最小的零件参数的标定值 和
15、容差等级 ，则问题的数学0ixijd模型即为 min732310(09)ijicdp，.stiiaxb， ，31ijd0,1,2,.7iji其中 是零件参数标定值的下界和上界，在问题中已经给出。,iab（4） 模型的求解上述模型有较多的算法求解，在这我们仍然采用随机模拟的方法来求解。步骤一：枚举 种容差等级组合；1082312步骤二：在每种组合下，产生 上均匀分布的随机数 ，并形成一个初始的可,iabix行点 ，127(.,)iixx1,2.N步骤三：将每一个可行点代入目标函数，进行比较，并输出目标函数值最小的对应可行点及容差等级组合。具体结果为：最优标定值,*0.75.30.195.81.2
16、07.590.62x最优容差等级,*dBCB目标函数值为 42.14 万元。习 题1 某报童以每份 0.03 元的价格买进报纸，以 0.05 元的价格出售，根据长期统计，报纸每天的销售量及百分率为已知当天销售不出去的报纸，将以每份 0.02 元的价格退还报社，试用模拟的方法确定报童每天买进报纸的数量，使报童的平均总收入最大？参考资料 1 徐钟济. 蒙特卡洛方法M ，上海：上海科学技术出版社， 19852 赵静，但琦. 数学建模和数学实验M ，北京：高等教育出版社， 20033 朱道元等. 数学建模案例精选M ，北京：科学出版社， 2003销售量 200 210 220 230 240 250百分率 0.10 0.20 0.40 0.15 0.10 0.05
