计算机模拟-排队论

1、word可编辑当排队系统的到达间隔时间和效劳时间的概率分布很复杂时，或不能用公式给出 时，那么就不能用解析法求解。这就需用随机模拟法求解，现举例说明。例 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

3、',sum(100*clock);n=50000;m=2 a1=rand(n,1);a2=a1; %a2初始化a2(find(a1<0.23)=0; a2(find(0.23<=a1&a1<0.53)=1; a2(find(0.53<=a1&a1<0.83)=2; a2(find(0.83<=a1&a1<0.93),1)=3; a2(find(0.93<=a1&a1<0.98),1)=4; a2(find(a1>=0.98)=5;a3=zeros(n,1);a4=zeros(n,1);a5=ze

4、ros(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); 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型机的性能平

5、均效劳率也是B型机的2倍，问应该购置1台 A 型机还是2台 B 型机。为了通过模拟答复这类问题，作如下具体假设，顾客平均每分钟到达1位， A 型 机的平均效劳时间为0.9分钟，B 型机为1.8分钟，顾客到达间隔和效劳时间都服从指数分布，2台B型机采取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，

6、它们很容易根据已有的到达间隔ik和效劳时间sk按照以下的递推关系得到w1 = 0，设c1,g1：ck=ck1+ik，gk=max(ck,gk1)+ sk，wk=max(0, gk1 ck)， k=2,3,L。在模拟A型机时，用cspan表示到达间隔时间，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); ssp

7、an=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(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 se

8、conds. 类似地，模拟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(mu2,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

9、(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_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

10、=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),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

11、seconds. 黄河小浪底调水调沙问题5.1问题的提出2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功。整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正常供水结束。小浪底水利工程按设计拦沙量为75.5 亿m3，在这之前，小浪底共积泥沙达14.15亿t。这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29

12、日先后到达小浪底，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:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10时间8:0020:008:0020:008:0020:008:00

13、20:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085现在，根据试验数据建立数学模型研究下面的问题：1给出估计任意时刻的排沙量及总排沙量的方法；2确定排沙量与水流量的关系。5.2模型的建立与求解给定的观测时刻是等间距的，以6月29日零时刻开始计时，那么各次观测时刻离开始时刻6月29日零时刻的时间分别为ti =3600(12i4)，i=1,2,L,24，其中计时单位为秒 。第1次观测的时刻t1=28800，最后一次观测的时刻t24=1022400

14、。记第ii= 1,2,L,24次观测时水流量为vi，含沙量为ci，那么第i次观测时的排沙量为yi=ci*vi 。有关的数据见表8。表 8 插值数据对应关系 单位：排沙量为 kg节点12345678时刻2880072000115200158400202200244800288000331200排沙量57600114000157500187000207000235200250000265200节点910111213141516时刻374400417600460800504000547200590400633600676800排沙量2862003024003128003074003068003000

15、00271400231000节点1718192021222324时刻7200007632008064008496008928009360009792001022400排沙量1600001110009100054000455003000080004500对于问题1，根据所给问题的试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，我们采用三次样条函数进行插值。利用 MATLAB 函数，求出三次样条函数，得到排沙量y=y(t)与时间的关系，然后进行积分，就可以得到总的排沙量最后求得总的排沙量为1.84

16、4 ×109 t，计算的 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.*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,

17、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=liu(:);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

18、(pp,'pp')%把 B 样条函数转化为 pp 格式TL=quadl(tt)fnval(pp,tt),t1,t2)对于问题2，研究排沙量与水量的关系，从试验数据可以看出，开始排沙量是随着水流量的增加而增长，而后是随着水流量的减少而减少。显然，变化规律并非是线性 的关系，为此，把问题分为两局部，从开始水流量增加到最大值 2720m3/s即增长的过 程为第一阶段，从水流量的最大值到结束为第二阶段，分别来研究水流量与排沙量的关系。画出排沙量与水流量的散点图。画散点图的程序如下：load data.txt liu=data(1,3,:); liu=liu'liu=liu(:)

19、; 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),'*')从散点图可以看出，第一阶段根本上是线性关系，第二阶段准备依次用二次、三次、 四次曲线来拟合，看哪一个模型的剩余标准差小就选取哪一个模型。最后求得第一阶段排沙量y与水流量v之间的预测模型为y = 250.5655v 373384.4661第二阶段的预测模型为一个四次多项式。y =2.7

20、693×107v4+0.0018v34.092v2+3891.0441v1.32262749668×106计算的 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.*liu;%以下是第一阶段的拟合format long enihe1_1=polyfit(liu(1:11),y(1:11),1)%拟合一次多项式，系

21、数排列从高次幂到低次幂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); rmse1_2=sqrt(cha1_2/8)%以下是第二阶段的拟合for j=1:3 str1=char('nihe2_' int2str(j

22、) '=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 例:当实际问题可以用马尔