当排队系统的到达间隔时间和服务时间的概率分布很复杂时，或不能用公式给出时，那么就不能用解析法求解。

这就需用随机模拟法求解，现举例说明。

例 1 设某仓库前有一卸货场，货车一般是夜间到达，白天卸货，每天只能卸货2车，若一天内到达数超过2车，那么就推迟到次日卸货。

根据表3所示的数据，货车到达数的概率分布（相对频率）平均为1.5车/天，求每天推迟卸货的平均车数。

解服从指数分布（这是定长服务时间）。

随机模拟法首先要求事件能按历史的概率分布规律出现。

模拟时产生的随机数与事件的对应关系如表4。

表 4 到达车数的概率及其对应的随机数我们用 a1 表示产生的随机数，a2 表示到达的车数，a3 表示需要卸货车数，a4表示实际卸货车数，a5 表示推迟卸货车数。

编写程序如下：clearrand('state',sum(100*clock));n=50000;m=2a1=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=zeros(n,1); %a2初始化a3(1)=a2(1);if a3(1)<=ma4(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)<=ma4(i)=a3(i);a5(i)=0;elsea4(i)=m;a5(i)=a3(i)-m;endenda=[a1,a2,a3,a4,a5];sum(a)/nm =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型机采取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=ck−1+ik，gk=max(ck,gk−1)+ sk，wk=max(0,gk−1− 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:mcspan=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:nctime(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));endresult1(j)=sum(wtime)/n;endresult_1=sum(result1)/mtocresult_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:mcspan=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:nctime(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)];endresult2(j)=sum(wtime)/n;endresult_2=sum(result2)/mtocresult_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:mctime(1)=exprnd(mu1);gtime(1)=ctime(1)+exprnd(mu2);wtime(1)=0;for i=2:nctime(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));endresult(j)=sum(wtime)/n;endresult=sum(result)/mtocresult =4.2162Elapsed time is 3.854620 seconds.黄河小浪底调水调沙问题5.1 问题的提出2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功。

整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正常供水结束。

小浪底水利工程按设计拦沙量为75.5 亿m3，在这之前，小浪底共积泥沙达14.15亿t。

这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700m3/s，使小浪底水库的排沙量也不断地增加。

表7是由小浪底观测站从6月29日到7月10日检测到的试验数据。

表 7 观测数据（1）给出估计任意时刻的排沙量及总排沙量的方法；（2）确定排沙量与水流量的关系。

5.2 模型的建立与求解已知给定的观测时刻是等间距的，以6月29日零时刻开始计时，则各次观测时刻（离开始时刻6月29日零时刻的时间）分别为ti =3600(12i−4)，i=1,2,L,24，其中计时单位为秒。

第1次观测的时刻t1=28800，最后一次观测的时刻t24=1022400 。

记第i（i= 1,2,L,24）次观测时水流量为vi，含沙量为ci，则第i次观测时的排沙量为yi=ci*vi 。

有关的数据见表8。

对于问题（1），根据所给问题的试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。

考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，我们采用三次样条函数进行插值。

利用 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=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,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(pp,'pp') %把 B 样条函数转化为 pp 格式TL=quadl(@(tt)fnval(pp,tt),t1,t2)对于问题（2），研究排沙量与水量的关系，从试验数据可以看出，开始排沙量是随着水流量的增加而增长，而后是随着水流量的减少而减少。

显然，变化规律并非是线性的关系，为此，把问题分为两部分，从开始水流量增加到最大值 2720m3/s（即增长的过程）为第一阶段，从水流量的最大值到结束为第二阶段，分别来研究水流量与排沙量的关系。

画出排沙量与水流量的散点图。

画散点图的程序如下：load data.txtliu=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),'*')从散点图可以看出，第一阶段基本上是线性关系，第二阶段准备依次用二次、三次、四次曲线来拟合，看哪一个模型的剩余标准差小就选取哪一个模型。