Over the last 30 or so years, I've built numerous healthcare discrete event simulation (DES) models. I built models as a student, as an analyst for several healthcare systems and consulting firms, and as a teacher and researcher at Oakland University. I've built models of inpatient units, pharmacies, pneumatic tube systems, outpatient clinics, surgical and recovery suites, emergency rooms and a whole bunch more. I've built models in several different commercial discrete event simulation languages/packages such as SIMAN, Arena, MedModel, and Simio. While each of these are very fine products, one thing that's always bothered me was how difficult to it was to share, review, and improve upon the models built by the healthcare simulation modeling community. These packages are really expensive often have pretty restrictive licensing. So, academics and practitioners publish articles describing how a simulation model of some healthcare related system (e.g. inpatient beds, pharmacy, lab, surgical suites...) was built and used in some project. Great. Inspiring. Not super helpful for the person looking to build a similar model, or better yet, build upon that same model. Models were rarely if ever made available and even if they were, might be written in a package that you don't have, don't know, or can't afford.

As the years went on, I found myself more and more attracted to the world of open source software. I fell in love with Linux and R and Python. However, open source discrete event simulation options were not plentiful. I did have success building a pretty big pneumatic tube simulation model in Java using the Simkit library developed by Arnie Buss and others at the Naval Postgraduate School. While I really like Simkit, I'm not a huge fan of programming in Java. Most of my research work and teaching is focused on R and Python and ideally would love to have a Python based DES framework. So, I was pretty happy when I stumbled onto the SimPy project and followed it for a bit as it moved through its early phases of growth and change. With SimPy 3.x now available, I decided to take the plunge and try to build some DES models for an ongoing research project involving simulating a simplified inpatient obstetrical patient flow network. I originally built the model in Simio. I describe the project in a bit of detail in a previous post on using the R package caret for building metamodels of simulation output.

I decided to start by building SimPy models of overly simplified versions of this patient flow network. As my SimPy knowledge grew, I'd add more features to the model and overhaul its architecture. My goal is be able to build models of a similar complexity to the Java based pneumatic tube model. I'll also explore the ecosystem of Python tools that might be used alongside SimPy for things like model animation, logging and monitoring of simulation runs, simulation experiment management and management of simulation input and output data.

If you aren't familiar with SimPy, visit the main site and check out the documentation. It's got an overview, installation instructions, a basic tutorial, topic guide, example models and the API. Another site with some very nice network models and tutorials is available from Grotto Networking. Models in SimPy can range from purely process based with no object oriented (OO) features to almost fully OO based models.

Note: This series of blog posts assumes some knowledge of DES modeling, Python, and SimPy. I'm not trying to teach you how to use SimPy in general, but instead to share my journey of building progressively move complex DES patient flow models in SimPy. In the end, I hope to end up with a nice set of models and accompanying explanatory text that others can use and build on. I highly recommend the SimPy in 10 Minutes introduction in the official docs.

OB patient flow system to model

Here's a picture of the simple system I'll be modeling.

For this first set of simple models, features/assumptions of the system include:

  • Patients in labor arrive according to a poisson process
  • Assume c-section rate is zero and that all patients flow from OBS --> LDR --> PP units.
  • Length of stay at each unit is exponentially distributed with a given unit specific mean
  • Each unit has capacity level (i.e. the number of beds)
  • Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.

I'm going to start by building the simplest model possible and then start layering on complexity and features.

obflow_1: A Hello World level SimPy patient flow model

For my first ever SimPy model I decided to model the system as a sequence of delays. In other words, even though SimPy has features for modeling capacitated resources, I'm going to ignore them for now. This model is entirely process oriented. SimPy models processes with generator functions that can yield for events such as time durations or some system state condition. No classes are created in this first model. You'll also note that all the functions take a simulation environment instance as their first argurment. This is what you do in SimPy. It handles the clock and next event list logic as well as providing the means for running or single stepping through your model.

The model prints out a message as each patient enters each delay stage. No summary stats are computed.

In [2]:
import simpy
import numpy as np
from numpy.random import RandomState

"""
Simple OB patient flow model - NOT OO

Details:

- Generate arrivals via Poisson process
- Simple delays modeling stay in OBS, LDR, and PP units. No capacitated resource contention modeled.
- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.

"""

