随机模拟方法及习题

随机模拟方法

在用传统方法难以解决的问题中，某些问题含有不确定的随机因素，分析起来通常比确定性的问题困难。有的模型难做定量分析，得不到解析的结果或者是有解析结果，但计算代价太大以至不能使用，在这种情况下，可以考虑随机模拟的方法即Monte Carlo 方法。该方法是一类以概率统计理论为指导的非常重要的数值计算方法，也是一种用于解决数值问题的基于计算机的统计抽样方法。目前，随机模拟方法已广泛应用于诸如生物信息学、统计物理学、计算机科学、材料科学、金融学和经济学等领域。

基本知识

基本思想

为了求解物理、数学、工程技术以及生产管理等方面的问题，首先建立一个概率或者随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或者抽样实验来计算所求参数的统计特征，最后给出所求解的近似值。而解的精确度可用估计值的标准误差来表示。随机模拟方法是一种独具风格的数值计算方法，其优点大致有如下三方面：（A ）方法的程序结构简单；（B ）算法的概率性和问题的维数无关；（C ）方法的适应强。

随机数和伪随机数

用Monte Carlo 方法模拟某过程的时候，需要产生各种概率分布的随机变量。最基本、最简单、最重要的随机变量是在[0,1]上均匀分布的随机变量。为了方便，通常把[0,1]上均匀分布随机变量的抽样值称为随机数，其他分布随机变量的抽样都可以借助于随机数来实现，因此，随机数是随机抽样的基本工具。在计算机上用数学的方法产生随机数是目前广泛使用的方法，它的特点是占用内存少、产生速度快、又便于重复产生，比如说平方取中法、移位指令加法、同余法等等。然而这种随机数是根据确定的递推公式求得的，存在着周期现象，初值确定后所有随机的数便被唯一确定下来，不满足真正随机数的要求，所以通常称数学方法产生的随机数为伪随机数。在实际应用中，只要这些伪随机数序列通过一系列的统计检验，还是可以把它当称“真正”的随机数来使用。

产生随机数的命令

在Matlab 软件中，可以直接产生满足各种分布的随机数，相关命令如下：

(1)产生m n 阶[,]a b 均匀分布(,)U a b 的随机数矩阵：unifrnd (a,b,m, n)；

产生一个[,]a b 均匀分布的随机数：unifrnd (a,b)；

(2)产生m n ?阶[0,1]均匀分布的随机数矩阵：rand (m, n)； 产生一个[0,1]均匀分布的随机数：rand ；

(3) 产生m n ?阶均值为μ，方差为2

σ的正态分布的随机数矩阵：normrnd (,μσ,m,

n)；

(4) 产生m n ?阶期望值为μ的指数分布的随机数矩阵：exprnd(μ,m,n)

若连续型随机变量X 的概率密度函数为

0()00

x

e x

f x x λλ-?≥=?

其中0λ>为常数，则称X 服从参数为λ的指数分布。排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。

注：参数为λ的指数分布的期望值为1

λ，故在MATLAB 中产生参数为λ的指数分布的命令为：exprnd(1/λ)

（5）产生m n ?阶参数为λ的泊松分布的随机数矩阵：poissrnd(λ,m,n) 注：参数为λ的泊松分布的期望值为λ。 应用示例

例一：排队论中的模拟

排队论主要研究随机服务系统的工作过程。在排队系统中，服务对象的到达时

间和服务时间都是随机的。排队论通过对每个个别的随机服务现象的统计研究，找出反映这些随机现象平均特性的规律，从而为设计新的服务系统和改进现有服务系统的工作提供依据。对于排队服务系统, 顾客常常注意排队的人是否太多, 等候的时间是否长, 而服务员则关心他空闲的时间是否太短。 于是人们常用排队的长度、等待的时间及服务利用率等指标来衡量系统的性能。

单服务员的排队模型：

在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客。当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店。

设：

(1)顾客到来间隔时间服从参数为0.1的指数分布；

