计算机模拟_排队论

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

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

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

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

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

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

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

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

编写程序如下：
clear
rand('state',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=zeros(n,1); %a2初始化
a3(1)=a2(1);
if a3(1)<=m
a4(1)=a3(1);a5(1)=0;
else
a4(1)=m;a5(1)=a2(1)-m;
end
for i=2:n
a3(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;
end
end
a=[a1,a2,a3,a4,a5];
sum(a)/n
m =
2
ans =
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个顾客。

程序如下：
tic
rand('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(i)=max(ctime(i),gtime(i-1))+sspan(i); wtime(i)=max(0,gtime(i-1)-ctime(i));
end
result1(j)=sum(wtime)/n;
end
result_1=sum(result1)/m
toc
result_1 =
4.0467
Elapsed time is 0.445770 seconds.
类似地，模拟B型机的程序如下：
tic
rand('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(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.7368
Elapsed time is 1.453880 seconds.
可以用下面的程序与上面的程序比较了解编程的效率问题。

tic
clear
rand('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),gtime(i-1))+exprnd(mu2);
wtime(i)=max(0,gtime(i-1)-ctime(i));
end
result(j)=sum(wtime)/n;
end
result=sum(result)/m
toc
result =
4.2162
Elapsed 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,clear
load 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,clear
load 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.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),'*')
从散点图可以看出，第一阶段基本上是线性关系，第二阶段准备依次用二次、三次、四次曲线来拟合，看哪一个模型的剩余标准差小就选取哪一个模型。

最后求得第一阶
段排沙量y与水流量v之间的预测模型为y = 250.5655v − 373384.4661第二阶段的预测模型为一个四次多项式。

y =−2.7693×10−7v4+0.0018v3−4.092v2+3891.0441v−1.32262749668×106
计算的 Matlab 程序如下：
clc, clear
load 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 e
nihe1_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); rmse1_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表示正常状态，用0表示不正常状态，所得的数据序列如下：1110010011111110011110111111001111111110001101101 111011011010111101110111101111110011011111100111
解设 X (n = 1,L,97) n 为第n 个时段的计算机状态，可以认为它是一个时齐马氏链，状态空间E = {0,1}，编写如下Matlab 程序：
a1='1110010011111110011110111111001111111110001101101';
a2='111011011010111101110111101111110011011111100111';
a=[a1 a2];
f00=length(findstr('00',a))
f01=length(findstr('01',a))
f10=length(findstr('10',a))
f11=length(findstr('11',a))
f00 =
8
f01 =
18
f10 =
18
f11 =
52
或者把上述数据序列保存到纯文本文件data1.txt中，存放在Matlab下的work 子目录中，编写程序如下：
%clc,clear
format rat
%fid=fopen('data1.txt','r');
%a=[];
%while (~feof(fid))
%a=[a fgetl(fid)];
%end
for i=0:1
for j=0:1
s=[int2str(i),int2str(j)];
f(i+1,j+1)=length(findstr(s,a));
end
end
fs=sum(f');
for i=1:2
f(i,:)=f(i,:)/fs(i);
end
f
f =
4/13 9/13
9/35 26/35。