# Arrival rate and length of stay inputs.
ARR_RATE = 0.4
MEAN_LOS_OBS = 3
MEAN_LOS_LDR = 12
MEAN_LOS_PP = 48

RNG_SEED = 6353

def source(env, arr_rate, prng=RandomState(0)):
    """Source generates patients according to a simple Poisson process

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        arr_rate : float
            exponential arrival rate
        prng : RandomState object
            Seeded RandomState object for generating pseudo-random numbers.
            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html

    """

    patients_created = 0

    # Infinite loop for generatirng patients according to a poisson process.
    while True:

        # Generate next interarrival time
        iat = prng.exponential(1.0 / arr_rate)

        # Generate length of stay in each unit for this patient
        los_obs = prng.exponential(MEAN_LOS_OBS)
        los_ldr = prng.exponential(MEAN_LOS_LDR)
        los_pp = prng.exponential(MEAN_LOS_PP)

        # Update counter of patients
        patients_created += 1

        # Create a new patient flow process.
        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),
                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)

        # Register the process with the simulation environment
        env.process(obp)

        # This process will now yield to a 'timeout' event. This process will resume after iat time units.
        yield env.timeout(iat)

# Define an obpatient_flow "process"
def obpatient_flow(env, name, los_obs, los_ldr, los_pp):
    """Process function modeling how a patient flows through system.

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        name : str
            process instance id
        los_obs : float
            length of stay in OBS unit
        los_ldr : float
            length of stay in LDR unit
        los_pp : float
            length of stay in PP unit
    """

    # Note how we are simply modeling each stay as a delay. There
    # is NO contention for any resources.
    print("{} entering OBS at {:.4f}".format(name, env.now))
    yield env.timeout(los_obs)

    print("{} entering LDR at {:.4f}".format(name, env.now))
    yield env.timeout(los_ldr)

    print("{} entering PP at {:.4f}".format(name, env.now))
    yield env.timeout(los_pp)


# Initialize a simulation environment
env = simpy.Environment()

# Initialize a random number generator.
# See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html
prng = RandomState(RNG_SEED)

# Create a process generator and start it and add it to the env
# Calling obpatient(env) creates the generator.
# env.process() starts and adds it to env
runtime = 25
env.process(source(env, ARR_RATE, prng))

# Run the simulation
env.run(until=runtime)
Patient1 entering OBS at 0.0000
Patient1 entering LDR at 0.3475
Patient2 entering OBS at 3.7668
Patient3 entering OBS at 5.6439
Patient1 entering PP at 6.1449
Patient4 entering OBS at 7.2625
Patient4 entering LDR at 7.3175
Patient4 entering PP at 7.8836
Patient5 entering OBS at 8.9829
Patient3 entering LDR at 10.8140
Patient2 entering LDR at 12.9060
Patient6 entering OBS at 16.8153
Patient5 entering LDR at 16.9444
Patient7 entering OBS at 17.8737
Patient7 entering LDR at 18.0987
Patient7 entering PP at 18.1332
Patient8 entering OBS at 18.7157
Patient6 entering LDR at 19.4866
Patient9 entering OBS at 20.4913
Patient8 entering LDR at 22.8456
Patient8 entering PP at 23.1882
Patient5 entering PP at 24.5646

This first simple model illustrates the basic structure of a SimPy model with a single process function. We see how the simulation environment is needed for orchestrating the process as it yields and resumes. Now let's move on to a slightly more complex model in which the OBS, LDR, and PP units are modeled as Resource objects.

obflow_2: OBS unit modeled as a Resource

To start, let's model just the OBS unit as a Resource - LDR and PP units will still be modeled with delays. By doing this we can make sure we get the request and release logic working for one unit before duplicating it for the other two units. Like obflow_1, this model is process oriented.

In [5]:
import simpy
import numpy as np
from numpy.random import RandomState

"""
Simple OB patient flow model 2 - NOT OO

Details:

- Generate arrivals via Poisson process
- Uses one Resource objects to model OBS, the other units are modeled as simple delays.
- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.

"""
# Arrival rate and length of stay inputs.
ARR_RATE = 0.4
MEAN_LOS_OBS = 3
MEAN_LOS_LDR = 12
MEAN_LOS_PP = 48

CAPACITY_OBS = 2

RNG_SEED = 6353

