系统模拟实验的三个案例

1、系统模拟实验的三个案例实验案例   赶上火车的概率    1实验案例 1.1    赶上火车的概率 1.1.1     问题描述 如图，一列火车从A站开往B站，某人每天赶往B站上这趟火车。他已了解到： （1）   火车从A站到B站的运行时间是均值为30分钟，标准差为2分钟的随机变量； （2）   火车在下午大约1点离开A站，离开时刻的频率分布如下：   出发时刻 午后1:00 午后1:05 午后1:10 频率 0.7 0.2 0.1 &

2、#160; 此人到达B站的时刻频率分布为： 时刻 午后1:28 午后1:30 午后1:32 午后1:34 频率 0.3 0.4 0.2 0.1 问他能赶上火车的概率是多少？   1.1.2    变量说明 ：火车从A站出发的时刻； ：火车从A站到B站的运行时间；单位：分钟 ：他到达B站的时刻 1.1.3    问题分析与假设 此问题包含多个随机因素。 这里假设，都是随机变量，其中服从正态分布。 1.1.4    模型建立 很显然，他能及时赶上火车的条件是：。为了简化计算，将下午1点记为初始时刻。

3、和的分布律如下： /min 0 5 10 0.7 0.2 0.1   /min 28 30 32 34 0.3 0.4 0.2 0.1 为了模拟随机变量。如果为在均匀分布的随机数，为了模拟随机变量，可以通过如下方法。 ，。 其中，和分别用来模拟随机变量和。   1.1.5    模拟算法 变量说明： k    临时变量，存储当前累计模拟次数 count  存储赶上火车的次数   第 1 步       输入模拟次数n 第 2 步

4、0;      k=1，count=0 第 3 步       当k<=n，执行第4步，否则执行第12步 第 4 步       生成均匀分布随机数赋给r 第 5 步       由r及公式确定T1模拟火车出发时刻 第 6 步       生成均匀分布随机数赋给r； 第 7 步  &

5、#160;    由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

6、.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       

7、60; T1=5;    else         T1=10;    end       %模拟随机变量t2（火车的运行时间）      T2=30+randn*2;         %模拟随机变量t3（他到达B站的时刻）      rt3=rand;    if rt3<0.3&

8、#160;        T3=28;    elseif rt3>=0.3 & rt3<0.7          T3=30;    elseif rt3>=0.7 & rt3<0.9          T3=32;    else    

9、0;       T3=34;    end      if T3 < T1 + T2,%赶上了            count=count+1;      end end%for prob=count/total   1.1.7    模拟结果   命令行中输入以下语句： sim_trai

10、n   运行结果输出：   输入模拟次数：100   prob =   0.6302 此次运行结果显示赶上火车的近似概率为0.6左右。下面列表给出多次运行模拟程序的结果。   序号 模拟次数 近似概率p 1 500 0.6280 2 500 0.6920 3 1000 0.6530 4 1000 0.6490 5 5000 0.6260 6 5000 0.6288   1.1.8    评价与改进方向 为了计算赶上火车的概率，本文采用了随机系统模拟的方法。如果能够从模型出发，对赶上火车的概率进行近似计算，

11、然后与模拟结果进行对比，这样模拟会更有说明力。 1.1.9    思考题 （1）   请思考用其它方法计算赶上火车的概率或近似概率。 （2）   如果要使得他赶上火车的概率大于95%，你有什么办法？结合上面的数学模型及模拟程序来思考。 （3）   通过该问题的建模求解，你能归纳出一般系统模拟的方法步骤么？   实验案例       理发店模拟     1     &#

12、160; 实验案例 1.1    案例：理发店模拟 例子：一个理发店有两位服务员A和B，顾客们随机到达店内，其中60的顾客仅需剪发，每位花5分钟时间，另外40顾客既要剪发又要洗发，每位用时8分钟。 理发店是个含有多种随机因素的系统，请对该系统进行模拟，并对其进行评判。（准备怎么做）可供参考内容“排队论”，“系统模拟”，“离散系统模拟”，“事件调度法” 1.1.1    问题分析 理发店系统包含诸多随机因素，为了对其进行评判就是要研究其运行效率，从理发店自身利益来说，要看服务员工作负荷是否合理，是否需要增加员工等考虑。从顾客角度讲，还要

13、看顾客的等待时间，顾客的等待队长，如等待时间过长或者等待的人过多，则顾客会离开。理发店系统是一个典型的排队系统，可以用排队论有关知识来研究。  1.1.2    模型假设 1            60的顾客只需剪发，40的顾客既要剪发，又要洗发； 2            每个服务员剪发需要的时间均为5分钟，既剪发又洗发则花8分钟； 3

