计算机模拟法相关知识 什么是“计算机模拟”?一个事例 题:开张营业时,待售自行车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天 输出结果 算法: (L、M需人为给定,x由机器产生) 第i天初,有C(i)+E(i)辆 (其中,C(i)是前夜剩余,E(i)是当天进货) 若C(i)+E(i)>x+L x≤C(i)+E(i)≤x+L若C(i)+(i)<x 则兽出y=x 兽出y=x 则兽出y=C(i)+E(i) 缺货Q=0 缺货Q=C 缺货Q=x-(C(i)+E(i) 不订货E(i+3)2=0订货E(i+3)=M 订货E(i+3)=M 当天净收F(i)=10y-当天净收F(i)=10y- 当天净收F(i)=10y- 0.8(c(i)+E(i)-y)0.8c(i)+E(i)-y)-75 2Q-75 交接C(i+1)= 交接C(i+1)= 交接C(i+1)=0 C(i)+E(i)-y C(i)+E(i)-y
计算机模拟法相关知识 一．什么是“计算机模拟”？ 一个事例 题：开张营业时，待售自行车 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 若 x C(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=JS JMNI(L, M) forj=1:800 f=zeros(1, 365) x=fix(rand(1,365)*100) c(1)=115;e(1)=0;e(2)=0;e(3)=0 fori=1:365 if c(i)+e(i)>x(i)+L f(i)=10*y-0.8*(c(i)+e(i)-y) c(i+1)=c(i)+e(i)-y; leif c(i)+e (i)<x (i) =c(i)+e(i);q=x(i)-y;e(i+3)=M; f(i)=10*y-2*q-75;c(i+1)=0; y-X f(i)=10*y-0.8*(c(i)+e(i)-y)-75; c(i+1)=c(i)+e(i)-y; end ff(j)=sum(f) end mean(ff 数据实验:把上述程序输入机器后,执行 Matlab命令L=150;M=180; JSJMN(L,M) 就的结果9.4万元。一般地,令L、M取其它一些值,结果见下表: LM507090110130150170 306.36.46.56.56.56.56.5 5010.511.311.812.112.212.212.1 7010.911.6|12112401244123|120 9010.711.31.712.012.011.916 110[10.5|11.011411.511.5111.4111 13010.210.610.911.011.0 15010.010.310.5 可见,当L=70,M=130时,有最大值12.44. 答:剩70辆时就立即发订货单,订货量130辆/次,这样可以使全年利润达到 最大12.44万元
模拟计算：对于给定数据 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 语句κ=fix(rand(1,1)*100)来产生;在实习题第3题中,每个轴承的寿命是随机数, 该随即机数必须服从题给的概率分布,怎样产生这个寿命数? Matlab中,rand命令产生开区间(0,1)上均匀分布随机数,在此基础上 灵活应用rand命令可以产生许多其它性质的随机数据。下面三个例题,介绍了 种性质的随机数据的产生方法。学生要认真学习、深刻体会。 例1:y是一个随机数,它的取值及相应的概率分布如下 y的取值 10 概率:0.130.260.5 请建立函数文件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=s jhs a=rand(1, 1) ifa<0.13 elseif a<0. 39 elseif a<0. 95 end 例2:随机组队问题 某校派30名学生组成10个队外出参加一个知识竞赛。该竞赛的规则是每3人 组成一个队。竞赛涉及的知识分4大类:自然科学、社会科学、文学、艺术。这4 类知识在题目中所占比重相同。已知这30人的成绩: *半*水水半*水水水水水*水水冰半冰*水冰冰*水冰冰*水**水水水**水冰**水*水冰**水水冰水**水水水半*水水冰*水*半*冰冰 学生代号:0102030405060708………27282930 自然科学:9059767088826950………66857180 社会科学:8569768489815767……73697976
二．怎样产生随机数？ 用计算机模拟法解题时，通常需要用机器模拟产生服从题给概率分布的随机 数。如：顾客对车的需求量是 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
文学:6788526584505558… 80668686 艺术:6990706886725878… 71888577 冰*水冰冰水水冰水*水*冰水冰*水冰冰*水水冰*水*水*水冰水*冰*水水冰水水冰冰水水水*水半**水水*水*冰*水水水水水冰水*冰 根据经验,假定在一个队内,某类知识3人的成绩较高的两个成绩的平均值代表该 队的该类知识水平。一个队的4类知识水平之和为该队的竞赛水平。问怎样组队, 可使10个队的竞赛水平之和达到最大? 分析说明:把30人等分成10组,共有几种分赔方案?答共有30种。 (3)°·10 若要对全部方案都进行计算检验,则短时间内算不完。现在介绍一种比较粗糙但行 之有效的方法来求出近似最佳方案:随机产生一个方案,将30人分成10个队,计 算出该方案的竞赛水平;这样的工作重复多次(如:数万次、数十万次),输出竞 赛水平最高的方案。此算法的核心是“怎样将30人随机地分成10个队” 将30人随机地分成10个队:所有人都编号,分别记为1号、2号、3号、 30号。构造一个2行30列型的矩阵,第1行为1234……30,第2行元 素为30个在开区间(0,1)上均匀分布的随机数。再把该矩阵的30个列调序,使 得第2行元素为从小到大排列,此时,令第1行的30元素每相邻的3个组成一队 即可。 A=[1:30;rand(1,30)];B= sorrows(A',2)’;b=B(1,:) for i=0: 9 G(i+1,)=[b(3*i+1),b(3*1+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.62皿,第2—-11每层楼髙3.9l皿,电梯起动、停止 的加速度均为1.22m/s.s,开、关门时间分别为3s,每位乘客进电梯需要ls,出 电梯需要0.8s.请用计算机模拟运行,电梯环行一周平均需要多长时间? 分析:核心是要随机产生每个层出电梯的人数,方法如下:总共10层,每层 的上班人数在总人数中的比例已知,把开区间(0,1)分割为10个小区间,小区 间的宽度分别对应那10个比例数:用命令rand(1,20)产生的20个随机数,有几个 落在第i个小区间就表示在第2+i层出电梯的人数是几。 随机产生每个层出电梯的人数:程序如 rs=[230,170,320,280,180,180,170,120,90,60];bl=s/sum(rs) fgd(1)=0; for j=1: 10
文 学：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 组，共有几种分赔方案？ 答共有 (3!) 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) s js=rand (1, 20 cdt=zeros(1, 10) for i=1: 20 for j=1: 10 if sjs(i)>fgd (j)&sjs(i)<fgd(+1) (=cdt ()+1 end end en [3:12;cdt] 运行该程序,得到这样一个输出:(第一行是楼层,第二行是出电梯人数) 1011 4342 42 再运行一次,结果当然会有变化 7072 8480 91011 12 6 4 0 完毕
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 完毕