simulation.m

排队论中M/M/1系统的MATLAB仿真主要是阻塞率的理论值与仿真数据的比较

%==================================
%blocking probability
%=========================================================================
%顾客到达的时间间隔服从lamda指数分布；其服务时间间隔服从μ指数分布（用u表示）；
%==========================================================================
%仿真程序
clear all;
format long;
s=10;                          %服务窗口数量为10
lamda=0.1:0.1:0.9;
u=0.1;             %话务量A=lamda/u,固定u，变化lamda
k=length(lamda);
peo_num=50000;                    %测试顾客人数为50000
num_p=zeros(1,k);             %阻塞率数组
temp=zeros(1,s);                %窗口数组
n=zeros(1,k);                 %阻塞人数计数
t=1;                        
serve_time(1,:)=exprnd(1/u,1,peo_num);%顾客的服务时间间隔
for lamda=0.1:0.1:0.9;                       
    arr_interval(1,:)=exprnd(1/lamda,1,peo_num);%顾客的到达时间间隔
        for i=1:s                         %初始化服务窗口
            temp(1,i)=serve_time(1,i);
        end;
        for j=s+1:peo_num-1           %循环开始并计算阻塞人数
            min=temp(1,1);
            pos=1;
                for i=2:s         %找出服务时间最小的窗口
                    if min>temp(1,i)
                        min=temp(1,i);
                        pos=i;
                    end;
                end;      
                if arr_interval(1,j)<min    %判断是否被阻塞并计数
                    n(1,t)=n(1,t)+1;              
                    arr_interval(1,j+1)=arr_interval(1,j)+arr_interval(1,j+1)-min;%更新下一顾客的达到时间其中涵盖了顾客的等待时间
                    flag=1; 
                else
                    flag=0;
                end;
                for jj=1:s               %更新窗口中的服务时间
                       if flag==0   
                            temp(1,jj)=temp(1,jj)-arr_interval(1,j);
                       else
                            temp(1,jj)=temp(1,jj)-min;  
                       end;
                       if temp(1,jj)<0
                            temp(1,jj)=0;          
                       end;
                end;                                   
             temp(1,pos)=serve_time(1,j); %更新新服务窗口的服务时间
        end;
     num_p(1,t)=n(1,t)/peo_num; %计算阻塞率
     t=t+1;
end;
%====================================================
%理论推导(erlang_C)
A=1:9;                  %A=λ/u话务量,S窗口数量
block_p=zeros(1,length(A));   %block_p理论阻塞率数组
for i=1:length(A)
    v=0;
    m=0;
        for j=0:s-1                 %分母上的求和公式
            m=m+power(i,j)/factorial(j);
        end;
      v=power(i,s)+factorial(s)*(1-i./s)*m;  %Erlang_C的分母部分
      block_p(1,i)=power(i,s)/v;
end;
semilogy(1:k,num_p,'r-');
hold on;
semilogy(1:k,block_p,'p-');
grid on;
title('M/M/m排队模型阻塞率的理论值和仿真值比较');
legend('仿真值','理论值');
xlabel('话务量');ylabel('阻塞率');