14、            顾客的到达间隔时间服从指数分布；4            服务中服务员不休息。1.1.3    变量说明 u：剪发时间（单位：分钟），u=5m； v:  既剪发又理发花的时间（单位：分钟），v=8m； T： 顾客到达的间隔时间，是随机变量，服从参数为 的指数分布，（单位：分钟）T0：顾客到达的平均间隔时间（单位：秒

15、），T0 ；  1.1.4    模型建立 由于该系统包含诸多随机因素，很难给出解析的结果，因此可以借助计算机模拟对该系统进行模拟。考虑一般理发店的工作模式，一般是上午9：00开始营业，晚上10：00左右结束，且一般是连续工作的，因此一般营业时间为13小时左右。 这里以每天运行12小时为例，进行模拟。 这里假定顾客到达的平均间隔时间T0服从均值3分钟的指数分布， 则有3小时到达人数约为 人，6小时到达人数约为 人，10小时到达人数约为 人，这里模拟顾客到达数为60人的情况。 （如何选择模拟的总人数或模拟总时间）1.1.5   

16、; 系统模拟 根据系统模拟的一般方法，需要考虑系统的如下数据、参数。1.        状态（变量）    （1）     等待服务的顾客数；（2）     A 是否正在服务； （3）     B 是否正在服务； 2.        实体：两名服务员、顾客们3.      

17、60; 事件： （1）     一名新顾客的到达；（2）     A 开始服务； （3）     A 结束服务； （4）     B 开始服务； （5）     B 结束服务； 4.        活动：（1）     顾客排队时间（2）     顾客们到达的间隔时间（3）     A 的

18、服务时间 （4）     B 的服务时间；  在系统模拟时，为了研究系统的整体情况，这里考虑顾客到达后不离开，且等待队长不限。要考虑如果服务员均空闲时，顾客先选择谁服务？要考虑模拟的时间设置还有顾客数目。模拟终止条件是根据顾客数目还是根据营业时间终止？ 1.1.6    计算机模拟算法设计 自行设计finished=0; 初始化运行时钟 while finished=0              if &#

19、160; 产生的顾客数不到规定数目时 then， 产生该顾客的有关数据； 将顾客加入等待队列； else 运行时钟继续； endif 处理服务员的状态（包括工作状态，空闲时间）； 获得服务员的服务优先顺序； 根据服务员优先顺序从等待队列中安排服务； endwhile 有无参考算法？  离散系统仿真算法：事件调度法 1.1.7    计算机模拟程序 顾客到达的间隔时间T的计算机产生方法，利用T=，   %理发店系统的模拟（案例分析之一） %关键词：面向事件的计算机模拟技术 clear all   curclock=0;%当前时刻，动态变

20、化   totalcustomer=0;%总共服务的顾客数 numsrv=2; srvstatus=zeros(numsrv,5);%服务员有关数据 %srvstatus 第1列：服务状态（0空闲，1正在服务）；第2列：当前服务顾客编号； %   第3列：当前服务结束时刻；第4列：服务员空闲时间；第5列：服务的顾客总数   endtime =0;%结束时间 waiting=;%等待队列数据 %waiting 第1列：顾客编号；第2列：顾客到达时刻；第3列：顾客开始接受服务时刻； %      第4列：接受服务时间；

21、第5列：顾客结束服务时刻；第6列：间隔时间 cur=zeros(1,6);%当前产生顾客的数据，对应关系同waiting avgwaitlen=;%平均等待队长 avgwaittime =;% 平均等待时间 ujiange=5;%平均间隔时间 finished=0; numsimucustumer=yesinput('输入等待模拟的顾客数：',10,10 1000); while   finished=0,       if totalcustomer < numsimucustumer   

22、    %产生一个顾客的到达及其有关性质的数据            totalcustomer = totalcustomer+1;            jiange= -log(rand)*ujiange;%与上一个顾客的到达的间隔时间       curclock = curclock + jiange; &#

23、160;     cur(1)= totalcustomer ;% 第1列：顾客编号       cur(2) = curclock;%第2列：顾客到达时刻 cur(6) = jiange; 第6列：间隔时间 %下面产生接受服务时间（可改进模型）       if rand<0.6, %产生顾客有关性质：这里是产生接受服务时间          cur(4) = 5; &

24、#160;     else          cur(4) = 8;       end       %放入等待队列       if isempty(waiting),          waiting= cur;   

25、60;   else          m,n=size(waiting);          waiting(m+1,:)= cur;       end    else       curclock = curclock + (-log(rand)*ujiange);   