def patient_generator(env, arr_rate, prng=RandomState(0)):
    """Generates patients according to a simple Poisson process

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        arr_rate : float
            exponential arrival rate
        prng : RandomState object
            Seeded RandomState object for generating pseudo-random numbers.
            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html

    """

    patients_created = 0

    # Infinite loop for generatirng patients according to a poisson process.
    while True:

        # Generate next interarrival time
        iat = prng.exponential(1.0 / arr_rate)

        # Generate length of stay in each unit for this patient
        los_obs = prng.exponential(MEAN_LOS_OBS)
        los_ldr = prng.exponential(MEAN_LOS_LDR)
        los_pp = prng.exponential(MEAN_LOS_PP)

        # Update counter of patients
        patients_created += 1

        # Create a new patient flow process.
        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),
                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)

        # Register the process with the simulation environment
        env.process(obp)

        # This process will now yield to a 'timeout' event. This process will resume after iat time units.
        yield env.timeout(iat)


def obpatient_flow(env, name, los_obs, los_ldr, los_pp):
    """Process function modeling how a patient flows through system.

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        name : str
            process instance id
        los_obs : float
            length of stay in OBS unit
        los_ldr : float
            length of stay in LDR unit
        los_pp : float
            length of stay in PP unit
    """

    print("{} trying to get OBS at {:.4f}".format(name, env.now))

    # Timestamp when patient tried to get OBS bed
    bed_request_ts = env.now
    # Request an obs bed
    bed_request = obs_unit.request()
    # Yield this process until a bed is available
    yield bed_request

    # We got an OBS bed
    print("{} entering OBS at {:.4f}".format(name, env.now))
    # Let's see if we had to wait to get the bed.
    if env.now > bed_request_ts:
        print("{} waited {:.4f} time units for OBS bed".format(name, env.now - bed_request_ts))

    # Yield this process again. Now wait until our length of stay elapses.
    # This is the actual stay in the bed
    yield env.timeout(los_obs)

    # All done with OBS, release the bed. Note that we pass the bed_request object
    # to the release() function so that the correct unit of the resource is released.
    obs_unit.release(bed_request)
    print("{} leaving OBS at {:.4f}".format(name, env.now))

    # Continue on through LDR and PP; modeled as simple delays for now.

    print("{} entering LDR at {:.4f}".format(name, env.now))
    yield env.timeout(los_ldr)

    print("{} entering PP at {:.4f}".format(name, env.now))
    yield env.timeout(los_pp)


# Initialize a simulation environment
env = simpy.Environment()

# Initialize a random number generator.
# See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html
prng = RandomState(RNG_SEED)

# Declare a Resource to model OBS unit. Default capacity is 1, we pass in desired capacity.
obs_unit = simpy.Resource(env, CAPACITY_OBS)

