有点标题党了,主要想介绍一下M/M/1排队系统（又叫单服务台排队系统）。意思是说只有一个服务员，服务员空闲的时候，到来一个顾客，那么该顾客就接受服务，否则，到来的顾客要排队等候，等待服务。 所以说嘛，类比下，女孩相当于服务员，追求者相当于顾客。。下面的MM1系统仿真是我大学的时候上计算机仿真课的时候编写的一个程序，核心思想有点类似于“微元法”。即动态再瞬时之中可以看成静态。
由于时间久远，我已经找不到当初写的实验报告和开发文档，就直接把代码贴上面吧，用的MATLAB。会把这篇文章发到首页。我知道会被撤下来，也会挨骂。但是我这么做主要是想试探下园子里有多少人用过MATLAB，有多少人是在读硕士，再做研究，而不是搞项目开发。
最后说一句：MATLAB对于通信类，计算机算法研究累的研究生来说，应该是一个非常奇妙的工具。它以矩阵处理和数学运算见长。正像他的开发者所言“一张非常好用的电子草纸”。利用MATLAB实现一些常规算法，要远比C，C++，Java这些来的快。
下面开始贴代码啦。
计算队长的函数
function AveQuLen=CountQueueLen(EventQueue,SimCurNo)
%计算系统的平均等待队长
m=find(EventQueue(14,:)==SimCurNo);%由于数据表的最后几列有异变
k=m(length(m));%找出最后一个已服务人数为SimCurNo的数据项
TotalPeriod=EventQueue(1,k);%该项的时间为系统的总的仿真时间
TimeSpan=EventQueue(1,1:k);%时间跨度
sum=0;
for jjj=2:k-1
sum=sum+EventQueue(7,jjj)*(EventQueue(1,jjj+1)-EventQueue(1,jjj));
end;
figure(1)
set(1,'position',[10,50,200,200]);
plot(TimeSpan,EventQueue(7,1:k));
title('系统队长分布');
xlabel('时间（min）');
ylabel('队长');
come=find(EventQueue(2,:)==1);%从事件记录表中找出到达事件
ReachTime=EventQueue(4,come);
Number=EventQueue(3,come);
DepartTime=EventQueue(12,come);
figure(2)
set(2,'position',[350,50,200,200]);
plot(Number,ReachTime,'r',Number,DepartTime,'b');
text(500,170,'（注：红色代表到达时间，蓝色代表离去时间）')
title('各顾客到达时间和离去时间')
StayTime=EventQueue(13,come);
WaitTime=EventQueue(10,come);
figure(3)
set(3,'position',[600,50,200,200]);
plot(Number,StayTime,'r',Number,WaitTime,'b');
title('各顾客在系统中的停留时间和等待时间');
text(700,100,'(注：红色代表停留时间，蓝色代表等待时间)')
AveQuLen=sum/TotalPeriod;
%计算系统的平均等待队长
m=find(EventQueue(14,:)==SimCurNo);%由于数据表的最后几列有异变
k=m(length(m));%找出最后一个已服务人数为SimCurNo的数据项
TotalPeriod=EventQueue(1,k);%该项的时间为系统的总的仿真时间
TimeSpan=EventQueue(1,1:k);%时间跨度
sum=0;
for jjj=2:k-1
sum=sum+EventQueue(7,jjj)*(EventQueue(1,jjj+1)-EventQueue(1,jjj));
end;
figure(1)
set(1,'position',[10,50,200,200]);
plot(TimeSpan,EventQueue(7,1:k));
title('系统队长分布');
xlabel('时间（min）');
ylabel('队长');
come=find(EventQueue(2,:)==1);%从事件记录表中找出到达事件
ReachTime=EventQueue(4,come);
Number=EventQueue(3,come);
DepartTime=EventQueue(12,come);
figure(2)
set(2,'position',[350,50,200,200]);
plot(Number,ReachTime,'r',Number,DepartTime,'b');
text(500,170,'（注：红色代表到达时间，蓝色代表离去时间）')
title('各顾客到达时间和离去时间')
StayTime=EventQueue(13,come);
WaitTime=EventQueue(10,come);
figure(3)
set(3,'position',[600,50,200,200]);
plot(Number,StayTime,'r',Number,WaitTime,'b');
title('各顾客在系统中的停留时间和等待时间');
text(700,100,'(注：红色代表停留时间，蓝色代表等待时间)')
AveQuLen=sum/TotalPeriod;
计算等待时间的函数
代码
unction AveWaitTime=CountWaitingTime(EventQueue,SimCusNo)
sum=0;
k=size(EventQueue);
for kk=2:k(2)
if EventQueue(2,kk)==1
sum=sum+EventQueue(10,kk);
end;
end;
AveWaitTime=sum/SimCusNo;
sum=0;
k=size(EventQueue);
for kk=2:k(2)
if EventQueue(2,kk)==1
sum=sum+EventQueue(10,kk);
end;
end;
AveWaitTime=sum/SimCusNo;
保存信息的函数
代码
function AssistArray=Reserve(EventQueue,t_CurTime,t_ReachTime);
j=1;
AssistArray=[];
k=size(EventQueue);
%将要处理的离去事件的下标保存到辅助数组AssistArray中，待下一步处理
for ii=k(2):-1:2
if EventQueue(2,ii)==2%判断是否是离开事件
if EventQueue(8,ii)~=Inf
break;%如果已经更新，则代表前面的离开事件都已更新，退出循环
end;
if (EventQueue(1,ii)>=t_CurTime)&(EventQueue(1,ii)<=t_ReachTime)%选出在指定范围内的事件
AssistArray(j)=ii;
j=j+1;
end;
end;
end;
j=1;
AssistArray=[];
k=size(EventQueue);
%将要处理的离去事件的下标保存到辅助数组AssistArray中，待下一步处理
for ii=k(2):-1:2
if EventQueue(2,ii)==2%判断是否是离开事件
if EventQueue(8,ii)~=Inf
break;%如果已经更新，则代表前面的离开事件都已更新，退出循环
end;
if (EventQueue(1,ii)>=t_CurTime)&(EventQueue(1,ii)<=t_ReachTime)%选出在指定范围内的事件
AssistArray(j)=ii;
j=j+1;
end;
end;
end;
负责调度的函数
clear all
close all
%MM1排队系统计算机仿真
SimCusNo=input('请设定MM1仿真系统的被服务顾客数SimCusNo=[1000]');
if isempty(SimCusNo)
SimCusNo=1000;
end;
%声明全局变量并将全局变量初始化
Ba=5;%到达时间间隔的均值
Bs=4;%服务时间的均值
EventQueue=[];%事件记录列表，它是一张表即二维矩阵，每一个列向量是一个事件
IfOutOfNum = 0; %是否超出指定的仿真人数范围
NextTime = 0; %下一到达时间
CustomerNum = 0; %系统中顾客数
EndTime = 0; %离去时间
HasServedNum = 0; %已服务人数
CurStatus =0; %服务台状态
i=1; %事件记录表的指针
EventInit=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%Event表示事件记录表中的一项完整记录具有十四项属性Event(1)表示仿真钟；Event（2）表示事件类型，共有两类事件1表到达事件，2表离开事
%件；Event（3）表示顾客编号；Event（4）表示到达时间对于到达事件它的值和仿真钟一样；Event(5)表示下一事件的到达时间，它服从负指数分
%布；Event(6)表示服务台状态0表示空闲，1表示繁忙，它的值取决于系统中的顾客数；Event(7)表示队长，它的值为系统中的顾客数减一;Even
%t(8)表示系统中的顾客数；Event(9)表示该顾客的服务开始时间;Event(10)表示该顾客的等待时间；Event(11)表示该顾客的服务时间它也服从负
%指数分布；Event(12)表示该顾客的离去时间；Event(13)表示该顾客的逗留时间；Event（14）表示系统的已服务人数
EventQueue(:,i)=EventInit;%事件纪录表初始化
i=i+1;%事件记录表指针后移
while (IfOutOfNum==0)
if (HasServedNum >=SimCusNo)
IfOutOfNum = 1; %超出仿真范围
end ;
%if CustomerNum==SimCusNo;
% IfOutOfNum=1;
% end;
% ****************************首先产生到达事件*************************************
t_CurTime = NextTime; %当前仿真钟
t_ReachTime = t_CurTime+(-Ba*log(rand())) ; %产生随机到达时间
%求顾客的开始服务时间
if EndTime<t_CurTime
t_StartTime=t_CurTime;%如果他到来时他的前一个顾客已经离开，他到来后就能开始服务
else
t_StartTime=EndTime;%如果他到来时他的前一个顾客还在服务中，那么他只能在前面的人走后才能开始服务
end ;
%求顾客的等待时间
if EndTime-NextTime>0
t_WaitTime=EndTime-NextTime;%如果他到来的时候前面的顾客没有离开，则他需要等待
else
t_WaitTime=0;%如果他到来的时候他前面的顾客已经离开则他的等待时间为0；
end;
%产生随机服务时间
t_ServeTime = -Bs*log(rand());
%当前服务台状态
if CustomerNum>0 %如果顾客到来时原系统中的顾客数大于零则表示当前服务台处于繁忙状态
CurStatus=1;
else
CurStatus=0;%反之则表示当前服务台处于空闲状态
end
%为到来事件属性赋值
Come_Event=[t_CurTime,... %1:仿真钟
1,... % 2:事件类型 1-到达
CustomerNum+HasServedNum+1,...% 3:顾客编号= 系统中的顾客数+已服务的顾客数+1
t_CurTime,...%4:到达时间 = 仿真钟
t_ReachTime,...%5:下一到达时间 = 到达时间+随机间隔时间
CurStatus,...%6:服务台状态
CustomerNum,...%7:队长=该顾客到来时原系统中的顾客数
CustomerNum+1,...%8:当前系统中的顾客数（包括新到来的顾客）
t_StartTime,...%9:服务开始时间 = 上一顾客离去时间 或 到达时间(当离去时间<到达时间)
t_WaitTime,...%10:等待时间 = 上一顾客离去时间-到达时间或为0
t_ServeTime,...%11:服务时间 = 随机时间长度
t_StartTime+t_ServeTime,...%12:离去时间 = 服务开始时间+服务时间
t_WaitTime+t_ServeTime,...%13:逗留时间 = 等待时间+服务时间
HasServedNum...%14:已服务人数
];
%*******************将事件投入事件记录表中，并和它前面的事件发生时间进行比较使系统的仿真钟向前推移***********
EventQueue(:,i)=Come_Event;
%找出它的前一到达事件（通过顾客编号查找）
trace=find(EventQueue(2,1:i-1)==1);%找出它前面的所有到达事件
%如果该事件是第一个到达事件及该事件前面没有其他事件setpos为一空矩阵
if length(trace)>0 %如果该到达事件不是第一个到达事件
setpos=trace(length(trace));%找到它的前一个到达事件
%min_w=find(min(EventQueue(1,setpos+1:i-1)));
for mm=setpos+1:i-1%用选择法进行排序
min_s=mm;
for(nn=mm+1:i)
if(EventQueue(1,nn)<EventQueue(1,min_s))
min_s=nn;
end;
if(min_s~=mm)
temp=EventQueue(:,min_s);
EventQueue(:,min_s)=EventQueue(:,mm);
EventQueue(:,mm)=temp;
end;
end;
end;
end;
i=i+1;% 事件记录表指针后移
%************************************************产生此事件对应的离开事件************
Leave_Event=[Come_Event(12),...%1:仿真钟 = 离去时间
2,... %2:事件类型 2-离开
Come_Event(3),...%3:顾客编号
Inf,...%4：到达时间
Inf,...%5:下一到达时间
Inf,...%6:服务台状态,在下面的模块中会更新
Inf,...%7:队长 在下面的模块中稍候会更新
Inf,...%8：系统中的顾客数，稍候作更新
Inf,...%9:服务开始时间
Inf,...%10：等待时间
Inf,...%11:服务时间
Come_Event(12),...%12：离去时间=当前仿真钟
Inf,...%13:逗留时间
Inf...%14:已服务人数，稍候作更新
];
%*******将此事件写入事件记录表****
EventQueue(:,i)=Leave_Event;
i=i+1;
%*********更新系统全局变量***************************************
NextTime = t_ReachTime; %下一顾客到达时间
CustomerNum=CustomerNum+1; %系统中顾客数加 1
EndTime = t_StartTime+t_ServeTime; %离去时间
%*****更新更新当前到达事件与下一到达事件间的所有离开事件 (系统中的顾客数, 已服务人数等)**********
%size(EventQueue);
AssistArray=Reserve(EventQueue,t_CurTime,t_ReachTime);
%对选中的离去事件进行更新
for jj=length(AssistArray):-1:1
%对系统全局变量进行更新
CustomerNum=CustomerNum-1;%系统中顾客数减 1
HasServedNum=HasServedNum+1;%已服务人数加 1
%对事件记录表中的离开事件项属性进行更新
EventQueue(8,AssistArray(jj))=CustomerNum;%对此时的系统顾客数进行更新
EventQueue(14,AssistArray(jj))=HasServedNum;%对已服务顾客数进行更新
if CustomerNum>0
CurStatus=1;
else
CurStatus=0;
end ;
%tmp=0;%辅助变量使用之前应把初值置零以免产生干扰
EventQueue(6,AssistArray(jj))=CurStatus;%服务台状态
tmp=CustomerNum;
if tmp-1>0
EventQueue(7,AssistArray(jj))=tmp-1;
else
EventQueue(7,AssistArray(jj))=0;
end;
end;
end;
%计算系统的平均等待时长
disp('平均等待时间如下')
AveWaitTime=CountWaitingTime(EventQueue,SimCusNo)
disp('平均等待队长如下')
AveQuLen=CountQueueLen(EventQueue,SimCusNo)
disp('仿真结束');
close all
%MM1排队系统计算机仿真
SimCusNo=input('请设定MM1仿真系统的被服务顾客数SimCusNo=[1000]');
if isempty(SimCusNo)
SimCusNo=1000;
end;
%声明全局变量并将全局变量初始化
Ba=5;%到达时间间隔的均值
Bs=4;%服务时间的均值
EventQueue=[];%事件记录列表，它是一张表即二维矩阵，每一个列向量是一个事件
IfOutOfNum = 0; %是否超出指定的仿真人数范围
NextTime = 0; %下一到达时间
CustomerNum = 0; %系统中顾客数
EndTime = 0; %离去时间
HasServedNum = 0; %已服务人数
CurStatus =0; %服务台状态
i=1; %事件记录表的指针
EventInit=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%Event表示事件记录表中的一项完整记录具有十四项属性Event(1)表示仿真钟；Event（2）表示事件类型，共有两类事件1表到达事件，2表离开事
%件；Event（3）表示顾客编号；Event（4）表示到达时间对于到达事件它的值和仿真钟一样；Event(5)表示下一事件的到达时间，它服从负指数分
%布；Event(6)表示服务台状态0表示空闲，1表示繁忙，它的值取决于系统中的顾客数；Event(7)表示队长，它的值为系统中的顾客数减一;Even
%t(8)表示系统中的顾客数；Event(9)表示该顾客的服务开始时间;Event(10)表示该顾客的等待时间；Event(11)表示该顾客的服务时间它也服从负
%指数分布；Event(12)表示该顾客的离去时间；Event(13)表示该顾客的逗留时间；Event（14）表示系统的已服务人数
EventQueue(:,i)=EventInit;%事件纪录表初始化
i=i+1;%事件记录表指针后移
while (IfOutOfNum==0)
if (HasServedNum >=SimCusNo)
IfOutOfNum = 1; %超出仿真范围
end ;
%if CustomerNum==SimCusNo;
% IfOutOfNum=1;
% end;
% ****************************首先产生到达事件*************************************
t_CurTime = NextTime; %当前仿真钟
t_ReachTime = t_CurTime+(-Ba*log(rand())) ; %产生随机到达时间
%求顾客的开始服务时间
if EndTime<t_CurTime
t_StartTime=t_CurTime;%如果他到来时他的前一个顾客已经离开，他到来后就能开始服务
else
t_StartTime=EndTime;%如果他到来时他的前一个顾客还在服务中，那么他只能在前面的人走后才能开始服务
end ;
%求顾客的等待时间
if EndTime-NextTime>0
t_WaitTime=EndTime-NextTime;%如果他到来的时候前面的顾客没有离开，则他需要等待
else
t_WaitTime=0;%如果他到来的时候他前面的顾客已经离开则他的等待时间为0；
end;
%产生随机服务时间
t_ServeTime = -Bs*log(rand());
%当前服务台状态
if CustomerNum>0 %如果顾客到来时原系统中的顾客数大于零则表示当前服务台处于繁忙状态
CurStatus=1;
else
CurStatus=0;%反之则表示当前服务台处于空闲状态
end
%为到来事件属性赋值
Come_Event=[t_CurTime,... %1:仿真钟
1,... % 2:事件类型 1-到达
CustomerNum+HasServedNum+1,...% 3:顾客编号= 系统中的顾客数+已服务的顾客数+1
t_CurTime,...%4:到达时间 = 仿真钟
t_ReachTime,...%5:下一到达时间 = 到达时间+随机间隔时间
CurStatus,...%6:服务台状态
CustomerNum,...%7:队长=该顾客到来时原系统中的顾客数
CustomerNum+1,...%8:当前系统中的顾客数（包括新到来的顾客）
t_StartTime,...%9:服务开始时间 = 上一顾客离去时间 或 到达时间(当离去时间<到达时间)
t_WaitTime,...%10:等待时间 = 上一顾客离去时间-到达时间或为0
t_ServeTime,...%11:服务时间 = 随机时间长度
t_StartTime+t_ServeTime,...%12:离去时间 = 服务开始时间+服务时间
t_WaitTime+t_ServeTime,...%13:逗留时间 = 等待时间+服务时间
HasServedNum...%14:已服务人数
];
%*******************将事件投入事件记录表中，并和它前面的事件发生时间进行比较使系统的仿真钟向前推移***********
EventQueue(:,i)=Come_Event;
%找出它的前一到达事件（通过顾客编号查找）
trace=find(EventQueue(2,1:i-1)==1);%找出它前面的所有到达事件
%如果该事件是第一个到达事件及该事件前面没有其他事件setpos为一空矩阵
if length(trace)>0 %如果该到达事件不是第一个到达事件
setpos=trace(length(trace));%找到它的前一个到达事件
%min_w=find(min(EventQueue(1,setpos+1:i-1)));
for mm=setpos+1:i-1%用选择法进行排序
min_s=mm;
for(nn=mm+1:i)
if(EventQueue(1,nn)<EventQueue(1,min_s))
min_s=nn;
end;
if(min_s~=mm)
temp=EventQueue(:,min_s);
EventQueue(:,min_s)=EventQueue(:,mm);
EventQueue(:,mm)=temp;
end;
end;
end;
end;
i=i+1;% 事件记录表指针后移
%************************************************产生此事件对应的离开事件************
Leave_Event=[Come_Event(12),...%1:仿真钟 = 离去时间
2,... %2:事件类型 2-离开
Come_Event(3),...%3:顾客编号
Inf,...%4：到达时间
Inf,...%5:下一到达时间
Inf,...%6:服务台状态,在下面的模块中会更新
Inf,...%7:队长 在下面的模块中稍候会更新
Inf,...%8：系统中的顾客数，稍候作更新
Inf,...%9:服务开始时间
Inf,...%10：等待时间
Inf,...%11:服务时间
Come_Event(12),...%12：离去时间=当前仿真钟
Inf,...%13:逗留时间
Inf...%14:已服务人数，稍候作更新
];
%*******将此事件写入事件记录表****
EventQueue(:,i)=Leave_Event;
i=i+1;
%*********更新系统全局变量***************************************
NextTime = t_ReachTime; %下一顾客到达时间
CustomerNum=CustomerNum+1; %系统中顾客数加 1
EndTime = t_StartTime+t_ServeTime; %离去时间
%*****更新更新当前到达事件与下一到达事件间的所有离开事件 (系统中的顾客数, 已服务人数等)**********
%size(EventQueue);
AssistArray=Reserve(EventQueue,t_CurTime,t_ReachTime);
%对选中的离去事件进行更新
for jj=length(AssistArray):-1:1
%对系统全局变量进行更新
CustomerNum=CustomerNum-1;%系统中顾客数减 1
HasServedNum=HasServedNum+1;%已服务人数加 1
%对事件记录表中的离开事件项属性进行更新
EventQueue(8,AssistArray(jj))=CustomerNum;%对此时的系统顾客数进行更新
EventQueue(14,AssistArray(jj))=HasServedNum;%对已服务顾客数进行更新
if CustomerNum>0
CurStatus=1;
else
CurStatus=0;
end ;
%tmp=0;%辅助变量使用之前应把初值置零以免产生干扰
EventQueue(6,AssistArray(jj))=CurStatus;%服务台状态
tmp=CustomerNum;
if tmp-1>0
EventQueue(7,AssistArray(jj))=tmp-1;
else
EventQueue(7,AssistArray(jj))=0;
end;
end;
end;
%计算系统的平均等待时长
disp('平均等待时间如下')
AveWaitTime=CountWaitingTime(EventQueue,SimCusNo)
disp('平均等待队长如下')
AveQuLen=CountQueueLen(EventQueue,SimCusNo)
disp('仿真结束');
由于时间久远，大学毕业我就把实验报告和文档删除了。仅留了一份代码，现在看来：其实应该都保存才对的。如果带进对MM1系统感兴趣，就自己Google些理论吧。另外如果正在做这个程序的XDXM需要源码的话可以到http://download.csdn.net/source/1030295，进行下载。
等什么时候有时间，在把MATLAB版的prim和kruskal算法贴上来。由于和Matlab语言结合的紧切，使得实现非常有特色。如果您正在学习MATLAB。希望我的代码能对你的学习带来启迪。
大家轻点骂我吧。可能对一大部分人来讲，我这篇文章是没有用的，对大学在读的XDXM做实验还是有些帮助的。