26、end%if totalcustomer<          %分配等待队列（看是否有服务员空闲，如果有则分配；否则继续执行）     %处理服务员的服务状态    for i=1:numsrv,       if srvstatus(i,1)=1 & srvstatus(i,3) <= curclock,          srvstatus(i

27、,1)=0;%设置为空闲状态          srvstatus(i,4)= curclock-srvstatus(i,3);%目前已经空闲的时间       elseif srvstatus(i,1)=1 & srvstatus(i,3) > curclock,          srvstatus(i,4)= 0;%没有休息（正在忙）   

28、0;   else          srvstatus(i,4)= curclock-srvstatus(i,3);%目前已经空闲的时间       end         end    %处理服务员服务的先后顺序（依据空闲时间）（精细处理）    tmp=srvstatus(:,4);    for i=1:num

29、srv,       value,id=max(tmp);       b(i)=id;       tmp(id)=0;%已经排序了         end             %此时等待队列必然不为空    for j=1:numsrv,   

30、;    i=b(j);%确定服务员的序号       if(srvstatus(i,1)=0)          %找一个顾客开始服务，同时计算该顾客什么时候接受服务，结束服务；          m,n=size(waiting);          if m=0, 

31、0;           break;          end                        if waiting(1,5)=0,%还没有开始接受服务    &#

32、160;        waiting(1,3)= curclock;             waiting(1,5)= waiting(1,3)+waiting(1,4);%结束时刻             srvstatus(i,1)=1;%设置为忙状态     

33、        srvstatus(i,2)=waiting(1,1);%顾客编号             srvstatus(i,3)= waiting(1,5);%结束时刻             srvstatus(i,5)=srvstatus(i,5)+1;%又服务了一个顾客    

34、;           %计算等待时间             avgwaittime(end+1) = waiting(1,3)-waiting(1,2);                    &

35、#160;   disp(sprintf('间隔时间(%8.2f) 顾客编号：%5d 接受服务员（%4d）服务(到达时刻%10.2f)',waiting(1,6),waiting(1,1),i,waiting(1,2)             endtime=max(endtime,waiting(1,5)             waiti

36、ng(1,:)=;%从等待队列中离开                      end       end%if         end%for        m,n=size(waiting);    %计算队长（这里

37、的计算式子可以参考排队论有关术语进行确定）    if totalcustomer < numsimucustumer       avgwaitlen(end+1)=m;    end    if sum(srvstatus(:,5)>=numsimucustumer,%队列为空，结束       finished=1;    end      end%while 

38、0; 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 his

39、t(avgwaitlen)   title('平均队长')   figure hist(avgwaittime) title('平均等待时间');    1.1.8    思考题 请运行模拟程序，并分析运行结果。 实验案例 一个修理厂的模拟     1       实验案例   1.1    案例：一个修理厂的模拟 1.1.1  

40、60; 问题描述 某修理厂设有3个停车位置，其中一个位置供正在修理的汽车停放。现以一天为一个时段，每天最多修好一辆车，每天到达修理站的汽车数有如下概率分布：   到达数 0 1 2 概率 06 02 02 假定在一个时段内一辆汽车能够修好的概率为0.7，本时段内未能完成修理的汽车于正在等待修理的汽车一起进入下一时段。试问：该停车厂有无必要增加停车位置，并说明理由。   1.1.2    模拟模型        这种排队论方面的问题采用固定时间增量法模拟。模拟以一天为一个时段，模拟纵时间

41、最好在1000天以上。        模拟汽车到达数量，根据概率分布：产生在0,1上均匀分布的随机数t，如果，则认为当天到达的车辆数为0辆;如果，则认为当天到达的车辆数为1辆，如果，则认为当天到达的车辆数为2辆。        模拟修理情况：由于一天最多修好一辆，而一个时段内一辆汽车修好的概率为0.7，则模拟每两车的修理情况，如果这些车所能修好数目大于等于1辆，则以当天修好1辆计。 1.1.3    模拟程序    

42、    本模拟程序编写了一个主函数queue，另外在函数queue中编写了2个子函数： getcome：模拟车辆到来情况，返回当天到来的车辆数目 getrepaired：模拟修理情况，返回修好的车辆数目 整个模拟程序如下：（2005/6/6新版本） function queue %2005-6-6 %排队模拟主程序 %排队问题模拟 % numdays=input('请输入模拟天数:') numstay=0;%假定最初修理站还没有待修理的汽车   LEN=6;%定义常量 matfrequence=zeros(1,LEN);%第i个元素表示当

43、天末还有i-1辆车在没有修好的时段频数   leave_norepair=0;%存储来到，但没有停车位置而离开的车辆数   for days=1:numdays%主循环，模拟numdays个时段    temp= getcome;    if numstay + temp>3 ,       leave_norepair = leave_norepair + (numstay + temp - 3);    end     

44、60; %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