Queueing_model

E. Quinn 3/3/2022

M/M/m queueing model for ICU capacity

Import standard python datascience packages

import sys
import math
import re
import copy as cp
import numpy as np
import scipy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pickle
%matplotlib inline

from datetime import datetime, timedelta, date
from datascience import *
import uuid
import random

from scipy.stats import nbinom
import numpy.random as npr

Overview

A major goal of Covid-19 mitigation strategies is to prevent the healthcare system from being overwhelmed with seriously ill Covid-19 patients.

According to the Johns Hopkins Coronavirus database, Rhode Island has 214 ICU beds.

Published numbers suggest the average ICU stay is abut 15.5 days: ICU LENGTH OF STAY IN CRITICALLY ILL PATIENTS WITH COVID-19: DOES RACE MATTER? CHEST Journal Volume 160, ISSUE 4, SUPPLEMENT, A1164, October 1, 202102517-4/fulltext) .

This notebook does a back-of-the-envelope queueing model (M/M/214 in Kendall's notation) to get an idea of the risk of exceeding capacity given the characteristics of the current pandemic in Rhode Island.

Formulation of the M/M/m model

The important input parameters are:

• $\lambda$ The arrival rate of seriously ill patients
• $\mu$ The departure rate of seriously ill patients
• $m$ The number of ICU beds available (214 in this case)

The important derived quantities are:

The traffic intensity:

$$\rho = \frac{\lambda}{m\mu}$$

The probability of the state with zero patients $P_0$:

$$P_0 = \left[\sum_{i=1}^{m-1}\frac{(m\rho)^i}{i!} + \frac{(m\rho)^m}{m!(1-\rho)}\right]^{-1} $$

The probability of the state with $i$ patients in the system:

$$P_i = \left\{\begin{array}{lcl}P_0\frac{(m\rho)^i}{i!}&for&i\leq m\\P_0\frac{m^m\rho^i}{m!}&for&i>m\end{array}\right.$$

The probability $P_Q$ that an arriving patient has to wait in the queue is (this is called the Erlang C formula):

$$P_Q = \frac{P_0(m\rho)^m}{m!(1-\rho)}$$

Compute and normalize queue length probabilities

def compute_P(m,mu,L,maxbeds):         #maxbeds,avg stay,arrival rate, highest number of beds to compute
    P = np.zeros(maxbeds)
    P[0] = 1.0
    for i in np.arange(1,maxbeds):
        if (i <= m):
            P[i] = P[i-1]*L/(mu*i)
        else:
            P[i] = P[i-1]*L/(mu*m)
    Ptot = 0.0
    for i in np.arange(maxbeds-1,-1,-1):
        Ptot += P[i]
    for i in np.arange(maxbeds-1,-1,-1):
        P[i] = P[i]/Ptot
        
    return(P)

Compute M/M/214 queueing system model

def M_M_m2(arrival=10.,severe_stay=17,fname='fig1'):
    severe_case_length_of_stay = severe_stay

    m = 214
    maxbeds = 300
    mu = 1/severe_case_length_of_stay

    Lambda = arrival
    rho = Lambda/(m*mu)
    print("Lambda: ",Lambda,' mu: ',mu,' rho: ',rho)
    
    title_string ='Arrival rate: ' + str(Lambda) + '    Average stay: ' + str(severe_stay) + ' days'
    
    P = compute_P(m,mu,Lambda,maxbeds)
    
    avgQ = 0.0
    for i in np.arange(len(P)):
        avgQ+= i*P[i]
    print('avgQ: ',avgQ)
    
    PQ = 0.0
    for i in np.arange(len(P)):
        if (i > m):
            PQ+=P[i]
    print('PQ: ',PQ )

    ICU = np.arange(len(P))
    
    fig, ax = plt.subplots()

    ax.set_title(title_string, fontsize = 18)
    ax.set_ylabel('Probability', fontsize = 18)
    ax.set_xlabel('Number of ICU Beds', fontsize = 18)

    plt.rcdefaults()
    plt.rcParams['figure.figsize'] = [10, 6]

    ax.bar(ICU,P)
    ax.set_ylabel('Probability',fontsize=24)
    ax.set_xlabel('Number of ICU Beds', fontsize = 18)
    ax.set_title(title_string,fontsize=18)
    plt.axvline(x=214)

    fig.tight_layout()
    pngfile = '../ICU_model2_MMn_' + fname + '.png'
    
    plt.savefig(pngfile)

M/M/214 queueing model: arrival rate=1

• Average length of stay: 15.5 days

M_M_m2(arrival=1.0,severe_stay=15.5,fname='fig1')

Lambda:  1.0  mu:  0.06451612903225806  rho:  0.07242990654205607
avgQ:  15.499999999999998
PQ:  3.61406382544511e-161

M/M/214 queueing model: arrival rate=6

• Average length of stay: 15.5 days

M_M_m2(arrival=6.0,severe_stay=15.5,fname='fig2')

Lambda:  6.0  mu:  0.06451612903225806  rho:  0.43457943925233644
avgQ:  93.00000000000007
PQ:  2.6161789310739812e-27

M/M/214 queueing model: arrival rate=10

• Average length of stay: 15.5 days

M_M_m2(arrival=10.0,severe_stay=15.5,fname='fig3')

Lambda:  10.0  mu:  0.06451612903225806  rho:  0.7242990654205608
avgQ:  155.00001149173968
PQ:  3.1682833698711115e-06

M/M/214 queueing model: arrival rate=12

• Average length of stay: 15.5 days
• The probability that there are no beds available when a patient arrives is 2.4%

M_M_m2(arrival=12.0,severe_stay=15.5,fname='fig4')

Lambda:  12.0  mu:  0.06451612903225806  rho:  0.8691588785046729
avgQ:  186.1848240008622
PQ:  0.02418496077745914

M/M/214 queueing model: arrival rate=13

• Average length of stay: 15.5 days
• There is a 27% chance that there will be no beds available when a patient arrives.

M_M_m2(arrival=13.0,severe_stay=15.5,fname='fig5')

Lambda:  13.0  mu:  0.06451612903225806  rho:  0.9415887850467289
avgQ:  205.91175373489764
PQ:  0.2668713007771805

M/M/214 queueing model: arrival rate=13.5

• Average length of stay: 15.5 days
• The system is just on the verge of being overwhelmed
• There is a 60% chance that there will be no beds available when a patient arrives

M_M_m2(arrival=13.5,severe_stay=15.5,fname='fig6')

Lambda:  13.5  mu:  0.06451612903225806  rho:  0.977803738317757
avgQ:  226.90833944859776
PQ:  0.6003283786196211

M/M/214 queueing model: arrival rate=14

• Average length of stay: 15.5 days
• The system is overwhelmed
• There is no equilibrium
• The queue length just grows without bound

M_M_m2(arrival=14.0,severe_stay=15.5,fname='fig7')

Lambda:  14.0  mu:  0.06451612903225806  rho:  1.014018691588785
avgQ:  259.6106837580506
PQ:  0.9091842777890705