%M/M/C模型的阻塞概率仿真
clc;clear all;close all;
N=2000;%用户数
H=30;%平均服务时间
lamda=0;%到达速率
mu=1/H;%指数分布参数
C=30;
for i=1:1:C
leave(i)=0;
end
LOSS(1)=0;
for k = 1:N;
lamda=lamda+0.01;
for i=1:1:C
leave(i)=0;
end
n=0;
%产生服从Poisson到达过程的用户到达时间
U1=rand(N,1);
temp=0;
for i=1:1:N
arrive(i)=temp-log(U1(i))/lamda;
temp=arrive(i);
end
%产生服从指数分布的用户服务时间
U2 = rand(N,1);
for i = 1:1:N
service(i) = -log(U2(i))/mu;
end
%产生用户离开时间
for i = 1:1:N
depart(i) = arrive(i)+service(i);
end
%计算阻塞概率
for i = 1 : 1 : N
flag = 0; %标志信道是否被阻塞
for j = 1 : 1 : C
if leave(j) < arrive(i)% 若果第i个用户到达的时间小于某一个用户在j个信道中离开的时间，则说明该信道空闲,可接入第i个用户。
leave(j) = depart(i);%则此时第j个信道中用户的离开时间记为depart(i)。
flag = 1;
break;
end
end
if flag == 0 % 阻塞
n=n+1;
end
end
LOSS(k)=n/N;
end
lamda_temp = 0.01 : 0.01 : N*0.01;
A = lamda_temp * H;%呼叫强度*保持时间=业务量
plot(A,LOSS,'g');%绘制呼损率的图形
hold on;
%计算阻塞概率的理论值
A = lamda_temp * H;
for j = 1 : length(A)
sum=0.0;
for i=1:1:C
temp = (A(j)^i) /factorial(i);
sum = sum + temp;
end
Pr(j) = (A(j)^C) / (factorial(C) * sum);%呼损率公式
end
plot(A, Pr, 'r');
title('呼叫过程拥塞概率的仿真结果')
xlabel('话务量(Elangs)')
ylabel('拥塞概率')
legend('仿真曲线','理论曲线')
grid on
