$M/M/2/4$ simulation of the probability that the queue gets full during first $10$ time units.

Question

Let $X(t)$ denote the total number of customers at time $t \geq 0$ in an $M/M/2/4$ queuing system in steady-state (/started according to its stationary distribution) with Poisson arrival process with rate $\lambda = 1$ and with exponentially distributed service times with mean $1$.

Describe a numerical procedure (/simulation procedure) that gives an approximate numerical value for the probability $$P\left( \max_{0\leq t \leq 10} X(t) \right) = 4$$ that the queuing system gets full during the first $10$ time units.

Expected answer: $0.46679$

My attempt at a solution in R code:

rm(list=ls())
N=10^6

lambda = 1
mu = 1
count = 0
change = 0


for (i in 1:N) {

   time = 0

   succ = 0

   x0 = runif(1)            

   if (x0<(8/23)) {                   % Here i determine where to start
   q = 0                              % This part should be correct i think
   }

   else if ((8/23)<=x0 && x0<(16/23)) {
   q = 1
   }

   else if ((16/23)<=x0 && x0<(20/23)) {
   q = 2
   }

   else if ((20/23)<=x0 && x0<(22/23)) {
   q = 3
   }

   else if (x0>=(22/23)) {
   q = 4
   succ = 1
   }

   while ((time<=10) && (succ==0)) {

      t1 = rexp(1,lambda)                % arrival times
      s1 = rexp(1,mu)                    % service times
      s2 = rexp(1,mu)

      if (q==0) {                        % here we can only have arrivals

         time = time + t1
         change = +1
       }

       if (q==1) {

           if (t1<s1) {                  % when q=n i check whether arrival
                                         % times or service times are 
           time = time + t1              % greatest and from there i 
           change = +1                   % determine which time i add to                               
        }                                % time and which change i make to
                                         % to the queue.
        else if (t1>s1) {

           time = time + s1
           change = -1
      }
    }

      if (q==2) {

          if ((t1<s1) && (t1<s2)) {

             time = time + t1
             change = +1
          }

          else if (t1>s1) {

             time = time + s1
             change = -1
          }

          else if (t1>s2) {

             time = time + s2
             change = -1
          }
      }

        if (q==3) {

            if ((t1<s1) && (t1<s2)) {

               time = time + t1
               change = +1
            }

            else if (t1>s1) {

               time = time + s1
               change = -1
            }

             else if (t1>s2) {

               time = time + s2
               change = -1
             }
        }

        if (q==4) {
          succ = succ + 1                  % Here i register the successes
        }

q = q + change                             % Here i change the queue
}

count = count + succ                       % Here i add up the successes

}
p = count/N                              % and this frequency should be the 
                                         % sought after probability

cat("P{max 0<=t<=10 X(t)=4} = ", p ,"\n")

Answer 1

I think there's a bug in your time comparisons for cases q==2 and q==3.

Try the changes on the lines marked:

      if (q==2) {

          if ((t1<s1) && (t1<s2)) {

             time = time + t1
             change = +1
          }

          else if (s2>s1) {              %%% This line              
             time = time + s1
             change = -1
          }

          else if (s1>s2) {              %%% This line

             time = time + s2
             change = -1
          }
      }

Same change needed for q == 3 case.

The first IF determines whether $t1$ is the smallest so the next two IFs must determine which of $s1$ and $s2$ is the smallest.