Calculation and Simulation | Queuing Theory

Calculation and Simulation | Queuing Theory

Firstly, we denote some variables in the process:

Variable	Notation
$\lambda$	The arrival rate
$\mu$	The service rate
$L$	The long-term average number of customers in the system
$L_q$	The long-term average length of queue
$W$	The long-term average total waiting time of each customers
$W_q$	The long-term average waiting time of each customers for queuing
${\pi }_n$	The margin probability of state n (n represents the number of customers in the system)
$q_{ij}$	The transition rate from state $i$ to state $j$
$k$	The number of dependent queues
$\rho$	$\rho =\lambda /\mu$

I. Numerical Calculation

In order to study the queuing model, we firstly study the $\text{birth and death process}$.

In this process, we assume that ${\lambda }_i=\lambda$, ${\mu }_i=\mu$, for $\forall i$. Therefore, we can solve the transition rate:
$q_{ij}=\{\begin{array}{ll}\lambda & j=i+1\\ \mu & j=i-1\end{array}$

In $\text{birth and death process}$, according to detailed balance equation, we have that

$$\begin{gathered} {\pi }_i{q}_{ij}={\pi }_j{q}_{ji} \\ {\sum }_i{\pi }_i=1 \end{gathered}$$

Denote $\rho =\lambda /\mu$. To make the process convergent, we need to make sure that $\rho =\lambda /\mu <1$. Then, we can solve the ${\pi }_n(n\ge 0)$:

${\pi }_n=(1-\rho ){\rho }^n$

And then, we can solve the number of customers in system $L$ and lengths of queues $L_q$ (number of dependent queues $k=1$):

$$\begin{gathered} L={\sum }_{n=0}^{\infty }n{\pi }_n=\frac{\rho }{1-\rho } \\ {L}_q={\sum }_{n=k}^{\infty }(n-k){\pi }_n={\sum }_{n=1}^{\infty }(n-1){\pi }_n=\frac{{\rho }^2}{1-\rho } \end{gathered}$$

Then, according to the $\text{Little’s Law}$, we can solve the average waiting time for the whole system $W$ and only for the queuing $W_q$:

$$\begin{gathered} W=L/\lambda =\frac{\rho }{(1-\rho )\lambda }=\frac{1}{\mu -\lambda } \\ {W}_q=L_q/\lambda =\frac{{\rho }^2}{(1-\rho )\lambda }=\frac{\lambda }{\mu (\mu -\lambda )} \end{gathered}$$

Also, we can find that:

$W=W_q+\frac{1}{\mu }$

Next, we apply the result for the $\text{Birth and Death Process}$ to the $\text{Three Server Organizations}$. In all models, we apply $\text{FCFS}$ strategy:

$\text{1. Frequency-division Multiplexing (FDM):}$

In Frequency-division Multiplexing model, the system is divided into $k$ independent Birth and Death Process, and the arrival rate ${\lambda }^{\prime }=\lambda /k$, the service rate ${\mu }^{\prime }=\mu$. Denote $\rho ={\lambda }^{\prime }/{\mu }^{\prime }=\lambda /k\mu$. To make the process convergent, we need to make sure that $\rho =\lambda /k\mu <1$.

So, according to the result in the $\text{Birth and Death Process}$, we have:

$$\begin{gathered} L=k\times \frac{\rho }{1-\rho }=k\times \frac{\frac{\lambda }{k\mu }}{1-\frac{\lambda }{k\mu }}=\frac{k\lambda }{k\mu -\lambda } \\ {L}_q=k\times \frac{{\rho }^2}{1-\rho }=k\times \frac{(\frac{\lambda }{k\mu }{)}^2}{1-\frac{\lambda }{k\mu }}=\frac{{\lambda }^2}{k{\mu }^2-\lambda \mu } \end{gathered}$$
$$\begin{gathered} W=\frac{1}{\mu -\lambda /k}=\frac{k}{k\mu -\lambda } \\ {W}_q=\frac{\lambda /k}{\mu (\mu -\lambda /k)}=\frac{\lambda }{k{\mu }^2-\lambda \mu } \end{gathered}$$

$\text{2. M/M/1:}$

In M/M/1 model, the system is constructed by a single Birth and Death Process with arrival rate ${\lambda }^{\prime }=\lambda$, service rate ${\mu }^{\prime }=k\mu$. Denote $\rho ={\lambda }^{\prime }/{\mu }^{\prime }=\lambda /k\mu$. To make the process convergent, we need to make sure that $\rho =\lambda /k\mu <1$.

Therefore, according to the result in the Birth and Death Process, we have:

$$\begin{gathered} L=\frac{\rho }{1-\rho }=\frac{\frac{\lambda }{k\mu }}{1-\frac{\lambda }{k\mu }}=\frac{\lambda }{k\mu -\lambda } \\ {L}_q=\frac{{\rho }^2}{1-\rho }=\frac{(\frac{\lambda }{k\mu }{)}^2}{1-\frac{\lambda }{k\mu }}=\frac{{\lambda }^2}{k^2{\mu }^2-k\lambda \mu } \end{gathered}$$
$$\begin{gathered} W=\frac{L}{\lambda }=\frac{1}{k\mu -\lambda } \\ {W}_q=\frac{L_q}{\lambda }=\frac{\lambda }{k^2{\mu }^2-k\lambda \mu } \end{gathered}$$

$\text{3. M/M/k:}$

Compared with the forming two models, the $M/M/k$ is more complicated. It is different with the naive Birth and Death Process we solve as above. The new State Transition Diagram is as below:

In this process, we assume that ${\lambda }_i=\lambda$, ${\mu }_i=\mu$, for $\forall i$. Denote $\rho ={\lambda }^{\prime }/{\mu }^{\prime }=\lambda /k\mu$. To make the process convergent, we need to make sure that $\rho =\lambda /k\mu <1$.

Therefore, we can solve the transition rate:

$q_{ij}=\{\begin{array}{ll}\lambda & j=i+1\\ i\mu & j=i-1,i\le k\\ s\mu & j=i-1,i>k\end{array}$

In $M/M/k$ model, according to detailed balance equation, we have that

$$\begin{gathered} {\pi }_i{q}_{ij}={\pi }_j{q}_{ji} \\ {\sum }_i{\pi }_i=1 \end{gathered}$$

Denote $\rho =\lambda /k\mu$, we can solve the ${\pi }_n(n\ge 0)$:

$$\begin{gathered} {\pi }_0=[{\sum }_{n=0}^{k-1}\frac{(\lambda /\mu {)}^n}{n!}+\frac{(\lambda /\mu {)}^k}{k!}\times \frac{1}{1-\lambda /k\mu }{]}^{-1} \\ {\pi }_n=\{\begin{array}{ll}\frac{(\lambda /\mu {)}^n}{n!}\times {\pi }_0& 0\le n\le k\\ \frac{(\lambda /\mu {)}^n}{k!}\times (\frac{\lambda }{k\mu }{)}^{n-k}=\frac{(\lambda /\mu {)}^n}{k!k^{n-k}}\times {\pi }_0& n\ge k\end{array} \end{gathered}$$

Therefore, according to the result above, we can solve that:

$L_q={\sum }_{n=k}^{\infty }(n-k){\pi }_n=\frac{{\pi }_0(\lambda /\mu {)}^k\rho }{k!(1-\rho {)}^2}$
$W_q=\frac{L_q}{\lambda }=\frac{{\pi }_0\rho {\lambda }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$
$W=W_q+\frac{1}{\mu }=\frac{{\pi }_0\rho {\lambda }^{k-1}+k!(1-\rho {)}^2{\mu }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$
$L=\lambda W=\frac{{\pi }_0\rho {\lambda }^k+\lambda k!(1-\rho {)}^2{\mu }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$

After the calculation, we can also find that when $k$ gets larger, the parameters ($L,L_q,W,W_q$) of $M/M/k$ model will decrease.

I. Theoretical Calculation:

Then we can conclude the theoretical value of all of the parameters (number of customers and waiting time) of each model:

	$FDM$	$M/M/1$	$M/M/k$
$L$	$\frac{k\lambda }{k\mu -\lambda }$	$\frac{\lambda }{k\mu -\lambda }$	$\frac{{\pi }_0\rho {\lambda }^k+\lambda k!(1-\rho {)}^2{\mu }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$
$L_q$	$\frac{{\lambda }^2}{k{\mu }^2-\lambda \mu }$	$\frac{{\lambda }^2}{k^2{\mu }^2-k\lambda \mu }$	$\frac{{\pi }_0(\lambda /\mu {)}^k\rho }{k!(1-\rho {)}^2}$
$W$	$\frac{k}{k\mu -\lambda }$	$\frac{1}{k\mu -\lambda }$	$\frac{{\pi }_0\rho {\lambda }^{k-1}+k!(1-\rho {)}^2{\mu }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$
$W_q$	$\frac{\lambda }{k{\mu }^2-\lambda \mu }$	$\frac{\lambda }{k^2{\mu }^2-k\lambda \mu }$	$\frac{{\pi }_0\rho {\lambda }^{k-1}}{k!(1-\rho {)}^2{\mu }^k}$

${\pi }_0=[{\sum }_{n=0}^{k-1}\frac{(\lambda /\mu {)}^n}{n!}+\frac{(\lambda /\mu {)}^k}{k!}\times \frac{1}{1-\lambda /k\mu }{]}^{-1}$

According to the theroetical calculation, we can find that the waiting time $W,W_q$ of Frequency-division Multiplexing is always $k$ times as long as that of $M/M/1$, and the number of customers $L,L_q$ is always $k$ times as larger as that of $M/M/1$ as well when $k>1$.

Thus, the performance of $M/M/1$ is better than $FDM$.

While the parameters of $M/M/k$ model is complicated, therefore, we need to use diagrams to compare their performance.

Next, we will use diagrams of parmeters $L,L_q,W,W_q$ (with respect to the number of queues $\text{k}$) to display the theoretical calculation of each model:

from matplotlib import pyplot as plt
import numpy as np
import math

def FDM(lam, mu, K):
    l, l_q, w, w_q = 0,0,0,0
    L,L_q,W,W_q = np.array([]),np.array([]),np.array([]),np.array([])

    for i in range(len(K)):
        k = K[i]

        l = (k*lam)/(k*mu-lam)
        l_q = pow(lam,2)/(k*pow(mu,2)-lam*mu)
        w = k/(k*mu-lam)
        w_q = lam/(k*pow(mu,2)-lam*mu)

        L = np.append(L,l)
        L_q = np.append(L_q,l_q)
        W = np.append(W, w)
        W_q = np.append(W_q,w_q)

    return L,L_q,W,W_q

def M_M_1(lam, mu, K):
    l, l_q, w, w_q = 0,0,0,0
    L,L_q,W,W_q = np.array([]),np.array([]),np.array([]),np.array([])
    
    for i in range(len(K)):
        k = K[i]

        l = lam/(k*mu-lam)
        l_q = pow(lam,2)/(pow(k,2)*pow(mu,2)-k*lam*mu)
        w = 1/(k*mu-lam)
        w_q = lam/(pow(k,2)*pow(mu,2)-k*lam*mu)

        L = np.append(L,l)
        L_q = np.append(L_q,l_q)
        W = np.append(W, w)
        W_q = np.append(W_q,w_q)

    return L,L_q,W,W_q
    
def M_M_k(lam, mu, K):

    sum_1,rho,pi_0 = 0,0,0
    L,L_q,W,W_q = np.array([]),np.array([]),np.array([]),np.array([])

    for i in range(len(K)):
        k = K[i]
        # rho = lambda/(k*mu)
        rho = lam/(k*mu)

        sum_1 = sum([(pow(lam/mu,n)/math.factorial(n))  for n in range(0,k)])

        # the margin probability of state 0: pi_0
        pi_0 = 1/(sum_1+pow(lam/mu,k)/math.factorial(k)*(1/(1-lam/(k*mu))))

        # number of customers in the system
        l = ((pi_0*rho*pow(lam,k)+lam*math.factorial(k)*pow(1-rho,2)*pow(mu,k-1)) \
            /(math.factorial(k)*pow(1-rho,2)*pow(mu,k)))
        # length of queues
        l_q = ((pi_0*pow(lam/mu,k)*rho)/(math.factorial(k)*pow(1-rho,2)))
        # waiting time of the customers in the system
        w = (pi_0*rho*pow(lam,k-1)+math.factorial(k)*pow(1-rho,2)*pow(mu,k-1)) \
            /(math.factorial(k)*pow(1-rho,2)*pow(mu,k))
        # waiting time of the customers for queuing
        w_q = (pi_0*rho*pow(lam,k-1))/(math.factorial(k)*pow(1-rho,2)*pow(mu,k))

        L = np.append(L,l)
        L_q = np.append(L_q,l_q)
        W = np.append(W, w)
        W_q = np.append(W_q,w_q)

    return L,L_q,W,W_q

def plot_performance(Model, lam, mu, K):

    L = [np.array([])]*len(Model)
    L_q = [np.array([])]*len(Model)
    W = [np.array([])]*len(Model)
    W_q = [np.array([])]*len(Model)

    i = 0
    for model in Model:

        if model == "FDM":
            l, l_q, w, w_q = FDM(lam, mu, K)

        if model == "MM1":
            l, l_q, w, w_q = M_M_1(lam, mu, K)
        
        if model == "MMk":
            l, l_q, w, w_q = M_M_k(lam, mu, K)

        L[i] = np.append(L[i], l)
        L_q[i] = np.append(L_q[i], l_q)
        W[i] = np.append(W[i], w)
        W_q[i] = np.append(W_q[i], w_q)

        i += 1

    plt.figure(figsize=[16,8],dpi=80)
    plt.subplot(2,2,1)
    for i in range(len(Model)):
        plt.plot(K,L[i])
    plt.legend(Model)
    plt.title("Number of customers L (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,2)
    for i in range(len(Model)):
        plt.plot(K,L_q[i])
    plt.legend(Model)
    plt.title("Length of queues L_q (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,3)
    for i in range(len(Model)):
        plt.plot(K,W[i])
    plt.legend(Model)
    plt.title("Average Waiting Time W (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,4)
    for i in range(len(Model)):
        plt.plot(K,W_q[i])
    plt.legend(Model)
    plt.title("Average Waiting Time for Queuing W_q (lambda=%.1f,mu=%.1f)"%(lam,mu))

# model name
Model = ["FDM","MM1","MMk"]
# number of dependent queues k
K = np.arange(1,30,1)

# arrival rate lambda
lam = 0.2
# service rate mu
mu = 0.3

plot_performance(Model, lam, mu, K)

From the result above, we can find that:

According to the diagram, fixing the parameters $\lambda$ and $\mu$, we can find that as the value of k increases, the parameters (number of customers $L,L_q$ and waiting time $W,W_q$) of each model will decrease. And when $k$ is large, the paramerters will converge to a fixed value, which means that the performance of the model wll be improved in the long term.

When it comes to the comparison of each model, the conclusion has also been shown in the image:

All of the parameters of $M/M/1$ model are smaller than those of another two models, thus, it is obvious that $M/M/1$ model performs better than others. The average waiting time of customers $L$ and average waiting time $W$ of $M/M/1$ model is much smaller than the other two models. Although when $k$ is small, the length of queues $L_q$ and the average waiting time for queuing $W_q$ of $M/M/1$ model is slightly larger than $M/M/k$ model, but when $k$ gets larger, the parameters of them will converges to a similar value. Therefore, on the whole, the $M/M/1$ model performs better than others.

Also, when $k$ is small, the parameters of $M/M/k$ model are smaller than those of $FDM$ model, and when $k$ gets larger, the parameters of $M/M/k$ models are still slightly smaller than those of $FDM$. Therefore, we can conclude that the performance of $M/M/k$ is better than that of $FDM$.

To sumarize, we can make the conclusion of the performance of the three models:

$M/M/1>M/M/k>FDM$

II. Simulation

Next, we will simulate each model according to the state transition diagram. The basics of simulation is the CTMC perspective 2.

import numpy as np
from numpy import random
from matplotlib import pyplot as plt

def Transition(model, start, lam, mu, k):
    state,time = start,0
    arrival,service = 0,0

    if model == "FDM":

        arrival = random.exponential(1/(lam/k))
        service = random.exponential(1/mu)

    if model == "MM1":

        arrival = random.exponential(1/lam)
        service = random.exponential(1/(k*mu))
    
    if model == "MMk":

        arrival = random.exponential(1/lam)
        if start != 0:
            if start < k:
                service = random.exponential(1/(start*mu))
            else:
                service = random.exponential(1/(k*mu))

    if arrival < service or start == 0:
        state = start + 1
        time = arrival
    elif arrival >= service and start > 0:
        state = start - 1
        time = service

    return state, time

def CTMC_simulation(model, N, lam, mu, k):

    start,time = 0,0
    total_customer,total_time = 0,0
    num_customer, wait_time = np.array([]),np.array([])
    average_customer,average_time = 0,0
    start_state = np.array([])

    # FDM model
    if model == "FDM":
        for i in range(k):
            num_customer, wait_time = np.array([]),np.array([])
            start_state = np.array([])

            state, time = Transition(model, start, lam, mu, k)
            num_customer = np.append(num_customer,state)
            wait_time = np.append(wait_time,time)
            start_state = np.append(start_state,start)

            for i in range(N-1):
                start = state
                state, time = Transition(model, start, lam, mu, k)
                num_customer = np.append(num_customer,state)
                wait_time = np.append(wait_time,time)
                start_state = np.append(start_state,start)

                if state > start:
                    total_customer += 1

            max_n = int(np.max(num_customer)+1)
            margin_prob = np.zeros(max_n)

            for i in range(max_n):
                num_weight = np.zeros(max_n)
                for j in range(len(start_state)):
                    if i == start_state[j]:
                        num_weight[i] += wait_time[j]
                margin_prob[i] = num_weight[i]/np.sum(wait_time)

            for i in range(max_n):
                average_customer += i*margin_prob[i]
            total_time += np.sum(wait_time*start_state)

        average_time = total_time/total_customer
        L_q = average_customer-lam/mu
        W_q = average_time - 1/mu

    # M/M/1 M/M/S model
    else:
        state, time = Transition(model, start, lam, mu, k)
        num_customer = np.append(num_customer,state)
        wait_time = np.append(wait_time,time)
        start_state = np.append(start_state,start)

        for i in range(N-1):
            start = state
            state, time = Transition(model, start, lam, mu, k)
            num_customer = np.append(num_customer,state)
            wait_time = np.append(wait_time,time)
            start_state = np.append(start_state,start)
        
            if state > start:
                total_customer += 1

        max_n = int(np.max(num_customer)+1)
        margin_prob = np.zeros(max_n)

        for i in range(max_n):
            num_weight = np.zeros(max_n)
            for j in range(len(start_state)):
                if i == start_state[j]:
                    num_weight[i] += wait_time[j]
            margin_prob[i] = num_weight[i]/np.sum(wait_time)

        for i in range(max_n):
            average_customer += i*margin_prob[i]

        average_time = np.sum(wait_time*start_state)/total_customer

        if model == "MM1":
            L_q = average_customer - lam/(k*mu)
            W_q = average_time - 1/(k*mu)
        elif model == "MMk":
            L_q = average_customer - lam/mu
            W_q = average_time - 1/mu

    L = average_customer
    W = average_time
    
    return L, L_q, W, W_q

def plot_CTMC(Model, N, lam, mu, K, num_exper):

    average_L = np.zeros([len(Model),len(K)])
    average_W = np.zeros([len(Model),len(K)])
    average_L_q = np.zeros([len(Model),len(K)])
    average_W_q = np.zeros([len(Model),len(K)])

    for n in range(num_exper):

        l,l_q,w,w_q = 0,0,0,0

        # number of customers in the system
        L = np.zeros([len(Model),len(K)])
        # length of queues
        L_q = np.zeros([len(Model),len(K)])
        # average waiting time of customers in the system
        W = np.zeros([len(Model),len(K)])
        # average waiting time of customers for queuing
        W_q = np.zeros([len(Model),len(K)])

        j = 0
        for model in Model:

            for i in range(len(K)):
                k = K[i]

                l,l_q,w,w_q = CTMC_simulation(model, N, lam, mu, k)
                L[j][i] += l
                L_q[j][i] += l_q
                W[j][i] += w
                W_q[j][i] += w_q

                # # Theorical value of W
                # W[j] = L[j]/lam

            j += 1

    for i in range(len(Model)):
        average_L[i] = L[i]/num_exper
        average_W[i] = W[i]/num_exper
        average_L_q[i] = L_q[i]/num_exper
        average_W_q[i] = W_q[i]/num_exper

    plt.figure(figsize=[16,8],dpi=80)
    plt.subplot(2,2,1)
    for i in range(len(Model)):
        plt.plot(K,average_L[i])
    plt.legend(Model)
    plt.title("Number of customers L (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,2)
    for i in range(len(Model)):
        plt.plot(K,average_L_q[i])
    plt.legend(Model)
    plt.title("Length of queues L_q (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,3)
    for i in range(len(Model)):
        plt.plot(K,average_W[i])
    plt.legend(Model)
    plt.title("Average Waiting Time W (lambda=%.1f,mu=%.1f)"%(lam,mu))

    plt.subplot(2,2,4)
    for i in range(len(Model)):
        plt.plot(K,average_W_q[i])
    plt.legend(Model)
    plt.title("Average Waiting Time for Queuing W_q (lambda=%.1f,mu=%.1f)"%(lam,mu))

    return average_L, average_W

# number of experiment
num_exper = 100
# total times of CTMC experiment
N = 600
# model name
Model = ["FDM","MM1","MMk"]

# arrival rate
lam = 0.2
# service rate
mu = 0.3
# number of queues
K = np.arange(1,31,1)

L,W = plot_CTMC(Model, N, lam, mu, K, num_exper)

From the simulation above, we can find that the results of simulation is similar to the theoretical calculation. And according to the results shown in the diagrams, the performance of models is:

$M/M/1>M/M/k>FDM$

III Pros and Cons of Each System

1. FDM

Strengths:

In FDM system, the entire process is divided into $k$ independent parts, which makes the FDM system is most robust. This feature allows any part of the model to be abnormal without affecting other processes, which is beneficial in reality.

Weaknesses:

The parameters of FDM ($L,L_q,W,W_q$) of FDM are always the largest among the three models, and thus, FDM model doesn’t perform well.

2. M/M/1

Strengths:

$M/M/1$ system has the best performance among the three systems. The parameters of it are relatively smaller. Therefore, the system can serve the most customers within a certain time and has the shortest queues.

Weaknesses:

The system is centralized to one process and less robust, which will lead to the most prone to failure in reality.

3. M/M/k

Strengths:

$M/M/k$ system has relatively better performances. When $k$ is small, the length of queues $L_q$ and the average waiting time for queuing of $M/M/1$ is the smallest. And it divides its service process into $k$ parts, which makes it more robust.

Weaknesses:

Compared with the $M/M/1$ system, $M/M/k$ system behaves worse in performance. And $M/M/k$ system is not as robust as FDM system.

IV. Explanation

In the process of simulation, it can be found that due to the dispersion of service resources, although the FDM and $M/M/1$ model have more robustness, they will also lead to the waste of certain service resources. Different from the MM1 model, these two models split the service resources into k parts, which leads to the fact that when a customer enters, the $M/M/1$ model will complete the service at the fastest speed (rate = $k\mu$), while the other two The speed (rate = $\mu$) of the model is slower in comparison.

Also, we can explain the phenomenon that the performance of $M/M/1$ is better than $M/M/k$ according to the $\text{State Transition Diagram}$:
We can find that when the number of customer in the system $n<k$, the service rate of $M/M/k$ is $n\mu$, while the service rate of $M/M/1$ is $k\mu$. And when $n\ge k$, the service rate of them are the same: $k\mu$. Also, the arrival rate of each model in all states are the same $\lambda$. Therefore, the performance of $M/M/1$ is better than $M/M/k$.

Therefore, it is intuitive that the performance of $\text{M/M/1}$ model is the best among the three models. $\text{M/M/k}$ is the second and $\text{FDM}$ is the worst.