(2)对顾客的服务时间服从[4,15]上的均匀分布； (3)排队按先到先服务规则，队长无限制。

假定一个工作日为8小时，时间以分钟为单位。要求：

[1] 模拟一个工作日内完成服务的个数及顾客平均等待时间t 。

[2] 模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间。 假设：

（1）顾客源是无穷的； （2）排队的长度没有限制；

（3）到达系统的顾客按先后顺序依次进入服务， 即“先到先服务”。 符号说明

w ：总等待时间； i c ：第i 个顾客的到达时刻； i

b ：第i 个顾客开始服务时刻；

i

e ：第i 个顾客服务结束时刻；

i x ：第i 个顾客与第1i +个顾客之间到达的间隔时间；

i

y ：对第i 个顾客的服务时间。

显然，如下关系式成立：

1i i i

c c x -=+；

i i i

e b y =+；

1m ax(,)

i i i b c e -=

模拟框图：

图8-1 程序框图

求解

经过编写程序计算，模拟出一个工作日内完成服务的顾客数34m =人，平均等待时间约42t =分钟。注意每次运行的结果会不一样，甚至结果的偏差较大。 类似地，模拟出100个工作日内完成服务的顾客数43m =人，平均等待时间约27t =分钟。注意每次运行的结果会不一样，但是结果的偏差较小，总的来说，数值结果

比较稳定。

例二：零件的参数设计

一件产品由若干零件组装而成，标志产品性能的某个参数取决于这些零件的参数。零件参数包括标定值和容差两部分。进行成批生产时，标定值表示一批零件该参数的平均值，容差则给出了参数偏离标定值的容许范围。若将零件参数视为随机变量，则标定值代表期望值，在生产部门无特殊要求时，容差通常规定为均方差的3倍。 进行零件参数设计，就是要确定其标定值和容差。这时要考虑两个方面因素： 当各零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大。

零件容差的大小决定了其制造成本，容差设计得越小，成本越高。 试通过如下的具体问题给出一般的零件参数设计方法。

粒子分离器某参数（记作y ）由7个零件的参数（记作127,,,x x x ）决定，经验公式为：

7

616

.1242

3

56

.02485

.012

35136.0162.2142.174x x x x x x x x x x x y ???

?

?????

?

???

??

???

??--?

???

??-???? ??=-．

y

的目标值（记作

y ）为1.50。当y 偏离0y

± 0.1时，产品为次品，质量损失为1000

元；当当y 偏离

y ± 0.3时，产品为废品，质量损失为9000元。

零件参数的标定值有一定的容许变化范围；容差分为A 、B 、C 三个等级，用与标定值的相对值表示，A 等为 ±1%，B 等为 ±5%，C 等为 ±10%。7个零件参数标定值的容许范围及不同容差等级零件的成本（元）如下表（符号 / 表示无此等级零件）：

表8-1 零件参数的相关数据

现成批生产，每批产量1000个。在原设计中，7个零件参数的标定值为：1x = 0.1，

2

x = 0.3，3x = 0.1，4x = 0.1，5x = 1.5，6x = 16，7x = 0.75；容差均取最便宜的等级。 请你综合考虑y 偏离0y 造成的损失和零件成本，重新设计零件参数（包括标定值

和容差），并与原设计比较，总费用降低了多少。

（1） 关于零件参数的假设

由于零件在加工制造过程中存在许多随机因素，因此，由中心极限定理知道零件的参数可以看成是服从正态分布的随机变量。设七个零件的加工是独立的，则七个零件参数可以视作是相互独立的正态随机变量。即

2

0(,),

1,2,...,7i i i x N x i σ= ，

其中，

i x 是第i 个零件参数i x 的标定值，也就是随机变量i x

的均值。

另一方面，零件参数的容差的选择均取最便宜的等级（C 等：r =10%），又根据概率统计中的3σ原则，即容差通常是均方差的3倍，因此，

容差=03i i i x r σ=，

所以，

0/3,1,2,...,7