# Run the simulation for a while
runtime = 75
env.process(patient_generator(env, ARR_RATE, prng))
env.run(until=runtime)
Patient1 trying to get OBS at 0.0000
Patient1 entering OBS at 0.0000
Patient1 leaving OBS at 0.3475
Patient1 entering LDR at 0.3475
Patient2 trying to get OBS at 3.7668
Patient2 entering OBS at 3.7668
Patient3 trying to get OBS at 5.6439
Patient3 entering OBS at 5.6439
Patient1 entering PP at 6.1449
Patient4 trying to get OBS at 7.2625
Patient5 trying to get OBS at 8.9829
Patient3 leaving OBS at 10.8140
Patient3 entering LDR at 10.8140
Patient4 entering OBS at 10.8140
Patient4 waited 3.5515 time units for OBS bed
Patient4 leaving OBS at 10.8690
Patient4 entering LDR at 10.8690
Patient5 entering OBS at 10.8690
Patient5 waited 1.8861 time units for OBS bed
Patient4 entering PP at 11.4351
Patient2 leaving OBS at 12.9060
Patient2 entering LDR at 12.9060
Patient6 trying to get OBS at 16.8153
Patient6 entering OBS at 16.8153
Patient7 trying to get OBS at 17.8737
Patient8 trying to get OBS at 18.7157
Patient5 leaving OBS at 18.8306
Patient5 entering LDR at 18.8306
Patient7 entering OBS at 18.8306
Patient7 waited 0.9569 time units for OBS bed
Patient7 leaving OBS at 19.0556
Patient7 entering LDR at 19.0556
Patient8 entering OBS at 19.0556
Patient8 waited 0.3400 time units for OBS bed
Patient7 entering PP at 19.0901
Patient6 leaving OBS at 19.4866
Patient6 entering LDR at 19.4866
Patient9 trying to get OBS at 20.4913
Patient9 entering OBS at 20.4913
Patient8 leaving OBS at 23.1856
Patient8 entering LDR at 23.1856
Patient8 entering PP at 23.5282
Patient2 entering PP at 25.7847
Patient5 entering PP at 26.4507
Patient10 trying to get OBS at 26.8687
Patient10 entering OBS at 26.8687
Patient10 leaving OBS at 27.6307
Patient10 entering LDR at 27.6307
Patient11 trying to get OBS at 27.6876
Patient11 entering OBS at 27.6876
Patient6 entering PP at 28.3114
Patient9 leaving OBS at 30.4208
Patient9 entering LDR at 30.4208
Patient11 leaving OBS at 30.4652
Patient11 entering LDR at 30.4652
Patient11 entering PP at 30.5128
Patient12 trying to get OBS at 31.4872
Patient12 entering OBS at 31.4872
Patient13 trying to get OBS at 31.9981
Patient13 entering OBS at 31.9981
Patient12 leaving OBS at 32.2156
Patient12 entering LDR at 32.2156
Patient13 leaving OBS at 32.5002
Patient13 entering LDR at 32.5002
Patient14 trying to get OBS at 33.5398
Patient14 entering OBS at 33.5398
Patient12 entering PP at 35.1766
Patient15 trying to get OBS at 35.6094
Patient15 entering OBS at 35.6094
Patient13 entering PP at 35.8941
Patient16 trying to get OBS at 36.5130
Patient17 trying to get OBS at 36.5369
Patient15 leaving OBS at 36.6000
Patient15 entering LDR at 36.6000
Patient16 entering OBS at 36.6000
Patient16 waited 0.0870 time units for OBS bed
Patient16 leaving OBS at 37.4484
Patient16 entering LDR at 37.4484
Patient17 entering OBS at 37.4484
Patient17 waited 0.9115 time units for OBS bed
Patient18 trying to get OBS at 39.6783
Patient16 entering PP at 39.7335
Patient17 leaving OBS at 39.9760
Patient17 entering LDR at 39.9760
Patient18 entering OBS at 39.9760
Patient18 waited 0.2977 time units for OBS bed
Patient9 entering PP at 40.7127
Patient14 leaving OBS at 41.0861
Patient14 entering LDR at 41.0861
Patient19 trying to get OBS at 42.1910
Patient19 entering OBS at 42.1910
Patient10 entering PP at 42.5130
Patient18 leaving OBS at 42.7172
Patient18 entering LDR at 42.7172
Patient17 entering PP at 44.1436
Patient18 entering PP at 44.4498
Patient15 entering PP at 44.6092
Patient20 trying to get OBS at 44.6790
Patient20 entering OBS at 44.6790
Patient3 entering PP at 44.9184
Patient20 leaving OBS at 45.2709
Patient20 entering LDR at 45.2709
Patient21 trying to get OBS at 46.6338
Patient21 entering OBS at 46.6338
Patient21 leaving OBS at 47.8777
Patient21 entering LDR at 47.8777
Patient19 leaving OBS at 48.7520
Patient19 entering LDR at 48.7520
Patient14 entering PP at 49.3856
Patient19 entering PP at 49.7737
Patient22 trying to get OBS at 49.9682
Patient22 entering OBS at 49.9682
Patient22 leaving OBS at 50.6934
Patient22 entering LDR at 50.6934
Patient21 entering PP at 50.6969
Patient23 trying to get OBS at 51.9160
Patient23 entering OBS at 51.9160
Patient24 trying to get OBS at 52.8406
Patient24 entering OBS at 52.8406
Patient24 leaving OBS at 53.3548
Patient24 entering LDR at 53.3548
Patient23 leaving OBS at 54.9689
Patient23 entering LDR at 54.9689
Patient25 trying to get OBS at 54.9893
Patient25 entering OBS at 54.9893
Patient24 entering PP at 55.3316
Patient20 entering PP at 55.7203
Patient23 entering PP at 57.8120
Patient22 entering PP at 58.5375
Patient26 trying to get OBS at 61.0166
Patient26 entering OBS at 61.0166
Patient26 leaving OBS at 61.1124
Patient26 entering LDR at 61.1124
Patient27 trying to get OBS at 62.3948
Patient27 entering OBS at 62.3948
Patient25 leaving OBS at 63.4644
Patient25 entering LDR at 63.4644
Patient27 leaving OBS at 63.5210
Patient27 entering LDR at 63.5210
Patient26 entering PP at 65.5079
Patient28 trying to get OBS at 68.1463
Patient28 entering OBS at 68.1463
Patient29 trying to get OBS at 68.4204
Patient29 entering OBS at 68.4204
Patient25 entering PP at 68.5573
Patient28 leaving OBS at 68.5577
Patient28 entering LDR at 68.5577
Patient30 trying to get OBS at 72.7648
Patient30 entering OBS at 72.7648
Patient30 leaving OBS at 72.7858
Patient30 entering LDR at 72.7858
Patient31 trying to get OBS at 72.8227
Patient31 entering OBS at 72.8227
Patient27 entering PP at 73.2677
Patient29 leaving OBS at 73.3031
Patient29 entering LDR at 73.3031
Patient31 leaving OBS at 73.5402
Patient31 entering LDR at 73.5402
Patient32 trying to get OBS at 73.7931
Patient32 entering OBS at 73.7931
Patient30 entering PP at 74.6737

