# Open Jackson Network Lab

In the lab you will learn

* How to represent arrival rates in a network as a system of linear equations in matrix form
* How to use `scipy.linalg` functions to solve for the arrival rates
* How to analyse the network

## Imports


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

## The example queueing network

**Question**: What is the expected number of customers in the entire queuing system?




## Step 1: Solve for $\lambda_i$

Given:

\begin{equation}
        \lambda_j = a_j + \sum_{i=1}^m \lambda_i p_{ij}
    \end{equation}
    
---

Plug the $a_j$ and $p_{ij}$ figures into the equation for $\lambda_j$ 

$$\lambda_1 = 1 + 0.1\lambda_2 + 0.4\lambda_3 $$
$$\lambda_2 = 4 + 0.6\lambda_1 + 0.4\lambda_3 $$
$$\lambda_3 = 3 + 0.3\lambda_1 + 0.3\lambda_2 $$

---
(optional) rearrange the equations with $a_j$ on the RHS

$$\lambda_1 - 0.1\lambda_2 - 0.4\lambda_3 = 1$$
$$-0.6\lambda_1 + \lambda_2 - 0.4\lambda_3 = 4$$
$$-0.3\lambda_1 - 0.3\lambda_2 + \lambda_3 = 3$$

---

Express the system of equations describing arrival rates in matrix form.

$$ \begin{bmatrix} 1 & -0.1 & -0.4\\ -0.6 & 1 & -0.4 \\ -0.3 & -0.3 & 1 \end{bmatrix} \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{bmatrix}= \begin{bmatrix} 1\\ 4\\ 3 \end{bmatrix}$$


```python
# create the matricies as numpy arrays
a = np.array([[1, -0.1, -0.4], [-0.6, 1, -0.4], [-0.3, -0.3, 1]])
b = np.array([1, 4, 3])
```


```python
a
```




    array([[ 1. , -0.1, -0.4],
           [-0.6,  1. , -0.4],
           [-0.3, -0.3,  1. ]])




```python
b
```




    array([1, 4, 3])



## Solve for the unknowns

To solve for the unknowns we need to take the dot product of the inverse of a and b

$$\begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{bmatrix} = \begin{bmatrix} 1 & -0.1 & -0.4\\ -0.6 & 1 & -0.4 \\ -0.3 & -0.3 & 1 \end{bmatrix}^{-1} \begin{bmatrix} 1\\ 4\\ 3 \end{bmatrix}$$


```python
np.dot(linalg.inv(a), b)
```




    array([ 5. , 10. ,  7.5])



An alternative way to solve for the unknowns is to use linalg.solve


```python
arrival_rate = linalg.solve(a, b)
arrival_rate
```




    array([ 5. , 10. ,  7.5])



## Step 2: Calculate the traffic intensity

$\rho_i = \dfrac{\lambda_i} {s_i \mu_i}$

Remember that you are treating each queuing system as an independent M/M/s

* Facility 1 is a M/M/1
* Facility 2 is a M/M/2
* Facility 3 is a M/M/1

#### Example solution


```python
servers = np.array([1, 2, 1])
service_rate = np.array([10, 10, 10])

rho = arrival_rate / (servers * service_rate)
rho
```




    array([0.5 , 0.5 , 0.75])



## Step 3: Calculate the number of people in the queue

For an M/M/1

$L_s = \dfrac{\rho}{1 - \rho}$

and for an M/M/s

$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}$

$L_q = \dfrac{p_o(\lambda/\mu)^s\rho}{s!(1-\rho)^2}$



$L_s = L_q + \dfrac{\lambda}{\mu}$



