计算机模拟_排队论

1、当排队系统的到达间隔时间和服务时间的概率分布很复杂时，或不能用公式给出 时，那么就不能用解析法求解。这就需用随机模拟法求解，现举例说明。例 1 设某仓库前有一卸货场，货车一般是夜间到达，白天卸货，每天只能卸货2车，若一天内到达数超过2车，那么就推迟到次日卸货。根据表3所示的数据，货车到达数的概率分布（相对频率）平均为1.5车/天，求每天推迟卸货的平均车数。表 3 到达车数的概率到达车数012345 6概率0.230.300.300.10.050.020.00解这是单服务台的排队系统，可验证到达车数不服从泊松分布，服务时间也不服从指数分布（这是定长服务时间）。随机模拟法首先要求事件能按历史的概率

2、分布规律出现。模拟时产生的随机数与事件的对应关系如表4。表 4 到达车数的概率及其对应的随机数到达车数概率累积概率对应的随机数0123450.230.300.300.10.050.020.230.530.830.930.981.000 x 0.230.23 x 0.530.53 x 0.830.83 x 0.930.93 x 0.980.98 x 1.00我们用 a1 表示产生的随机数，a2 表示到达的车数，a3 表示需要卸货车数，a4 表 示实际卸货车数，a5 表示推迟卸货车数。编写程序如下：clear rand(state,sum(100*clock);n=50000;m=2 a1=ran

3、d(n,1);a2=a1; %a2初始化a2(find(a10.23)=0; a2(find(0.23=a1&a10.53)=1; a2(find(0.53=a1&a10.83)=2; a2(find(0.83=a1&a10.93),1)=3; a2(find(0.93=a1&a1=0.98)=5;a3=zeros(n,1);a4=zeros(n,1);a5=zeros(n,1); %a2初始化a3(1)=a2(1);if a3(1)=m a4(1)=a3(1);a5(1)=0;elsea4(1)=m;a5(1)=a2(1)-m;endfor i=2:na3(i)=a2(i)+a5(i-1);

4、 if a3(i)=m a4(i)=a3(i);a5(i)=0; else a4(i)=m;a5(i)=a3(i)-m;endend a=a1,a2,a3,a4,a5;sum(a)/n m = 2ans = 0.4985 1.4909 2.3782 1.4909 0.8874 例2银行计划安置自动取款机，已知A型机的价格是B型机的2倍，而A型机的性能平均服务率也是B型机的2倍，问应该购置1台 A 型机还是2台 B 型机。为了通过模拟回答这类问题，作如下具体假设，顾客平均每分钟到达1位， A 型 机的平均服务时间为0.9分钟，B 型机为1.8分钟，顾客到达间隔和服务时间都服从指数分布，2台B型机

5、采取M/M/2模型（排一队），用前100名顾客（第 1 位顾客到达时取款机前为空）的平均等待时间为指标，对A型机和B型机分别作1000次模拟，进行比较。理论上已经得到，A型机和B型机前100名顾客的平均等待时间分别为1(100)=4.13，2(100)=3.70，即 B 型机优。对于M/M/1模型，记第k位顾客的到达时刻为ck，离开时刻为gk，等待时间为wk，它们很容易根据已有的到达间隔ik和服务时间sk按照以下的递推关系得到（w1 = 0，设c1,g1已知）：ck=ck1+ik，gk=max(ck,gk1)+ sk，wk=max(0, gk1 ck)， k=2,3,L。在模拟A型机时，用cs

6、pan表示到达间隔时间，sspan表示服务时间，ctime表示到达时间，gtime表示离开时间，wtime表示等待时间。我们总共模拟了m次，每次n个顾客。程序如下：ticrand(state,sum(100*clock);n=100;m=1000;mu1=1;mu2=0.9; for j=1:m cspan=exprnd(mu1,1,n); sspan=exprnd(mu2,1,n); ctime(1)=cspan(1); gtime(1)=ctime(1)+sspan(1); wtime(1)=0; for i=2:n ctime(i)=ctime(i-1)+cspan(i); gtime(

7、i)=max(ctime(i),gtime(i-1)+sspan(i); wtime(i)=max(0,gtime(i-1)-ctime(i); end result1(j)=sum(wtime)/n;endresult_1=sum(result1)/mtoc result_1 = 4.0467Elapsed time is 0.445770 seconds. 类似地，模拟B型机的程序如下：ticrand(state,sum(100*clock);n=100;m=1000;mu1=1;mu2=1.8;for j=1:m cspan=exprnd(mu1,1,n); sspan=exprnd(m

8、u2,1,n); ctime(1)=cspan(1); ctime(2)=ctime(1)+cspan(2); gtime(1:2)=ctime(1:2)+sspan(1:2); wtime(1:2)=0;flag=gtime(1:2); for i=3:n ctime(i)=ctime(i-1)+cspan(i); gtime(i)=max(ctime(i),min(flag)+sspan(i); wtime(i)=max(0,min(flag)-ctime(i); flag=max(flag),gtime(i); end result2(j)=sum(wtime)/n;end result

9、_2=sum(result2)/m toc result_2 = 3.7368Elapsed time is 1.453880 seconds. 可以用下面的程序与上面的程序比较了解编程的效率问题。ticclearrand(state,sum(100*clock);n=100;m=1000;mu1=1;mu2=0.9;for j=1:m ctime(1)=exprnd(mu1); gtime(1)=ctime(1)+exprnd(mu2); wtime(1)=0; for i=2:n ctime(i)=ctime(i-1)+exprnd(mu1); gtime(i)=max(ctime(i),

10、gtime(i-1)+exprnd(mu2); wtime(i)=max(0,gtime(i-1)-ctime(i); end result(j)=sum(wtime)/n;end result=sum(result)/mtoc result = 4.2162Elapsed time is 3.854620 seconds. 黄河小浪底调水调沙问题5.1问题的提出2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功。整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正

11、常供水结束。小浪底水利工程按设计拦沙量为75.5 亿m3，在这之前，小浪底共积泥沙达14.15亿t。这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700m3/s，使小浪底水库的排沙量也不断地增加。表7是由小浪底观测站从6月29日到7月10日检测到的试验数据。表 7 观测数据日期6.296.307.17.27.37.4时间8:0020:008:0020:008:0020:008:0020:008

12、:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10时间8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085现在，根据试验数据建立数学模型研究下面的问题：（1）给出估计任意时刻的排沙量及总排沙量

13、的方法；（2）确定排沙量与水流量的关系。5.2模型的建立与求解已知给定的观测时刻是等间距的，以6月29日零时刻开始计时，则各次观测时刻（离开始时刻6月29日零时刻的时间）分别为ti =3600(12i4)，i=1,2,L,24，其中计时单位为秒 。第1次观测的时刻t1=28800，最后一次观测的时刻t24=1022400 。记第i（i= 1,2,L,24）次观测时水流量为vi，含沙量为ci，则第i次观测时的排沙量为yi=ci*vi 。有关的数据见表8。表 8 插值数据对应关系 单位：排沙量为 kg节点12345678时刻2880072000115200158400201600244800288

14、000331200排沙量57600114000157500187000207000235200250000265200节点910111213141516时刻374400417600460800504000547200590400633600676800排沙量286200302400312800307400306800300000271400231000节点1718192021222324时刻7200007632008064008496008928009360009792001022400排沙量1600001110009100054000455003000080004500对于问题（1），根据所

15、给问题的试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，我们采用三次样条函数进行插值。利用 MATLAB 函数，求出三次样条函数，得到排沙量y=y(t)与时间的关系，然后进行积分，就可以得到总的排沙量最后求得总的排沙量为1.844 109 t，计算的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始数据格式把水流量和排沙量排成4行12列liu=data(1,3,:);liu=liu;liu=liu(:);sha=data(2,4,:);sha=s

16、ha;sha=sha(:);y=sha.*liu;y=y;i=1:24;t=(12*i-4)*3600;t1=t(1);t2=t(end);pp=csape(t,y);xsh=pp.coefs%求得插值多项式的系数矩阵，每一行是一个区间上多项式的系数。TL=quadl(tt)ppval(pp,tt),t1,t2)也可以利用 3 次 B 样条函数进行插值，求得总的排沙量也为1.844 109 t，计算 的 Matlab 程序如下：clc,clearload data.txt%data.txt 按照原始数据格式把水流量和排沙量排成4行12列liu=data(1,3,:);liu=liu;liu=l

17、iu(:);sha=data(2,4,:);sha=sha;sha=sha(:);y=sha.*liu;y=y;i=1:24;t=(12*i-4)*3600;t1=t(1);t2=t(end);pp=spapi(4,t,y)%三次 B 样条pp2=fn2fm(pp,pp)%把 B 样条函数转化为 pp 格式TL=quadl(tt)fnval(pp,tt),t1,t2)对于问题（2），研究排沙量与水量的关系，从试验数据可以看出，开始排沙量是随着水流量的增加而增长，而后是随着水流量的减少而减少。显然，变化规律并非是线性 的关系，为此，把问题分为两部分，从开始水流量增加到最大值 2720m3/s（即

18、增长的过 程）为第一阶段，从水流量的最大值到结束为第二阶段，分别来研究水流量与排沙量的关系。画出排沙量与水流量的散点图。画散点图的程序如下：load data.txt liu=data(1,3,:); liu=liu;liu=liu(:); sha=data(2,4,:); sha=sha;sha=sha(:); y=sha.*liu;subplot(1,2,1), plot(liu(1:11),y(1:11),*)subplot(1,2,2), plot(liu(12:24),y(12:24),*)从散点图可以看出，第一阶段基本上是线性关系，第二阶段准备依次用二次、三次、 四次曲线来拟合，看

19、哪一个模型的剩余标准差小就选取哪一个模型。最后求得第一阶段排沙量y与水流量v之间的预测模型为y = 250.5655v 373384.4661第二阶段的预测模型为一个四次多项式。y =2.7693107v4+0.0018v34.092v2+3891.0441v1.32262749668106计算的 Matlab 程序如下：clc, clearload data.txt%data.txt 按照原始数据格式把水流量和排沙量排成4行12列liu=data(1,3,:); liu=liu; liu=liu(:); sha=data(2,4,:); sha=sha; sha=sha(:); y=sha.

20、*liu;%以下是第一阶段的拟合format long enihe1_1=polyfit(liu(1:11),y(1:11),1)%拟合一次多项式，系数排列从高次幂到低次幂nihe1_2=polyfit(liu(1:11),y(1:11),2)yhat1_1=polyval(nihe1_1,liu(1:11);%求预测值yhat1_2=polyval(nihe1_2,liu(1:11);%以下求误差平凡和与剩余标准差cha1_1=sum(y(1:11)-yhat1_1).2); rmse1_1=sqrt(cha1_1/9)cha1_2=sum(y(1:11)-yhat1_2).2); rmse

21、1_2=sqrt(cha1_2/8)%以下是第二阶段的拟合for j=1:3 str1=char(nihe2_ int2str(j) =polyfit(liu(12:24),y(12:24), int2str(j+1) ); eval(str1) str2=char(yhat2_ int2str(j) =polyval(nihe2_ int2str(j) ,liu(12:24);); eval(str2) str3=char(cha2_ int2str(j) =sum(y(12:24)-yhat2_ int2str(j) ).2); rmse2_ int2str(j) =sqrt(cha2_ int2str(j) /(11-j); eval(str3)end 例:当实际问题可以用马尔可夫链来描述时，首先要确定它的状态空间及参数集合，然后确定它的一步转移概率。关于这一概率的确定，可以由问题的内在规律得到，也可以由过去经验给出，还可以根据观测数据来估计。例4 某计算机机房的一台计算机经常出故障，研究者每隔15分钟观察一次计算机的运行状态，收集了24小时的数据（共作97次观察）。用1表示正常状态，用