Seems to be working. Note that with an arrival rate of 0.4 and a runtime of 75, we expect about...

In [6]:
print("Expected number of arrivals = {:.2f}".format(ARR_RATE * runtime))
Expected number of arrivals = 30.00

We got 32 with this random number seed. No cause for alarm. Of course, we should run this thing longer and compute some stats so that we can make sure things are working fine. We'll do that shortly.

obflow_3: All units modeled as a Resources

Let's build on obflow_2 by modeling all three units as resources. Since these are hospital beds, there are no queues per se. If a patient finishes their stay in OBS and there are no beds available in the LDR, they wait in OBS. Likewise, when trying to get into PP after delivering the baby in the LDR, they may get blocked in the LDR waiting for an available PP bed. In more complex models we may actually do length of stay adjustments for those patients who are blocked. For now, let's just keep it simple and not do any such adjustments.

In SimPy terms, we have to be careful that we don't release a bed we are in until we have successfully reserved a bed in the next unit we are to visit. As you'll see, this leads to quite a bit of repetitive code in our process oriented approach. This will motivate us to start exploring object oriented models.

Finally, we'll improve our patient generator by adding the ability to model an initial delay before patient generation should start as well as to specify a stop time for the patient generation process (which can be different than the simulation model should stop).

In [9]:
import simpy
import numpy as np
from numpy.random import RandomState

"""
Simple OB patient flow model 3 - NOT OO

Details:

- Generate arrivals via Poisson process
- Uses one Resource objects to model OBS, LDR, and PP.
- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.
- Additional functionality added to arrival generator (initial delay and arrival stop time).

"""
ARR_RATE = 0.4
MEAN_LOS_OBS = 3
MEAN_LOS_LDR = 12
MEAN_LOS_PP = 48

CAPACITY_OBS = 2
CAPACITY_LDR = 6
CAPACITY_PP = 24

RNG_SEED = 6353

def patient_generator(env, arr_stream, arr_rate, initial_delay=0,
                      stoptime=simpy.core.Infinity, prng=RandomState(0)):
    """Generates patients according to a simple Poisson process

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        arr_rate : float
            exponential arrival rate
        initial_delay: float (default 0)
            time before arrival generation should begin
        stoptime: float (default Infinity)
            time after which no arrivals are generated
        prng : RandomState object
            Seeded RandomState object for generating pseudo-random numbers.
            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html

    """

    patients_created = 0

    # Yield for the initial delay
    yield env.timeout(initial_delay)

    # Generate arrivals as long as simulation time is before stoptime
    while env.now < stoptime:

        iat = prng.exponential(1.0 / arr_rate)

        # Sample los distributions
        los_obs = prng.exponential(MEAN_LOS_OBS)
        los_ldr = prng.exponential(MEAN_LOS_LDR)
        los_pp = prng.exponential(MEAN_LOS_PP)


        # Create new patient process instance
        patients_created += 1
        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),
                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)

        env.process(obp)

        # Compute next interarrival time

        yield env.timeout(iat)


