计算机模拟法相关知识
一．什么是“计算机模拟”？ 一个事例题：开张营业时，待售自行车115辆；(售价-进货价-营业税)=10元/辆；
保管费：0.8元/天辆；发生缺货(有顾客要购车而无货)时，损失费：2元/天辆；
顾客对车的需求量：x辆/天，x是从0至99的均匀分布随机数；
想要进货时，必须提前3天发订货单，订货费75元/次，如：2007年07月25日发出订货单，则2007年07月28日一上班就到货；
计划：剩L辆时就立即发订货单，订货量M辆/次，请确定L=？M=？ 题毕。
分析构造模拟模型：模拟运行一年365天，令L、M取各种可能的数据，分别计算总利润，使得利润达到最大的数据L、M即为所求。
初始状态
——————————|
|增加一天
(1),这天初：库存C=？，进货E=？
(2),这天当中：售出y=？，发订货单否？
缺货量Q=？
(3),这天末算帐：当天的净收入F=？
——————IF < 365天|
| IF 365天输出结果
算法,
( L、M需人为给定，x由机器产生)
第i天初，有C(i)+E(i)辆
(其中，C(i)是前夜剩余，E(i)是当天进货)
| | |
| | |
若C(i)+E(i)>x+L
则 兽出y=x
缺货Q=0
不订货E(i+3)=0
当天净收F(i)=10y-
0.8(C(i)+E(i)-y)
交接 C(i+1)=
C(i)+E(i)-y
若xC(i)+E(i) x+L
则 兽出y=x
缺货Q=0
订货E(i+3)=M
当天净收F(i)=10y-
0.8(C(i)+E(i)-y)-75
交接 C(i+1)=
C(i)+E(i)-y
若C(i)+E(i)<x
则兽出y=C(i)+E(i)
缺货Q=x-(C(i)+E(i))
订货E(i+3)=M
当天净收F(i)=10y-
2Q-75
交接C(i+1)=0
模拟计算：对于给定数据L、M，用机器随机产生365个数据（这些数在0至99间均匀分布）分别作为日需求量，计算出一年的总利润。由于日需求量数据是随机产生的，所以即使对固定的L、M，每次计算的结果也会不同，故要计算多次，对全部结果取平均值。下面的Matlab程序对固定的L、M计算了800次，800个结果取均值并输出该均值(程序名：JSJMN1.m),
function yy=JSJMN1(L,M)
for j=1:800
f=zeros(1,365);
x=fix(rand(1,365)*100);
c(1)=115;e(1)=0;e(2)=0;e(3)=0;
for i=1:365
if c(i)+e(i)>x(i)+L
y=x(i);e(i+3)=0;
f(i)=10*y-0.8*(c(i)+e(i)-y);
c(i+1)=c(i)+e(i)-y;
elseif c(i)+e(i)<x(i)
y=c(i)+e(i);q=x(i)-y;e(i+3)=M;
f(i)=10*y-2*q-75;c(i+1)=0;
else
y=x(i);e(i+3)=M;
f(i)=10*y-0.8*(c(i)+e(i)-y)-75;
c(i+1)=c(i)+e(i)-y;
end
end
ff(j)=sum(f);
end
mean(ff)
数据实验：把上述程序输入机器后，执行Matlab命令 L=150;M=180;JSJMN1(L,M) 就的结果9.4万元。 一般地，令L、M取其它一些值，结果见下表：
L\M
50
70
90
110
130
150
170
30
6.3
6.4
6.5
6.5
6.5
6.5
6.5
50
10.5
11.3
11.8
12.1
12.2
12.2
12.1
70
10.9
11.6
12.1
12.40
12.44
12.3
12.0
90
10.7
11.3
11.7
12.0
12.0
11.9
11.6
110
10.5
11.0
11.4
11.5
11.5
11.4
11.1
130
10.2
10.6
10.9
11.0
11.0
150
10.0
10.3
10.5
可见，当L=70,M=130时，有最大值12.44.
答：剩70辆时就立即发订货单，订货量130辆/次，这样可以使全年利润达到最大12.44万元.
二．怎样产生随机数？
用计算机模拟法解题时，通常需要用机器模拟产生服从题给概率分布的随机数。如：顾客对车的需求量是x辆/天，x是从0至99的均匀分布随机数，可用Matlab语句x=fix(rand(1,1)*100)来产生；在实习题第3题中，每个轴承的寿命是随机数，该随即机数必须服从题给的概率分布，怎样产生这个寿命数？
Matlab中，rand命令产生开区间（0，1）上均匀分布随机数，在此基础上，灵活应用rand命令可以产生许多其它性质的随机数据。下面三个例题，介绍了三种性质的随机数据的产生方法。学生要认真学习、深刻体会。
例1：y是一个随机数，它的取值及相应的概率分布如下：
============================================
y的取值,2 6 8 10
概 率,0.13 0.26 0.56 0.05
============================================
请建立函数文件 sjhs.m,产生数据y,
分析说明：用rand命令产生一个数a,则a属于开区间（0，1），现将这个大区间分割为4个小区间 (0,0.13),(0.13,0.39),(0.39,0.95),(0.95,1)，注意到这4个小区间之宽度分别就等于题给的4个概率。 问你，这个a落在第1个小区间的概率是几？落在第3个小区间的概率又是几？ 设计思想：a落在第几个小区间就令y取值为2、6、8、10中的第几个数。
函数文件（程序）：
function y=sjhs
a=rand(1,1);
if a<0.13
y=2;
elseif a<0.39
y=6;
elseif a<0.95
y=8;
else
y=10;
end
例2：随机组队问题
某校派30名学生组成10个队外出参加一个知识竞赛。该竞赛的规则是每3人组成一个队。竞赛涉及的知识分4大类：自然科学、社会科学、文学、艺术。这4类知识在题目中所占比重相同。已知这30人的成绩：
***********************************************************************
学生代号：01 02 03 04 05 06 07 08 ……… 27 28 29 30
自然科学：90 59 76 70 88 82 69 50 ……… 66 85 71 80
社会科学：85 69 76 84 89 81 57 67 ……… 73 69 79 76
文 学：67 88 52 65 84 50 55 58 ……… 80 66 86 86
艺 术：69 90 70 68 86 72 58 78 ……… 71 88 85 77
***********************************************************************
根据经验，假定在一个队内，某类知识3人的成绩较高的两个成绩的平均值代表该队的该类知识水平。一个队的4类知识水平之和为该队的竞赛水平。问怎样组队，可使10个队的竞赛水平之和达到最大？
分析说明：把30人等分成10组，共有几种分赔方案？ 答共有种。若要对全部方案都进行计算检验，则短时间内算不完。现在介绍一种比较粗糙但行之有效的方法来求出近似最佳方案：随机产生一个方案，将30人分成10个队，计算出该方案的竞赛水平；这样的工作重复多次（如：数万次、数十万次），输出竞赛水平最高的方案。此算法的核心是“怎样将30人随机地分成10个队”。
将30人随机地分成10个队：所有人都编号，分别记为1号、2号、3号、……、30号。构造一个2行30列型的矩阵，第1行为 1 2 3 4 …… 30，第2行元素为30个在开区间（0，1）上均匀分布的随机数。再把该矩阵的30个列调序，使得第2行元素为从小到大排列，此时，令第1行的30元素每相邻的3个组成一队即可。
A=[1:30;rand(1,30)];B=sortrows(A',2)';b=B(1,:);
for i=0:9
G(i+1,:)=[b(3*i+1),b(3*i+2),b(3*i+3)];
end
G
例3：电梯环行一周所需时间一个办公大楼，第3层至第12层的上班人数分别是
230，170，320，280，180，180，170，120，90，60.
一台电梯速度为200（单位：m/min）。 只关注上午上班高峰期的需求，电梯从1层满载20人只管送人上楼到4—12层，一旦变空就立刻下行到1层中途不停。已知如下参数：第一层楼高7.62m，第2----11每层楼高3.91m，电梯起动、停止的加速度均为1.22m/s.s,开、关门时间分别为3s，每位乘客进电梯需要1s，出电梯需要0.8s,请用计算机模拟运行，电梯环行一周平均需要多长时间？
分析：核心是要随机产生每个层出电梯的人数，方法如下：总共10层，每层的上班人数在总人数中的比例已知，把开区间（0，1）分割为10个小区间，小区间的宽度分别对应那10个比例数；用命令rand(1,20)产生的20个随机数，有几个落在第i个小区间就表示在第2+i层出电梯的人数是几。
随机产生每个层出电梯的人数：程序如下
clear
rs=[230,170,320,280,180,180,170,120,90,60];bl=rs/sum(rs);
fgd(1)=0;
for j=1:10
fgd(j+1)=fgd(j)+bl(j);
end
sjs=rand(1,20);
cdt=zeros(1,10);
for i=1:20
for j=1:10
if sjs(i)>fgd(j)&sjs(i)<fgd(j+1)
cdt(j)=cdt(j)+1;
end
end
end
[3:12;cdt]
运行该程序，得到这样一个输出：（第一行是楼层，第二行是出电梯人数）
3 4 5 6 7 8 9 10 11 12
1 3 4 2 0 4 1 1 3 1
再运行一次，结果当然会有变化：
3 4 5 6 7 8 9 10 11 12
6 2 4 2 2 0 2 2 0 0
完毕。
一．什么是“计算机模拟”？ 一个事例题：开张营业时，待售自行车115辆；(售价-进货价-营业税)=10元/辆；
保管费：0.8元/天辆；发生缺货(有顾客要购车而无货)时，损失费：2元/天辆；
顾客对车的需求量：x辆/天，x是从0至99的均匀分布随机数；
想要进货时，必须提前3天发订货单，订货费75元/次，如：2007年07月25日发出订货单，则2007年07月28日一上班就到货；
计划：剩L辆时就立即发订货单，订货量M辆/次，请确定L=？M=？ 题毕。
分析构造模拟模型：模拟运行一年365天，令L、M取各种可能的数据，分别计算总利润，使得利润达到最大的数据L、M即为所求。
初始状态
——————————|
|增加一天
(1),这天初：库存C=？，进货E=？
(2),这天当中：售出y=？，发订货单否？
缺货量Q=？
(3),这天末算帐：当天的净收入F=？
——————IF < 365天|
| IF 365天输出结果
算法,
( L、M需人为给定，x由机器产生)
第i天初，有C(i)+E(i)辆
(其中，C(i)是前夜剩余，E(i)是当天进货)
| | |
| | |
若C(i)+E(i)>x+L
则 兽出y=x
缺货Q=0
不订货E(i+3)=0
当天净收F(i)=10y-
0.8(C(i)+E(i)-y)
交接 C(i+1)=
C(i)+E(i)-y
若xC(i)+E(i) x+L
则 兽出y=x
缺货Q=0
订货E(i+3)=M
当天净收F(i)=10y-
0.8(C(i)+E(i)-y)-75
交接 C(i+1)=
C(i)+E(i)-y
若C(i)+E(i)<x
则兽出y=C(i)+E(i)
缺货Q=x-(C(i)+E(i))
订货E(i+3)=M
当天净收F(i)=10y-
2Q-75
交接C(i+1)=0
模拟计算：对于给定数据L、M，用机器随机产生365个数据（这些数在0至99间均匀分布）分别作为日需求量，计算出一年的总利润。由于日需求量数据是随机产生的，所以即使对固定的L、M，每次计算的结果也会不同，故要计算多次，对全部结果取平均值。下面的Matlab程序对固定的L、M计算了800次，800个结果取均值并输出该均值(程序名：JSJMN1.m),
function yy=JSJMN1(L,M)
for j=1:800
f=zeros(1,365);
x=fix(rand(1,365)*100);
c(1)=115;e(1)=0;e(2)=0;e(3)=0;
for i=1:365
if c(i)+e(i)>x(i)+L
y=x(i);e(i+3)=0;
f(i)=10*y-0.8*(c(i)+e(i)-y);
c(i+1)=c(i)+e(i)-y;
elseif c(i)+e(i)<x(i)
y=c(i)+e(i);q=x(i)-y;e(i+3)=M;
f(i)=10*y-2*q-75;c(i+1)=0;
else
y=x(i);e(i+3)=M;
f(i)=10*y-0.8*(c(i)+e(i)-y)-75;
c(i+1)=c(i)+e(i)-y;
end
end
ff(j)=sum(f);
end
mean(ff)
数据实验：把上述程序输入机器后，执行Matlab命令 L=150;M=180;JSJMN1(L,M) 就的结果9.4万元。 一般地，令L、M取其它一些值，结果见下表：
L\M
50
70
90
110
130
150
170
30
6.3
6.4
6.5
6.5
6.5
6.5
6.5
50
10.5
11.3
11.8
12.1
12.2
12.2
12.1
70
10.9
11.6
12.1
12.40
12.44
12.3
12.0
90
10.7
11.3
11.7
12.0
12.0
11.9
11.6
110
10.5
11.0
11.4
11.5
11.5
11.4
11.1
130
10.2
10.6
10.9
11.0
11.0
150
10.0
10.3
10.5
可见，当L=70,M=130时，有最大值12.44.
答：剩70辆时就立即发订货单，订货量130辆/次，这样可以使全年利润达到最大12.44万元.
二．怎样产生随机数？
用计算机模拟法解题时，通常需要用机器模拟产生服从题给概率分布的随机数。如：顾客对车的需求量是x辆/天，x是从0至99的均匀分布随机数，可用Matlab语句x=fix(rand(1,1)*100)来产生；在实习题第3题中，每个轴承的寿命是随机数，该随即机数必须服从题给的概率分布，怎样产生这个寿命数？
Matlab中，rand命令产生开区间（0，1）上均匀分布随机数，在此基础上，灵活应用rand命令可以产生许多其它性质的随机数据。下面三个例题，介绍了三种性质的随机数据的产生方法。学生要认真学习、深刻体会。
例1：y是一个随机数，它的取值及相应的概率分布如下：
============================================
y的取值,2 6 8 10
概 率,0.13 0.26 0.56 0.05
============================================
请建立函数文件 sjhs.m,产生数据y,
分析说明：用rand命令产生一个数a,则a属于开区间（0，1），现将这个大区间分割为4个小区间 (0,0.13),(0.13,0.39),(0.39,0.95),(0.95,1)，注意到这4个小区间之宽度分别就等于题给的4个概率。 问你，这个a落在第1个小区间的概率是几？落在第3个小区间的概率又是几？ 设计思想：a落在第几个小区间就令y取值为2、6、8、10中的第几个数。
函数文件（程序）：
function y=sjhs
a=rand(1,1);
if a<0.13
y=2;
elseif a<0.39
y=6;
elseif a<0.95
y=8;
else
y=10;
end
例2：随机组队问题
某校派30名学生组成10个队外出参加一个知识竞赛。该竞赛的规则是每3人组成一个队。竞赛涉及的知识分4大类：自然科学、社会科学、文学、艺术。这4类知识在题目中所占比重相同。已知这30人的成绩：
***********************************************************************
学生代号：01 02 03 04 05 06 07 08 ……… 27 28 29 30
自然科学：90 59 76 70 88 82 69 50 ……… 66 85 71 80
社会科学：85 69 76 84 89 81 57 67 ……… 73 69 79 76
文 学：67 88 52 65 84 50 55 58 ……… 80 66 86 86
艺 术：69 90 70 68 86 72 58 78 ……… 71 88 85 77
***********************************************************************
根据经验，假定在一个队内，某类知识3人的成绩较高的两个成绩的平均值代表该队的该类知识水平。一个队的4类知识水平之和为该队的竞赛水平。问怎样组队，可使10个队的竞赛水平之和达到最大？
分析说明：把30人等分成10组，共有几种分赔方案？ 答共有种。若要对全部方案都进行计算检验，则短时间内算不完。现在介绍一种比较粗糙但行之有效的方法来求出近似最佳方案：随机产生一个方案，将30人分成10个队，计算出该方案的竞赛水平；这样的工作重复多次（如：数万次、数十万次），输出竞赛水平最高的方案。此算法的核心是“怎样将30人随机地分成10个队”。
将30人随机地分成10个队：所有人都编号，分别记为1号、2号、3号、……、30号。构造一个2行30列型的矩阵，第1行为 1 2 3 4 …… 30，第2行元素为30个在开区间（0，1）上均匀分布的随机数。再把该矩阵的30个列调序，使得第2行元素为从小到大排列，此时，令第1行的30元素每相邻的3个组成一队即可。
A=[1:30;rand(1,30)];B=sortrows(A',2)';b=B(1,:);
for i=0:9
G(i+1,:)=[b(3*i+1),b(3*i+2),b(3*i+3)];
end
G
例3：电梯环行一周所需时间一个办公大楼，第3层至第12层的上班人数分别是
230，170，320，280，180，180，170，120，90，60.
一台电梯速度为200（单位：m/min）。 只关注上午上班高峰期的需求，电梯从1层满载20人只管送人上楼到4—12层，一旦变空就立刻下行到1层中途不停。已知如下参数：第一层楼高7.62m，第2----11每层楼高3.91m，电梯起动、停止的加速度均为1.22m/s.s,开、关门时间分别为3s，每位乘客进电梯需要1s，出电梯需要0.8s,请用计算机模拟运行，电梯环行一周平均需要多长时间？
分析：核心是要随机产生每个层出电梯的人数，方法如下：总共10层，每层的上班人数在总人数中的比例已知，把开区间（0，1）分割为10个小区间，小区间的宽度分别对应那10个比例数；用命令rand(1,20)产生的20个随机数，有几个落在第i个小区间就表示在第2+i层出电梯的人数是几。
随机产生每个层出电梯的人数：程序如下
clear
rs=[230,170,320,280,180,180,170,120,90,60];bl=rs/sum(rs);
fgd(1)=0;
for j=1:10
fgd(j+1)=fgd(j)+bl(j);
end
sjs=rand(1,20);
cdt=zeros(1,10);
for i=1:20
for j=1:10
if sjs(i)>fgd(j)&sjs(i)<fgd(j+1)
cdt(j)=cdt(j)+1;
end
end
end
[3:12;cdt]
运行该程序，得到这样一个输出：（第一行是楼层，第二行是出电梯人数）
3 4 5 6 7 8 9 10 11 12
1 3 4 2 0 4 1 1 3 1
再运行一次，结果当然会有变化：
3 4 5 6 7 8 9 10 11 12
6 2 4 2 2 0 2 2 0 0
完毕。
