# M/M/s/$\infty$ queues

**In this lab you will learn:**

* How to convert M/M/S queuing equations into python functions
* How to analyse M/M/s queuing systems to inform health system improvement and design.


```python
from scipy import math
import numpy as np
import pandas as pd
```

An **M/M/s** system is a queuing process having Poisson arrival pattern, $s$ servers with $s$ i.i.d expeonential service times.  Service times do not depend on the state of the system.  The system (i.e. queue + service) has infinite capacity and a FIFO queue discipline.

#### Traffic Intensity
\begin{equation*}
\rho = \frac{\lambda}{s\mu}
\label{eq:rho} \tag{1}
\end{equation*}

#### Inference about the number of patients in the system

\begin{equation*}
P_0 = \left[ \sum_{n=0}^{s-1} \frac{\left(\lambda/ \mu \right)^n}{n!} + \frac{\left( \lambda / \mu \right)^s}{s!\left(1-\rho\right)}  \right]^{-1}
\label{eq:p0} \tag{2}
\end{equation*}


\begin{equation}
  P_n = \left\{
    \begin{array}{l}
      \dfrac{\left( \lambda / \mu \right)^n}{n!}p_0,   \>\>\>\>  n \leq s\\ \\
      \dfrac{\left( \lambda / \mu \right)^n}{s!s^{n-s}}p_0,  \>\>\>\> n > s
    \end{array}
  \right.
\label{eq:pn} \tag{3} 
\end{equation}

#### Expected number of customers in the queue for service

\begin{equation}
L_q = \frac{p_0\left(\lambda / \mu \right)^s \rho}{s!\left(1 - \rho\right)}
\tag{4}
\end{equation}

#### Little's Formula

\begin{equation}
L_s = \lambda W_s  \\ L_q = \lambda W_q
\tag{5a, 5b}
\end{equation}

\begin{equation*}
    W_s = W_q + \dfrac{1}{\mu} 
    \tag{6}
\end{equation*}

\begin{equation*}
    L_s = L_q + \dfrac{\lambda}{\mu}
    \tag{7}
\end{equation*}

# Hospital Pharmacy example

During the afternoon, a pharmacy based in a large hospital has 3 trained pharmacists on duty to check and fullfill patient prescriptions for drugs to take home with them at discharge. They are are able to handle 15 transactions per hour. The service times are exponentially distributed. During this busy period, prescriptions arrive at the pharmacy according to a Possion process, at a mean rate of 40 per hour.  

**Questions**

1. What is the probability that there are more than 3 prescriptions in the pharmacy at any one time
2. Calculate the expected number of drug prescriptions waiting to be fulfilled
3. Calcluate the expected number of drug prescriptions in the system
4. Calculate the expected prescription turnaround time

### Example Solution:

This is a M/M/3 system with $\lambda=40$ and $\mu = 15$

#### Is the system in control?

Let's first check that steady state conditions hold by calculating the traffic intensity $\rho$.

\begin{equation}
\rho = \frac{\lambda}{s\mu}
\label{eq:rho} \tag{1}
\end{equation}

Steady state conditions hold if $\rho < 1$


```python
def traffic_intensity(_lambda, mu, s):
    '''
    calculate the traffic intensity (server utilisation)
    of an M/M/s queue
    '''
    return _lambda / (s * mu)
```


```python
#calculate traffic intensity
LAMBDA = 40
MU = 15
S = 3

rho = traffic_intensity(LAMBDA, S, MU)
rho
```




    0.8888888888888888



**Conclusion**: $\rho < 1$ steady state conditions will hold.

### 1. Calculate the probability that there are 3 drug orders in the pharmacy at any one time

Steady state probabilities are given by

\begin{equation*}
P_0 = \left[ \sum_{n=0}^{s-1} \frac{\left(\lambda/ \mu \right)^n}{n!} + \frac{\left( \lambda / \mu \right)^s}{s!\left(1-\rho\right)}  \right]^{-1}
\label{eq:p0} \tag{2}
\end{equation*}


\begin{equation}
  P_n = \left\{
    \begin{array}{l}
      \dfrac{\left( \lambda / \mu \right)^n}{n!}p_0,   \>\>\>\>  n \leq s\\ \\
      \dfrac{\left( \lambda / \mu \right)^n}{s!s^{n-s}}p_0,  \>\>\>\> n > s
    \end{array}
  \right.
\label{eq:pn} \tag{3} 
\end{equation}


```python
def prob_system_empty(_lambda, mu, s):
    '''
    The probability that a M/M/s/infinity queue is empty
    '''
    p0 = 0.0
    rho = traffic_intensity(_lambda, mu, s)
    
    for n in range(s):
        p0 += ((_lambda / mu) ** n) / math.factorial(n)

    p0 += ((_lambda / mu) ** s) / (math.factorial(s) * (1 - rho))
    return p0**-1
```


```python
p0 = prob_system_empty(LAMBDA, MU, S)
print(f'p0 = {p0:.2f}')
```

    p0 = 0.03



```python
def prob_n_in_system(n, _lambda, mu, s, return_all_solutions=True):
    '''
    Calculate the probability that n customers
    in the system (queuing + service)
    
    Parameters:
    --------
    n: int,
        Number of customers in the system
    
    _lambda: float
        Mean arrival rate to system
        
    mu: float
        Mean service rate
        
    s: int
        number of servers
        
    return_all_solutions: bool, optional (default=True)
        Returns all solutions for 0,1 ... n
        
    Returns:
    ------
        np.ndarray of solutions
    
    '''
    p0 = prob_system_empty(_lambda, mu, s)
    probs = [p0]
    
    #for n <= s
    for i in range(1, min(s+1, n+1)):
        pn = (((_lambda / mu)**i) / math.factorial(i)) * p0
        probs.append(pn)
        
    #for n > s
    for i in range(s+1, n+1):
        pn = (((_lambda / mu)**i) / (math.factorial(s) * (s**(n-s)))) * p0
        probs.append(pn)
    
    if return_all_solutions:
        return np.array(probs)
    else:
        return probs[:-1]
```


```python
prob = prob_n_in_system(3, LAMBDA, MU, S)

#returns: [p0, p1, p2, p3] => probabilities of 3 or less drug orders
prob.sum()
```




    0.29110418830045015




```python
#prob.sum() => p(X <=3)
more_than_three = 1 - prob.sum()
print(f'P(X > 3) = {more_than_three:.2f}')
```

    P(X > 3) = 0.71


### 2. Expected number of drug prescriptions waiting to be fullfilled


$L_q$ = Expected number of customers in the queue for service

\begin{equation}
L_q = \frac{p_0\left(\lambda / \mu \right)^s \rho}{s!\left(1 - \rho\right)^2}
\tag{4}
\end{equation}


```python
def mean_queue_length(_lambda, mu, s):
    '''
    Mean length of queue Lq
    '''
    p0 = prob_system_empty(_lambda, mu, s)
    rho = traffic_intensity(_lambda, mu, s)
    
    lq = (p0 * ((_lambda / mu)**s) * rho) / (math.factorial(s) * (1 - rho)**2)
    return lq
```


```python
lq = mean_queue_length(LAMBDA, MU, S)
```


```python
print(f'lq = {lq:.2f}')
```

    lq = 6.38


### 3. Expected number of drug prescriptions in the system

$L_s$ = Expected number of customers in the queue

We have already calculated $L_q$ therefore we will use

\begin{equation}
L_s = L_q + \frac{\lambda}{\mu}
\tag{5}
\end{equation}


```python
ls = lq + (LAMBDA / MU)
```


```python
print(f'Ls = {ls:.2f}')
```

    Ls = 9.05


### 4. Expected prescription turnaround time

Using:

\begin{equation}
L_s = \lambda W_s
\tag{5}
\end{equation}


\begin{equation}
\frac{L_s}{\lambda} = W_s
\end{equation}


```python
ws = ls / LAMBDA
```


```python
print(f'Ws = {ws:.2f}')
```

    Ws = 0.23


## MMsQueue Class

A somewhat cleaner way of analytic modelling of queues is to implement a class.  An example implementation is below.


```python
class MMsQueue(object):
    '''
    M/M/S/inf/inf/FIFO system
    '''
    def __init__(self, _lambda, mu, s):
        '''
        Constructor
        
        Parameters:
        -------
        _lambda: float
            The arrival rate of customers to the facility
            
        mu: float
            The service rate of the facility
            
        s: int
            The number of servers
        '''
        self._lambda = _lambda
        self.mu = mu
        self.s = s
        self.rho = self._get_traffic_intensity()
        
        #create a dict of performance metrics
        #solve for L_q then use little's law to calculate remaining KPIs
        self.metrics = {}
        self.metrics[r'$\rho$'] = self.rho
        self.metrics[r'$L_q$'] = self._get_mean_queue_length()
        self.metrics[r'$L_s$'] = self.metrics[r'$L_q$'] + (_lambda / mu)
        self.metrics[r'$W_s$'] = self.metrics[r'$L_s$'] / _lambda
        self.metrics[r'$W_q$'] = self.metrics[r'$W_s$'] - (1 / mu)
        
    def _get_traffic_intensity(self):
        '''
        calculate the traffic intensity (server utilisation)
        of an M/M/s queue
        '''
        return self._lambda / (self.s * self.mu)  
    
    def _get_mean_queue_length(self):
        '''
        Mean length of queue Lq
        '''
        p0 = self.prob_system_empty()
       
        lq = (p0 * ((self._lambda / self.mu)**self.s) * 
              self.rho) / (math.factorial(self.s) * (1 - self.rho)**2)
        return lq
        
    def prob_system_empty(self):
        '''
        The probability that a M/M/s/infinity queue is empty
        '''
        p0 = 0.0

        for n in range(self.s):
            p0 += ((self._lambda / self.mu) ** n) / math.factorial(n)

        p0 += ((self._lambda / self.mu) ** self.s) / (math.factorial(self.s) 
                                                      * (1 - self.rho))
        return p0**-1
    
    def prob_n_in_system(self, n, return_all_solutions=True, as_frame=True):
        '''
        Calculate the probability that n customers
        in the system (queuing + service)

        Parameters:
        --------
        n: int,
            Number of customers in the system

        return_all_solutions: bool, optional (default=True)
            Returns all solutions for 0,1 ... n
            
        as_frame: bool, optional (default=True)
            If True, returns all solutions in a pd.DataFrame
            else returns all solutions as np.ndarray
            has no effect is return_all_solutions == False

        Returns:
        ------
            np.ndarray of solutions

        '''
        p0 = self.prob_system_empty()
        probs = [p0]

        #for n <= s
        for i in range(1, min(self.s+1, n+1)):
            pn = (((self._lambda / self.mu)**i) / math.factorial(i)) * p0
            probs.append(pn)

        #for n > s
        for i in range(self.s+1, n+1):
            pn = (((self._lambda / self.mu)**i) / (math.factorial(self.s) 
                                                   * (self.s**(n-self.s)))) * p0
            probs.append(pn)

        if return_all_solutions:
            results = np.array(probs)
            if as_frame:
                return pd.DataFrame(results, columns=['P(X=n)'])
            else:
                return results
        else:
            return probs[:-1]
        
    def summary_frame(self):
        '''
        Return performance metrics
        
        Returns:
        ---------
        pd.DataFrame
        '''
        df = pd.Series(self.metrics).to_frame()
        df.columns = ['performance']
        return df
```


```python
model = MMsQueue(LAMBDA, MU, S)
model.summary_frame()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>performance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\rho$</th>
      <td>0.888889</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>6.380062</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>9.046729</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.226168</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.159502</td>
    </tr>
  </tbody>
</table>
</div>




```python
model.prob_n_in_system(5)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>P(X=n)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.028037</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.074766</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.099688</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.088612</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.026255</td>
    </tr>
    <tr>
      <th>5</th>
      <td>0.070014</td>
    </tr>
  </tbody>
</table>
</div>




```python
#county hospital example
model = MMsQueue(2, 3, 2)
model.summary_frame()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>performance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\rho$</th>
      <td>0.333333</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>0.083333</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>0.750000</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.375000</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.041667</td>
    </tr>
  </tbody>
</table>
</div>




```python
model.prob_n_in_system(2)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>P(X=n)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.500000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.333333</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.111111</td>
    </tr>
  </tbody>
</table>
</div>



4bbfd3b8a5027b00b5a5ad36f92a349e14da1e53
ipynb
Jupyter Notebook
qt_lab_mms.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
qt_lab_mms.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
qt_lab_mms.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
Qwen/Qwen-72B
1. YES
2. YES