def obpatient_flow(env, name, los_obs, los_ldr, los_pp):
    """Process function modeling how a patient flows through system.

        Parameters
        ----------
        env : simpy.Environment
            the simulation environment
        name : str
            process instance id
        los_obs : float
            length of stay in OBS unit
        los_ldr : float
            length of stay in LDR unit
        los_pp : float
            length of stay in PP unit
    """

    # Note the repetitive code and the use of separate request objects for each
    # stay in the different units.

    # OBS
    print("{} trying to get OBS at {:.4f}".format(name, env.now))
    bed_request_ts = env.now
    bed_request1 = obs_unit.request() # Request an OBS bed
    yield bed_request1
    print("{} entering OBS at {:.4f}".format(name, env.now))
    if env.now > bed_request_ts:
        print("{} waited {:.4f} time units for OBS bed".format(name, env.now-  bed_request_ts))
    yield env.timeout(los_obs) # Stay in obs bed

    print("{} trying to get LDR at {:.4f}".format(name, env.now))
    bed_request_ts = env.now
    bed_request2 = ldr_unit.request()  # Request an LDR bed
    yield bed_request2

    # Got LDR bed, release OBS bed
    obs_unit.release(bed_request1)  # Release the OBS bed
    print("{} leaving OBS at {}".format(name, env.now))

    # LDR stay
    print("{} entering LDR at {:.4f}".format(name, env.now))
    if env.now > bed_request_ts:
        print("{} waited {:.4f} time units for LDR bed".format(name, env.now - bed_request_ts))
    yield env.timeout(los_ldr) # Stay in LDR bed

    print("{} trying to get PP at {:.4f}".format(name, env.now))
    bed_request_ts = env.now
    bed_request3 = pp_unit.request()  # Request a PP bed
    yield bed_request3

    # Got PP bed, release LDR bed
    ldr_unit.release(bed_request2)  # Release the obs bed
    print("{} leaving LDR at {:.4f}".format(name, env.now))

    # PP stay
    print("{} entering PP at {:.4f}".format(name, env.now))
    if env.now > bed_request_ts:
        print("{} waited {:.4f} time units for PP bed".format(name, env.now - bed_request_ts))
    yield env.timeout(los_pp) # Stay in PP bed
    pp_unit.release(bed_request3)  # Release the PP bed

    print("{} leaving PP and system at {:.4f}".format(name, env.now))

# Initialize a simulation environment
env = simpy.Environment()

prng = RandomState(RNG_SEED)

rho_obs = ARR_RATE * MEAN_LOS_OBS / CAPACITY_OBS
rho_ldr = ARR_RATE * MEAN_LOS_LDR / CAPACITY_LDR
rho_pp = ARR_RATE * MEAN_LOS_PP / CAPACITY_PP

# Declare Resources to model all units
obs_unit = simpy.Resource(env, CAPACITY_OBS)
ldr_unit = simpy.Resource(env, CAPACITY_LDR)
pp_unit = simpy.Resource(env, CAPACITY_PP)