```python
class MMSQueue:
    '''
    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)
        
    @property
    def total_in_system(self):
        return self.metrics[r'$L_s$']
        
    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
queues = [MMSQueue(l, mu, s) for l, mu, s in zip(arrival_rate, service_rate, 
                                                 servers)]

queues[0].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.5</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>0.5</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>1.0</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.2</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.1</td>
    </tr>
  </tbody>
</table>
</div>




```python
queues[1].summary_frame().round(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>performance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\rho$</th>
      <td>0.50</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>0.33</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>1.33</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.13</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.03</td>
    </tr>
  </tbody>
</table>
</div>




```python
queues[2].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.75</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>2.25</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>3.00</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.40</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.30</td>
    </tr>
  </tbody>
</table>
</div>



#### Example solution


```python
# total in system
total_in_system = sum([q.total_in_system for q in queues])
print(f'Expected # of customers in queuing network = {total_in_system:.2f}')
```

    Expected # of customers in queuing network = 5.33


## Convenience class for Network problem


```python
class JacksonNetwork:
    '''
    Open Jackson Queuing Network
    
    A simple implementation where a system of linear equations is specified in
    matrix form by a user as opposed to manually adding nodes.
    
    Create and conveniently holds a collection of MMSQueues.
    '''
    def __init__(self, a, b, service_times, servers):
        '''
        Init the network
        
        Params:
        ------
        a: np.ndarray
            routing of entities through network 
            Note this rearranged for calculation.
            
        b: np.ndarray
            arrival rates to nodes from outside the network
            
        service_times: np.ndarray
            Service rate at each node
            
        servers: np.ndarray
            No. of servers at each node
                
        '''
        self.a = a
        self.b = b
        self.service_times = service_times
        self.servers = servers
        
        # number of queues
        self.n_nodes = len(service_times)
        
        # solve for arrival rates
        self.arrival_rate = linalg.solve(a, b)
        
        # create MMS queues
        self.queues = self._create_mms_queues()
           
    @property
    def total_in_system(self):
        '''
        The expected number of entities in the jackson network
        Sum of expected number of entities in each MMSQueue
        
        Returns:
        -------
        float
        '''
        return sum([q.total_in_system for q in queues])
    
    def summary_frame(self):
        '''
        Tabular summary the performance measures in each queuing node
        
        For each node returns Rho, L_q, L_s, W_q, W_s, 
        
        Returns:
        --------
        pd.DataFrame
        '''
        summary = pd.DataFrame()
        
        for node_id, q in zip(range(len(self.queues)), self.queues):
            node_summary = q.summary_frame()
            node_summary.columns = ['node_' + str(node_id)]
            summary = pd.concat([summary, node_summary], axis=1)
        
        return summary
        
    
    def _create_mms_queues(self):
        '''
        Create and parameterise the required number of MMSQueue objects
        
        Returns:
        --------
        List
        
        '''
        return [MMSQueue(l, mu, s) for l, mu, s in zip(self.arrival_rate, 
                                                       self.service_times, 
                                                       self.servers)]
    
```


```python
model = JacksonNetwork(a, b, service_rate, servers)
model.summary_frame().round(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>node_0</th>
      <th>node_1</th>
      <th>node_2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>$\rho$</th>
      <td>0.5</td>
      <td>0.50</td>
      <td>0.75</td>
    </tr>
    <tr>
      <th>$L_q$</th>
      <td>0.5</td>
      <td>0.33</td>
      <td>2.25</td>
    </tr>
    <tr>
      <th>$L_s$</th>
      <td>1.0</td>
      <td>1.33</td>
      <td>3.00</td>
    </tr>
    <tr>
      <th>$W_s$</th>
      <td>0.2</td>
      <td>0.13</td>
      <td>0.40</td>
    </tr>
    <tr>
      <th>$W_q$</th>
      <td>0.1</td>
      <td>0.03</td>
      <td>0.30</td>
    </tr>
  </tbody>
</table>
</div>




```python
print(f'Expected # of customers in queuing network = {model.total_in_system:.2f}')
```

    Expected # of customers in queuing network = 5.33


696fc2952df86c00a92fd4b934654a990b0ab2d8
ipynb
Jupyter Notebook
qt_lab_jackson_networks.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
qt_lab_jackson_networks.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
qt_lab_jackson_networks.ipynb
TomMonks/jackson-network
bbb6fce334803b771fa57cc055fe06c5e93e3cdb
MIT
Qwen/Qwen-72B
1. YES
2. YES