INFO4-EP TD-MMKK

This code simulates the continuous time Markov chain associated to a M/M/K/K queue and prints the proportion of time where all servers are occupied (the blocking probability). This is compared with the one obtained by the Erlang B formula.

knitr::opts_chunk$set(echo = TRUE)

set.seed(17); # Set seed for reproducibility
N=1e5; # number of events to simulate;

K=10;         # number of servers
lambda=0.6*K; # arrival rate
mu=1;         # service rate

p_blocking=0; # blocking probability
p_blocking_transient=rep(N, 0); # transient blocking probability
T=0; # simulation time

state=0;

for (i in 1:N) {

  U=lambda+state*mu;
  time_in_state=rexp(1,U)
  T=T+time_in_state;

  if (state==K) {
    # For the moment, let us find the overall time spent in state K
    p_blocking=p_blocking+time_in_state;
  }

  p_blocking_transient[i]=p_blocking/T;

  if (runif(1)<lambda/U){
      # An arrival occurred
      if (state<K){
        # The incoming job has room so it is accepted
        state=state+1;
      }
  } else {
      # A departure occurred
      state=state-1;    
  }
}

# The blocking probability by simulation
p_blocking=p_blocking/T;

# The blocking probability by theory (The Erlang B formula)
p_blocking_theory=0; 
normalizing_constant=0;
for (i in 0:K){
  normalizing_constant = normalizing_constant + (lambda/mu)**i / factorial(i);
}

p_blocking_theory = (lambda/mu)**K / factorial(K);
p_blocking_theory = p_blocking_theory / normalizing_constant;


# Print the blocking probabilities (simulation and theory)
events=1:N;
plot(events, p_blocking_transient, type = "l", lty = 1)
abline(h = p_blocking_theory, col="red")

print(p_blocking)

## [1] 0.0423121

print(p_blocking_theory)

## [1] 0.04314184

This code simulates again the continuous time Markov chain associated to a M/M/K/K queue as above but in aslightl different (equivalent) manner – code written in class on March 18, 2022.

knitr::opts_chunk$set(echo = TRUE)

set.seed(17); # Set seed for reproducibility
N=1e5; # number of events to simulate;
K=10;         # number of servers
lambda=0.6*K; # arrival rate
mu=1;         # service rate

p_blocking=0; # blocking probability with T
p_blocking_transient=rep(N, 0); # transient blocking probability
p_blocking_theory=0;
T=0; # simulation time
state=0;
for (i in 1:N) {

  A=rexp(1,lambda); # Time to next job arrival
  D=A+1;
  if (state>0){
      D=rexp(1,mu*state);  # Time to next service completion
  }
  
  T=T+min(A,D); # time after i events

  # state is updated here
  if (state==K){
    p_blocking=p_blocking+min(A,D);
  }
    p_blocking_transient[i]=p_blocking/T;

  if(A>D){
    if (state>0){
      state=state-1;
    }
  } else {
    if (state<K){
      state=state+1;
    }
  }

}
p_blocking=p_blocking/T;

rho=lambda/mu; normalizing_costant=0;
for (i in 0:K){
  normalizing_costant = normalizing_costant + rho**i/factorial(i);
}
p_blocking_theory=rho**K/factorial(K);
p_blocking_theory=p_blocking_theory/normalizing_costant;

print(p_blocking)

## [1] 0.04463967

print(p_blocking_theory)

## [1] 0.04314184

events=1:N;
plot(events, p_blocking_transient, type = "l", lty = 1)
abline(h = p_blocking_theory, col="red")