系统模拟实验的三个案例

系统模拟实验的三个案例
实验案例 赶上火车的概率 1实验案例
1.1 赶上火车的概率 1.1.1 问题描述
如图，一列火车从A 站开往B 站，某人每天赶往B 站上这趟火车。

他已了解到：
（1） 火车从A 站到B 站的运行时间是均值为30分钟，标准差为2分钟的随机变量； （2） 火车在下午大约1点离开A 站，离开时刻的频率分布如下：
问他能赶上火车的概率是多少？ 1.1.2 变量说明
1T ：火车从A 站出发的时刻；
2T ：火车从A 站到B 站的运行时间；单位：分钟 3T ：他到达B 站的时刻
1.1.3 问题分析与假设
此问题包含多个随机因素。

这里假设1T ，2T ，3T 都是随机变量，其中2T 服从正态分布。

1.1.4 模型建立
很显然，他能及时赶上火车的条件是：213T T T +<。

为了简化计算，将下午1点记为
初始时刻。

1T 和
3T 的分布律如下：
为了模拟随机变量。

如果r 为在)1,0(均匀分布的随机数，为了模拟随机变量31,T T ，可
以通过如下方法。

⎪⎩⎪
⎨⎧≤≤<≤<≤=19.0,109
.07.0,57.00,01r r r t ，
⎪⎪⎩⎪⎪⎨⎧≤≤<≤<≤<≤=0.19.0,349
.07.0,327.03.0303.00,283r r r r t ，。

其中，1t 和
3t 分别用来模拟随机变量1T 和3T 。

1.1.5 模拟算法
变量说明：
k 临时变量，存储当前累计模拟次数 count 存储赶上火车的次数
第 1 步 输入模拟次数n 第 2 步 k=1，count=0
第 3 步 当k<=n ，执行第4步，否则执行第12步 第 4 步 生成均匀分布随机数赋给r
第 5 步 由r 及公式确定T1模拟火车出发时刻 第 6 步 生成均匀分布随机数赋给r ；
第 7 步 由r 及公式确定T3模拟人达到时刻
第 8 步 生成正态分布随机数T2模拟火车运行时间 第 9 步 IF T1+T2 > T3, count=count+1，END 第 10 步 k = k + 1 第 11 步 执行第3步
第 12 步 输出赶上火车频率p=count/n
1.1.6 模拟程序 %sim_train.m
total=input('输入模拟次数：'); count=0;
for i=1:total,
%模拟随机变量t1（火车从A 站出发的时刻） rt1=rand; if rt1<0.7 T1=0;
elseif rt1>=0.7 & rt1<0.9 T1=5; else T1=10; end
%模拟随机变量t2（火车的运行时间） T2=30+randn*2;
%模拟随机变量t3（他到达B站的时刻）
rt3=rand;
if rt3<0.3
T3=28;
elseif rt3>=0.3 & rt3<0.7
T3=30;
elseif rt3>=0.7 & rt3<0.9
T3=32;
else
T3=34;
end
if T3 < T1 + T2,%赶上了
count=count+1;
end
end%for
prob=count/total
1.1.7 模拟结果
命令行中输入以下语句：
sim_train
运行结果输出：
输入模拟次数：100
prob =
0.6302
此次运行结果显示赶上火车的近似概率为0.6左右。

下面列表给出多次运行模拟程序的结果。

1.1.8 评价与改进方向
为了计算赶上火车的概率，本文采用了随机系统模拟的方法。

如果能够从模型出发，对赶上火车的概率进行近似计算，然后与模拟结果进行对比，这样模拟会更有说明力。

1.1.9 思考题
（1）请思考用其它方法计算赶上火车的概率或近似概率。

（2）如果要使得他赶上火车的概率大于95%，你有什么办法？结合上面的数学模型及模拟程序来思考。

（3）通过该问题的建模求解，你能归纳出一般系统模拟的方法步骤么？
实验案例理发店模拟
1 实验案例 1.1 案例：理发店模拟
例子：一个理发店有两位服务员A和B，顾客们随机到达店内，其中60％的顾客仅需剪发，每位花5分钟时间，另外40％顾客既要剪发又要洗发，每位用时8分钟。

理发店是个含有多种随机因素的系统，请对该系统进行模拟，并对其进行评判。

（准备怎么做）
可供参考内容
“排队论”，“系统模拟”，“离散系统模拟”，“事件调度法”
1.1.1 问题分析
理发店系统包含诸多随机因素，为了对其进行评判就是要研究其运行效率，从理发店自身利益来说，要看服务员工作负荷是否合理，是否需要增加员工等考虑。

从顾客角度讲，还要看顾客的等待时间，顾客的等待队长，如等待时间过长或者等待的人过多，则顾客会离开。

理发店系统是一个典型的排队系统，可以用排队论有关知识来研究。

1.1.2 模型假设
1． 60％的顾客只需剪发，40％的顾客既要剪发，又要洗发；
2． 每个服务员剪发需要的时间均为5分钟，既剪发又洗发则花8分钟； 3． 顾客的到达间隔时间服从指数分布； 4． 服务中服务员不休息。

1.1.3 变量说明
u ：剪发时间（单位：分钟），u=5m ；
v: 既剪发又理发花的时间（单位：分钟），v=8m ；
T ： 顾客到达的间隔时间，是随机变量，服从参数为 λ的指数分布，（单位：分钟）
T 0：顾客到达的平均间隔时间（单位：秒），T 0＝ λ1
；
1.1.4 模型建立
由于该系统包含诸多随机因素，很难给出解析的结果，因此可以借助计算机模拟对该系统进行模拟。

考虑一般理发店的工作模式，一般是上午9：00开始营业，晚上10：00左右结束，且一般是连续工作的，因此一般营业时间为13小时左右。

这里以每天运行12小时为例，进行模拟。

这里假定顾客到达的平均间隔时间T 0服从均值3分钟的指数分布， 则有
3小时到达人数约为 60360
3=⨯人， 6小时到达人数约为 120
360
6=⨯人， 10小时到达人数约为 200360
10=⨯人，
这里模拟顾客到达数为60人的情况。

（如何选择模拟的总人数或模拟总时间） 1.1.5 系统模拟
根据系统模拟的一般方法，需要考虑系统的如下数据、参数。

1.
2. 实体：两名服务员、顾客们
3. 事件：
4. 活动：
在系统模拟时，为了研究系统的整体情况，这里考虑顾客到达后不离开，且等待队长不限。

要考虑如果服务员均空闲时，顾客先选择谁服务？
要考虑模拟的时间设置还有顾客数目。

模拟终止条件是根据顾客数目还是根据营业时间终止？
1.1.6 计算机模拟算法设计
有无参考算法？
离散系统仿真算法：事件调度法
1.1.7 计算机模拟程序
顾客到达的间隔时间T的计算机产生方法，利用
T=
rand
T
rand ln
ln
1
-
=
-
λ，
%理发店系统的模拟（案例分析之一）
%关键词：面向事件的计算机模拟技术
clear all
curclock=0;%当前时刻，动态变化
totalcustomer=0;%总共服务的顾客数
numsrv=2;
srvstatus=zeros(numsrv,5);%服务员有关数据
%srvstatus 第1列：服务状态（0空闲，1正在服务）；第2列：当前服务顾客编号；
% 第3列：当前服务结束时刻；第4列：服务员空闲时间；第5列：服务的顾客总数
endtime =0;%结束时间
waiting=[];%等待队列数据
%waiting 第1列：顾客编号；第2列：顾客到达时刻；第3列：顾客开始接受服务时刻；% 第4列：接受服务时间；第5列：顾客结束服务时刻；第6列：间隔时间
cur=zeros(1,6);%当前产生顾客的数据，对应关系同waiting
avgwaitlen=[];%平均等待队长
avgwaittime =[];% 平均等待时间
ujiange=5;%平均间隔时间
finished=0;
numsimucustumer=yesinput('输入等待模拟的顾客数：',10,[10 1000]);
while finished==0,
if totalcustomer < numsimucustumer
%产生一个顾客的到达及其有关性质的数据
totalcustomer = totalcustomer+1;
jiange= -log(rand)*ujiange;%与上一个顾客的到达的间隔时间
curclock = curclock + jiange;
cur(1)= totalcustomer ;% 第1列：顾客编号
cur(2) = curclock;%第2列：顾客到达时刻
cur(6) = jiange;第6列：间隔时间
%下面产生接受服务时间（可改进模型）
if rand<0.6, %产生顾客有关性质：这里是产生接受服务时间
cur(4) = 5;
else
cur(4) = 8;
end
%放入等待队列
if isempty(waiting),
waiting= cur;
else
[m,n]=size(waiting);
waiting(m+1,:)= cur;
end
else
curclock = curclock + (-log(rand)*ujiange);
end%if totalcustomer<
%分配等待队列（看是否有服务员空闲，如果有则分配；否则继续执行）
%处理服务员的服务状态
for i=1:numsrv,
if srvstatus(i,1)==1 & srvstatus(i,3) <= curclock,
srvstatus(i,1)=0;%设置为空闲状态
srvstatus(i,4)= curclock-srvstatus(i,3);%目前已经空闲的时间elseif srvstatus(i,1)==1 & srvstatus(i,3) > curclock,
srvstatus(i,4)= 0;%没有休息（正在忙）
else
srvstatus(i,4)= curclock-srvstatus(i,3);%目前已经空闲的时间end
end
%处理服务员服务的先后顺序（依据空闲时间）（精细处理）
tmp=srvstatus(:,4);
for i=1:numsrv,
[value,id]=max(tmp);
b(i)=id;
tmp(id)=0;%已经排序了
end
%此时等待队列必然不为空
for j=1:numsrv,
i=b(j);%确定服务员的序号
if(srvstatus(i,1)==0)
%找一个顾客开始服务，同时计算该顾客什么时候接受服务，结束服务； [m,n]=size(waiting);
if m==0,
break;
end
if waiting(1,5)==0,%还没有开始接受服务
waiting(1,3)= curclock;
waiting(1,5)= waiting(1,3)+waiting(1,4);%结束时刻
srvstatus(i,1)=1;%设置为忙状态
srvstatus(i,2)=waiting(1,1);%顾客编号
srvstatus(i,3)= waiting(1,5);%结束时刻
srvstatus(i,5)=srvstatus(i,5)+1;%又服务了一个顾客
%计算等待时间
avgwaittime(end+1) = waiting(1,3)-waiting(1,2);
disp(sprintf('间隔时间(%8.2f) 顾客编号：%5d 接受服务员（%4d）服务(到达时刻%10.2f)',waiting(1,6),waiting(1,1),i,waiting(1,2)))
endtime=max(endtime,waiting(1,5))
waiting(1,:)=[];%从等待队列中离开
end
end%if
end%for
[m,n]=size(waiting);
%计算队长（这里的计算式子可以参考排队论有关术语进行确定）
if totalcustomer < numsimucustumer
avgwaitlen(end+1)=m;
end
if sum(srvstatus(:,5))>=numsimucustumer,%队列为空，结束
finished=1;
end
end%while
disp('服务顾客数：')
disp(srvstatus(:,5)')
disp('平均队长');
disp(mean(avgwaitlen));
disp('运行时间（分钟，小时）');
disp(sprintf('%8.f%8.f',curclock,curclock/60));
disp('平均等待时间（分钟）');
disp(mean(avgwaittime ));
disp('结束时间（分钟）');
disp(endtime );
figure
hist(avgwaitlen)
title('平均队长')
figure
hist(avgwaittime)
title('平均等待时间');
1.1.8 思考题
请运行模拟程序，并分析运行结果。

实验案例 一个修理厂的模拟
1 实验案例
1.1 案例：一个修理厂的模拟 1.1.1 问题描述
某修理厂设有3个停车位置，其中一个位置供正在修理的汽车停放。

现以一天为一个时段，每天最多修好一辆车，每天到达修理站的汽车数有如下概率分布：
假定在一个时段内一辆汽车能够修好的概率为0.7，本时段内未能完成修理的汽车于正在等待修理的汽车一起进入下一时段。

试问：该停车厂有无必要增加停车位置，并说明理由。

1.1.2 模拟模型
这种排队论方面的问题采用固定时间增量法模拟。

模拟以一天为一个时段，模拟纵时间最好在1000天以上。

模拟汽车到达数量，根据概率分布：产生在[0,1]上均匀分布的随机数t ，如果6.00<≥t t 且，则认为当天到达的车辆数为0辆;如果8.06.0<≥t t 且，则认为当天到达的车辆数为1辆，如果18.0<≥t t 且，则认为当天到达的车辆数为2辆。

模拟修理情况：由于一天最多修好一辆，而一个时段内一辆汽车修好的概率为0.7，则模拟每两车的修理情况，如果这些车所能修好数目大于等于1辆，则以当天修好1辆计。

1.1.3 模拟程序
本模拟程序编写了一个主函数queue ，另外在函数queue 中编写了2个子函数： getcome ：模拟车辆到来情况，返回当天到来的车辆数目
getrepaired ：模拟修理情况，返回修好的车辆数目
整个模拟程序如下：（2005/6/6新版本）
function queue
%2005-6-6
%排队模拟主程序
%排队问题模拟
%
numdays=input('请输入模拟天数:')
numstay=0;%假定最初修理站还没有待修理的汽车
LEN=6;%定义常量
matfrequence=zeros(1,LEN);%第i个元素表示当天末还有i-1辆车在没有修好的时段频数leave_norepair=0;%存储来到，但没有停车位置而离开的车辆数
for days=1:numdays%主循环，模拟numdays个时段
temp= getcome;
if numstay + temp>3 ,
leave_norepair = leave_norepair + (numstay + temp - 3);
end
%numcome=numstay+getcome;%2004-10-10：这里有问题，受限制与停车位置数量
numcome=min(3,numstay + temp);%
%头一天还没有修好的车辆数＋当天新到来的车辆数
%numstay表示当天末还没有修理好的车辆数目
numstay=max(0,numcome - getrepaired(numcome));%
matfrequence(numstay+1)=matfrequence(numstay+1) + 1 ;
end
matfrequence
prob=matfrequence/numdays
disp(sprintf('平均每天夜里停放在修理站的车辆数=%4.2f',...
sum(matfrequence/numdays.*[0:LEN-1])))
%sprintf('=%.4f',sum(matfrequence/numdays.*[0:LEN-1]))
disp(sprintf('平均每天因位置而未修理而离开修理站的车辆数=%4.2f',...
leave_norepair/numdays))
%sprintf('=%.4f',leave_norepair/numdays)
leave_norepair
function num=getcome
%模拟车辆到来情况，返回当天到来的车辆数目
t=rand;
if t>=0 & t<0.6
num= 0;%当天到来车俩数为0辆
elseif t>=0.6 & t<0.8
num=1;%当天到来车俩数为1辆
else
num=2;%当天到来车俩数为2辆
end
function r=getrepaired(num_cur)
%模拟修理情况，返回修好的车辆数目
%n为需要修理的车辆数目
%r为n辆车修好了r辆
%num_cur 当前（天）车辆数
r=0;
if num_cur<=0,%如果根本没有车，当然就没有修好车
return
end
%只考虑当前正在修的这辆车是否能够修好
if rand<0.7,%(0,0.7) 认为修好，[0.7,1)认为没有修好
r=1;
end
1.1.4 模拟结果
程序运行结果如下：
请输入模拟天数: (100): 10000
numdays =
10000
matfrequence =
4230 2729 2375 666 0 0
prob =
0.4230 0.2729 0.2375 0.0666 0 0
平均每天夜里停放在修理站的车辆数=0.95
平均每天因位置而未修理而离开修理站的车辆数=0.09
leave_norepair =
883
模拟
最后得到平均每天夜里停放在修理站的车辆数约为0.95辆。

因此从这个结果分析可以初步认为：“不必要增加车位”或“增加车位必要性不大”。