i i i x r i σ==。

（2） 产品参数y 的分布

为求产品质量的损失费用，必须求出产品的次品率和废品率，因而必须求出

产品参数y

的分布，设y 由七个零件的参数表示的经验公式为：

127(,,...,)y f x x x =，

从理论上，利用

i

x 的分布函数以及函数f 的形式，可推导出y 的近似分布。此处，采

用随机模拟的方法来假设、验证y 的分布。

步骤一：利用计算机命令产生n 组相互独立的正态分布随机数

127(,,...,)

j j j j x x x x =， 1,2,...,j n =。

步骤二：由

127(,,...,)

j j j j y f x x x =得到y 的一组样本值，画出频率直方图（如图），

用分布拟合的2

χ检验法检验y 服从正态分布的假设。

图8-2 产品参数y 的直方图

（3） 模型的建立

显然，目标函数应该是产品的总费用最小，而产品总费用=零件总成本+质量损失费用，因为问题带有随机因素，所以费用应该选择期望值。 引入

ij c ：第i 个参数取第j 个容差等级时所需的成本，1,2,...,7i =，1,2,3j =。 ij

d ：0-1变量，如果第i 个参数取第j 个容差等级则取1，否则取0。

根据（2）的分析结果，2

(,)y y N y

σ ，概率密度函数2

2()2()y y y f y σ

--=

，利用

一般正态分布和标准正态分布的关系，所以

正品的概率

1.50.1

1 1.50.1

()p f y dy

+-=

?

,

类似地，次品的概率

1.4

1.8

2 1.2

1.6

()()p f y dy f y dy =

+?

?

1.4 1.2 1.8 1.6(

)(

)(

)(

)

y

y

y

y

y

y

y

y

σσσσ----=Φ-Φ+Φ-Φ ,

()

Φ 表示标准正态分布的分布函数。 废品的概率

1.23 1.8

()()p f y dy f y dy

∞

-∞

=

+

?

?

1.8 1.21(

)(

)

y

y

y

y

σσ--=-Φ+Φ .

因此目标函数为

7

3

231

1

1000(10009000)

ij

ij i j c

d p p ==?+?+?∑

∑.

问题就是要使目标函数最小的零件参数的标定值0i x 和容差等级ij

d ，则问题的数学模

型即为

m in

7

3

231

1

1000(10009000)

ij

ij i j c

d p p ==?+?+?∑

∑

..s t

i i i

a x

b ≤≤，

3

1

1

ij

j d

==∑，

0,1,

,1,2,...,7

ij d i ==，

其中

,i i

a b 是零件参数标定值的下界和上界，在问题中已经给出。

（4） 模型的求解

上述模型有较多的算法求解，在这我们仍然采用随机模拟的方法来求解。 步骤一：枚举1081233132=??????种容差等级组合； 步骤二：在每种组合下，产生

[,]

i i a b 上均匀分布的随机数i x

，并形成一个初始的可行

点

127(,,...,)

i i i x x x x =，1,2,...,i N =

步骤三：将每一个可行点代入目标函数，进行比较，并输出目标函数值最小的对应

可行点及容差等级组合。 具体结果为： 最优标定值

*

[0.075

0.3750.11950.1198 1.20817.01590.5625]

x =,

最优容差等级

*

[]

d B

B B

C C B B =,

目标函数值为42.14万元。

习题

1．某报童以每份0.03元的价格买进报纸，以0.05元的价格出售，根据长期统计，报纸每天的销售量及百分率为

已知当天销售不出去的报纸，将以每份0.02元的价格退还报社，试用模拟的方法确定报童每天买进报纸的数量，使报童的平均总收入最大？

参考资料

[1] 徐钟济. 蒙特卡洛方法[M]，上海：上海科学技术出版社，1985

[2] 赵静，但琦.数学建模和数学实验[M]，北京：高等教育出版社，2003

[3] 朱道元等. 数学建模案例精选[M]，北京：科学出版社，2003