Comparing Numerical and Monte carlo simulation for Erlang B loss formular

I want to compare numerical evaluation(theoritical simulation) and the Monte carlo simulation for the Erlang B loss formular (E) for blocking probability. I nned someone to help me correct my codes because both results are supposed to be aligned.

; λ - Poisson arrival rate, μ - exponentially distributed service time, C- number of servers

The sellected values are: ; ;

My code for the numerical (theoritical simulation ) is :

 lambda=0.0:0.05:0.7;       %  arrival rate
 mu=1/180;                % Service rate 
 C=12;                      %No of severs
 a=lambda/mu;  

 x=(a.^C)/factorial(C);
 value=zeros(C+1,length(a));
for i=0:C
  value(i+1,:)=a.^i/factorial(i);
 end
 hold on
 M=sum(value,1);
 E=x.*M.^(-1);                     % Blocking probability in the Licensed Sprcturnm band
plot(lambda,E)

The Codes for the Monte Carlo simulation is:

% Declare and initialize the parameters
lambda=0.0:0.05:0.7;       %  arrival rate
 mu=1/180;                % Service rate 
 C=12;                      %No of severs
%Set number of iterations for Monte-Carlo simulation
Ns=1000;
E=zeros(Ns,15);
for k=1:Ns
   arv = poissrnd(lambda); % Generate random variable for Poisson arrival
   ex=exprnd(mu);        %Generate random variable for exponential service time
 % calculate E - Erlang blocking probability 
 a=arv./ex;
 x=(a.^C)/factorial(C);
 value=zeros(C+1,length(a));
for i=0:C
  value(i+1,:)=a.^i/factorial(i);
 end
 hold on
 M=sum(value,1);
 E(k+1,:)=x.*M.^(-1);    
end
M_C=mean(E);
plot(lambda,M_C)
xlabel('PU Arrival Rate \lambda');
ylabel('E')
legend({'Numerical','MC Simulation' })

The result obtained :

Answer 1

[vidyesh]

Hello Daniel,

The discrepancy in your results likely stems from generating random values for the Poisson arrival rate (λ) in your Monte Carlo simulation. For accurate comparison with the numerical (theoretical) simulation of the Erlang B formula, you should generate random values for the normalized traffic intensity (), not (λ) itself. This approach ensures that both simulations are aligned in terms of the offered load to the system, leading to comparable blocking probability results. Adjust your Monte Carlo simulation accordingly, and you should see the alignment between the two methods.

Please find the attached code below with the mentioned change.

lambda=0.0:0.05:0.7;       %  arrival rate
mu=1/180;                % Service rate 
C=12;                      %No of severs
a=lambda/mu;  

x=(a.^C)/factorial(C);
value=zeros(C+1,length(a));
for i=0:C
    value(i+1,:)=a.^i/factorial(i);
end
hold on
M=sum(value,1);
E=x.*M.^(-1);                     % Blocking probability in the Licensed Sprcturnm band
plot(lambda,E)
% Declare and initialize the parameters
%Set number of iterations for Monte-Carlo simulation
Ns=1000;
E=zeros(Ns,15);
for k=1:Ns
    %Generate random variable for exponential service time
    % calculate E - Erlang blocking probability 
    a=poissrnd (lambda./mu);
    x=(a.^C)/factorial(C);
    value=zeros(C+1,length(a));
    for i=0:C
        value(i+1,:)=a.^i/factorial(i);
    end
    hold on
    M=sum(value,1);
    E(k+1,:)=x.*M.^(-1);    
end
M_C=mean(E);
plot(lambda,M_C)
xlabel('PU Arrival Rate \lambda');
ylabel('E')
legend({'Numerical','MC Simulation' })

Hope this helps.