# Run the simulation for a while. Let's shut arrivals off after 50 time units.
runtime = 75
stop_arrivals = 50
env.process(patient_generator(env, "Type1", ARR_RATE, 0, stop_arrivals, prng))
env.run(until=runtime)
Patient1 trying to get OBS at 0.0000
Patient1 entering OBS at 0.0000
Patient1 trying to get LDR at 0.3475
Patient1 leaving OBS at 0.3474891089551544
Patient1 entering LDR at 0.3475
Patient2 trying to get OBS at 3.7668
Patient2 entering OBS at 3.7668
Patient3 trying to get OBS at 5.6439
Patient3 entering OBS at 5.6439
Patient1 trying to get PP at 6.1449
Patient1 leaving LDR at 6.1449
Patient1 entering PP at 6.1449
Patient4 trying to get OBS at 7.2625
Patient5 trying to get OBS at 8.9829
Patient3 trying to get LDR at 10.8140
Patient3 leaving OBS at 10.813985162144862
Patient3 entering LDR at 10.8140
Patient4 entering OBS at 10.8140
Patient4 waited 3.5515 time units for OBS bed
Patient4 trying to get LDR at 10.8690
Patient4 leaving OBS at 10.869006908518232
Patient4 entering LDR at 10.8690
Patient5 entering OBS at 10.8690
Patient5 waited 1.8861 time units for OBS bed
Patient4 trying to get PP at 11.4351
Patient4 leaving LDR at 11.4351
Patient4 entering PP at 11.4351
Patient2 trying to get LDR at 12.9060
Patient2 leaving OBS at 12.90600490911474
Patient2 entering LDR at 12.9060
Patient6 trying to get OBS at 16.8153
Patient6 entering OBS at 16.8153
Patient7 trying to get OBS at 17.8737
Patient8 trying to get OBS at 18.7157
Patient5 trying to get LDR at 18.8306
Patient5 leaving OBS at 18.830564537083035
Patient5 entering LDR at 18.8306
Patient7 entering OBS at 18.8306
Patient7 waited 0.9569 time units for OBS bed
Patient7 trying to get LDR at 19.0556
Patient7 leaving OBS at 19.05561983104953
Patient7 entering LDR at 19.0556
Patient8 entering OBS at 19.0556
Patient8 waited 0.3400 time units for OBS bed
Patient7 trying to get PP at 19.0901
Patient7 leaving LDR at 19.0901
Patient7 entering PP at 19.0901
Patient6 trying to get LDR at 19.4866
Patient6 leaving OBS at 19.48663625825889
Patient6 entering LDR at 19.4866
Patient9 trying to get OBS at 20.4913
Patient9 entering OBS at 20.4913
Patient8 trying to get LDR at 23.1856
Patient8 leaving OBS at 23.185576349161767
Patient8 entering LDR at 23.1856
Patient8 trying to get PP at 23.5282
Patient8 leaving LDR at 23.5282
Patient8 entering PP at 23.5282
Patient2 trying to get PP at 25.7847
Patient2 leaving LDR at 25.7847
Patient2 entering PP at 25.7847
Patient5 trying to get PP at 26.4507
Patient5 leaving LDR at 26.4507
Patient5 entering PP at 26.4507
Patient10 trying to get OBS at 26.8687
Patient10 entering OBS at 26.8687
Patient10 trying to get LDR at 27.6307
Patient10 leaving OBS at 27.63070937211375
Patient10 entering LDR at 27.6307
Patient11 trying to get OBS at 27.6876
Patient11 entering OBS at 27.6876
Patient6 trying to get PP at 28.3114
Patient6 leaving LDR at 28.3114
Patient6 entering PP at 28.3114
Patient9 trying to get LDR at 30.4208
Patient9 leaving OBS at 30.420840322655998
Patient9 entering LDR at 30.4208
Patient11 trying to get LDR at 30.4652
Patient11 leaving OBS at 30.465221449496795
Patient11 entering LDR at 30.4652
Patient11 trying to get PP at 30.5128
Patient11 leaving LDR at 30.5128
Patient11 entering PP at 30.5128
Patient12 trying to get OBS at 31.4872
Patient12 entering OBS at 31.4872
Patient13 trying to get OBS at 31.9981
Patient13 entering OBS at 31.9981
Patient12 trying to get LDR at 32.2156
Patient12 leaving OBS at 32.215634051591636
Patient12 entering LDR at 32.2156
Patient13 trying to get LDR at 32.5002
Patient13 leaving OBS at 32.500238016146994
Patient13 entering LDR at 32.5002
Patient14 trying to get OBS at 33.5398
Patient14 entering OBS at 33.5398
Patient7 leaving PP and system at 34.5316
Patient12 trying to get PP at 35.1766
Patient12 leaving LDR at 35.1766
Patient12 entering PP at 35.1766
Patient15 trying to get OBS at 35.6094
Patient15 entering OBS at 35.6094
Patient13 trying to get PP at 35.8941
Patient13 leaving LDR at 35.8941
Patient13 entering PP at 35.8941
Patient16 trying to get OBS at 36.5130
Patient17 trying to get OBS at 36.5369
Patient15 trying to get LDR at 36.6000
Patient15 leaving OBS at 36.59999969293714
Patient15 entering LDR at 36.6000
Patient16 entering OBS at 36.6000
Patient16 waited 0.0870 time units for OBS bed
Patient16 trying to get LDR at 37.4484
Patient16 leaving OBS at 37.44839668463982
Patient16 entering LDR at 37.4484
Patient17 entering OBS at 37.4484
Patient17 waited 0.9115 time units for OBS bed
Patient11 leaving PP and system at 38.1145
Patient18 trying to get OBS at 39.6783
Patient16 trying to get PP at 39.7335
Patient16 leaving LDR at 39.7335
Patient16 entering PP at 39.7335
Patient17 trying to get LDR at 39.9760
Patient17 leaving OBS at 39.97599817647182
Patient17 entering LDR at 39.9760
Patient18 entering OBS at 39.9760
Patient18 waited 0.2977 time units for OBS bed
Patient9 trying to get PP at 40.7127
Patient9 leaving LDR at 40.7127
Patient9 entering PP at 40.7127
Patient14 trying to get LDR at 41.0861
Patient14 leaving OBS at 41.086116528209494
Patient14 entering LDR at 41.0861
Patient19 trying to get OBS at 42.1910
Patient19 entering OBS at 42.1910
Patient10 trying to get PP at 42.5130
Patient10 leaving LDR at 42.5130
Patient10 entering PP at 42.5130
Patient18 trying to get LDR at 42.7172
Patient18 leaving OBS at 42.71719377930215
Patient18 entering LDR at 42.7172
Patient1 leaving PP and system at 42.8634
Patient17 trying to get PP at 44.1436
Patient17 leaving LDR at 44.1436
Patient17 entering PP at 44.1436
Patient18 trying to get PP at 44.4498
Patient18 leaving LDR at 44.4498
Patient18 entering PP at 44.4498
Patient15 trying to get PP at 44.6092
Patient15 leaving LDR at 44.6092
Patient15 entering PP at 44.6092
Patient20 trying to get OBS at 44.6790
Patient20 entering OBS at 44.6790
Patient3 trying to get PP at 44.9184
Patient3 leaving LDR at 44.9184
Patient3 entering PP at 44.9184
Patient20 trying to get LDR at 45.2709
Patient20 leaving OBS at 45.270857202928546
Patient20 entering LDR at 45.2709
Patient21 trying to get OBS at 46.6338
Patient21 entering OBS at 46.6338
Patient21 trying to get LDR at 47.8777
Patient21 leaving OBS at 47.8776934769442
Patient21 entering LDR at 47.8777
Patient19 trying to get LDR at 48.7520
Patient19 leaving OBS at 48.75197164898342
Patient19 entering LDR at 48.7520
Patient14 trying to get PP at 49.3856
Patient14 leaving LDR at 49.3856
Patient14 entering PP at 49.3856
Patient19 trying to get PP at 49.7737
Patient19 leaving LDR at 49.7737
Patient19 entering PP at 49.7737
Patient22 trying to get OBS at 49.9682
Patient22 entering OBS at 49.9682
Patient22 trying to get LDR at 50.6934
Patient22 leaving OBS at 50.69343889384862
Patient22 entering LDR at 50.6934
Patient21 trying to get PP at 50.6969
Patient21 leaving LDR at 50.6969
Patient21 entering PP at 50.6969
Patient20 trying to get PP at 55.7203
Patient20 leaving LDR at 55.7203
Patient20 entering PP at 55.7203
Patient16 leaving PP and system at 56.0439
Patient18 leaving PP and system at 57.0161
Patient12 leaving PP and system at 58.3965
Patient22 trying to get PP at 58.5375
Patient22 leaving LDR at 58.5375
Patient22 entering PP at 58.5375
Patient10 leaving PP and system at 59.8131
Patient9 leaving PP and system at 60.4868
Patient5 leaving PP and system at 63.4959
Patient17 leaving PP and system at 69.8364
Patient14 leaving PP and system at 69.9692
Patient15 leaving PP and system at 70.5547
Patient22 leaving PP and system at 72.2739
Patient20 leaving PP and system at 73.4076

Notice we get far fewer patients because we turned off the arrival spigot at time 50.

As a veteran simulation modeler, I'm not really liking this process oriented approach. It's not going to scale well as we add more units (e.g. if we were modeling all the inpatient beds) and more routing complexity. Having built a pretty big object oriented DES model (in Java with Simkit) for analyzing pneumatic tube systems, I know that that's the way to go. In the next post, I'll share my first object oriented SimPy model. Not only does an object oriented approach let us avoid so much code repitition, it also makes collecting stats and logging easier.