import simpy
from numpy.random import RandomState
"""
Simple OB patient flow model 5 - Very simple OO
Details:
- Generate arrivals via Poisson process
- Define an OBUnit class that contains a simpy.Resource object as a member.
Not subclassing Resource, just trying to use it as a member.
- Routing is done via setting ``out`` member of an OBUnit instance to
another OBUnit instance to which the OB patient flow instance should be
routed. The routing logic, for now, is in OBUnit object. Later,
we need some sort of router object and data driven routing.
- Trying to get patient flow working without a process function that
explicitly articulates the sequence of units and stays.
Key Lessons Learned:
- Any function that is a generator and might potentially yield for an event
must get registered as a process.
"""
# Arrival rate and length of stay inputs.
ARR_RATE = 0.4
MEAN_LOS_OBS = 3
MEAN_LOS_LDR = 12
MEAN_LOS_PP = 48
# Unit capacities
CAPACITY_OBS = 2
CAPACITY_LDR = 6
CAPACITY_PP = 24
# Random number seed
RNG_SEED = 6353An object oriented SimPy patient flow simulation model
In Part 2 of this series we build a simple object oriented SimPy model for simulating inpatient obsetrics patient flow.
python
simpy
simulation
healthcare
patient flow
In Part 1 of this series, we built a few simple patient flow simulation models using SimPy. They all were process oriented. In this post, I’ll share my first attempt at an object oriented (OO) SimPy model for the same system. As mentioned in my previous post, I’ve had good luck building large complex discete event simulation models (including patient flow) with the object oriented Java library Simkit. In the process of building these initial OO models, I learned a bunch from the posts done by the folks at Grotto Networking.
OB patient flow system to model
As a reminder, 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_5_oo_2: A simple object oriented SimPy patient flow model
Since I’m new to SimPy, this first object oriented model will be relatively simple. There are only four classes:
OBPatientGenerator- generates patientsOBPatient- a patient; generated byOBPatientGeneratorOBUnit- for modeling the OBS, LDR, and PP unitsExitFlow- sink to which patients are routed when ready to leave the system
My goal was to just figure out how to model a simple static route, OBS -> LDR -> PP -> Exit, using SimPy in an OO way. Let’s look at each class and then at the main driver script. You can get the code for this and the other models at https://github.com/misken/obsimpy.
Imports and constants
The OBPatientGenerator class
This class generates OB patients according to a simple poisson process (exponential interarrival times). Its run() method gets registered as a SimPy process. The main loop in run() yields until the time for the next arrival. Upon arrival, a new OBPatient instance is generated and a request is initiated for a bed in the first unit of the patient’s route - OBS in this case. We’ll see shortly how the routing information is included in the OBPatient class (a “hack” that we’ll remedy later with proper routing related classes).
A few niceties were added to the class based on things I know will be useful to have later. These include:
- arrival stream and patient type attributes; future versions of this model will support multiple patient types (e.g. spontaneous labor and induced labor) and arrival streams (e.g. random and scheduled),
- ability to shut off arrivals based on simulation time or number of patients generated.
class OBPatientGenerator(object):
""" Generates patients.
Set the "out" member variable to resource at which patient generated.
Parameters
----------
env : simpy.Environment
the simulation environment
arr_rate : float
Poisson arrival rate (expected number of arrivals per unit time)
arr_stream : int
Arrival stream id (default 1). Currently there is just one arrival
stream corresponding to the one patient generator class. In future,
likely to be be multiple generators for generating random and
scheduled arrivals
initial_delay : float
Starts generation after an initial delay. (default 0.0)
stoptime : float
Stops generation at the stoptime. (default Infinity)
max_arrivals : int
Stops generation after max_arrivals. (default Infinity)
debug: bool
If True, status message printed after
each patient created. (default False)
"""
def __init__(self, env, arr_rate, patient_type=1, arr_stream=1,
initial_delay=0,
stoptime=simpy.core.Infinity,
max_arrivals=simpy.core.Infinity, debug=False):
self.env = env
self.arr_rate = arr_rate
self.patient_type = patient_type
self.arr_stream = arr_stream
self.initial_delay = initial_delay
self.stoptime = stoptime
self.max_arrivals = max_arrivals
self.debug = debug
self.out = None
self.num_patients_created = 0
self.prng = RandomState(RNG_SEED)
self.action = env.process(self.run()) # starts the run() method as a SimPy process
def run(self):
"""The patient generator.
"""
# Delay for initial_delay
yield self.env.timeout(self.initial_delay)
# Main generator loop that terminates when stoptime reached
while self.env.now < self.stoptime and \
self.num_patients_created < self.max_arrivals:
# Delay until time for next arrival
# Compute next interarrival time
iat = self.prng.exponential(1.0 / self.arr_rate)
yield self.env.timeout(iat)
self.num_patients_created += 1
# Create new patient
obp = OBPatient(self.env.now, self.num_patients_created,
self.patient_type, self.arr_stream,
prng=self.prng)
if self.debug:
print("Patient {} created at {:.2f}.".format(
self.num_patients_created, self.env.now))
# Set out member to OBunit object representing next destination
self.out = obunits[obp.planned_route_stop[1]]
# Initiate process of requesting first bed in route
self.env.process(self.out.put(obp))The OBPatient class
These are the patients. As mentioned above, routing info (for now) is hard coded into this class. Each patient has same static route of OBS -> LDR -> PP. Lists are used to store routing information. Python uses 0 based indexing but we’ve chosen to use the [1] spot for OBS, [2] for LDR, and [3] for PP quantities. The [0] slot isn’t used. Just wanted to avoid a bunch of current_stop_num - 1 type of code fragments. Makes the code more readable, especially since we often have to refer to the previous stop and wish to avoid things like current_stop_num - 2 to mean the previous stop. Obviously this is simple to change if you’d rather use 0 for stop number 1.
Several lists are used to store bed request instances as well as useful timestamps (any variable name ending in _ts such as entry and exit times to each unit.
class OBPatient(object):
""" Models an OB patient
Parameters
----------
arr_time : float
Patient arrival time
patient_id : int
Unique patient id
patient_type : int
Patient type id (default 1). Currently just one patient type.
In our prior research work we used a scheme with 11 patient types.
arr_stream : int
Arrival stream id (default 1). Currently there is just one arrival
stream corresponding to the one patient generator class. In future,
likely to be be multiple generators for generating random and
scheduled arrivals.
"""
def __init__(self, arr_time, patient_id, patient_type=1, arr_stream=1,
prng=RandomState(0)):
self.arr_time = arr_time
self.patient_id = patient_id
self.patient_type = patient_type
self.arr_stream = arr_stream
self.name = 'Patient_{}'.format(patient_id)
# Hard coding route, los and bed requests for now
# Not sure how best to do routing related data structures.
# Hack for now using combination of lists here, the out member
# and the obunits dictionary.
self.current_stay_num = 0
self.route_length = 3
self.planned_route_stop = []
self.planned_route_stop.append(None)
self.planned_route_stop.append("OBS")
self.planned_route_stop.append("LDR")
self.planned_route_stop.append("PP")
self.planned_route_stop.append("EXIT")
self.planned_los = []
self.planned_los.append(None)
self.planned_los.append(prng.exponential(MEAN_LOS_OBS))
self.planned_los.append(prng.exponential(MEAN_LOS_LDR))
self.planned_los.append(prng.exponential(MEAN_LOS_PP))
# Since we have fixed route for now, just initialize full list to
# hold bed requests
self.bed_requests = [None for _ in range(self.route_length + 1)]
self.request_entry_ts = [None for _ in range(self.route_length + 1)]
self.entry_ts = [None for _ in range(self.route_length + 1)]
self.exit_ts = [None for _ in range(self.route_length + 1)]
def __repr__(self):
return "patientid: {}, arr_stream: {}, time: {}". \
format(self.patient_id, self.arr_stream, self.arr_time)The OBUnit class
This class is used to model the OBS, LDR and PP units. It does most of the heavy lifting related to patients flowing into and out of the various units. It can be used to model both finite and infinite capacity resources. A few basic statistal accumulators are included. The out member gets set to the next unit to be visited by the patient.
class OBunit(object):
""" Models an OB unit with fixed capacity.
Parameters
----------
env : simpy.Environment
the simulation environment
name : str
unit name
capacity : integer (or None)
Number of beds. Use None for infinite capacity.
"""
def __init__(self, env, name, capacity=None, debug=False):
if capacity is None:
self.capacity = simpy.core.Infinity
else:
self.capacity = capacity
self.env = env
self.name = name
self.debug = debug
# Use a simpy Resource as one of the class members
self.unit = simpy.Resource(env, capacity)
# Statistical accumulators
self.num_entries = 0
self.num_exits = 0
self.tot_occ_time = 0.0
# The out member will get set to destination unit
self.out = None
def put(self, obp):
""" A process method called when a bed is requested in the unit.
The logic of this method is reminiscent of the routing logic
in the process oriented obflow models 1-3. However, this method
is used for all of the units - no repeated logic.
Parameters
----------
env : simpy.Environment
the simulation environment
obp : OBPatient object
the patient requestion the bed
"""
if self.debug:
print("{} trying to get {} at {:.4f}".format(obp.name,
self.name, self.env.now))
# Increments patient's attribute number of units visited
obp.current_stay_num += 1
# Timestamp of request time
bed_request_ts = self.env.now
# Request a bed
bed_request = self.unit.request()
# Store bed request and timestamp in patient's request lists
obp.bed_requests[obp.current_stay_num] = bed_request
obp.request_entry_ts[obp.current_stay_num] = self.env.now
# Yield until we get a bed
yield bed_request
# Seized a bed.
obp.entry_ts[obp.current_stay_num] = self.env.now
# Check if we have a bed from a previous stay and release it.
# Update stats for previous unit.
if obp.bed_requests[obp.current_stay_num - 1] is not None:
previous_request = obp.bed_requests[obp.current_stay_num - 1]
previous_unit_name = \
obp.planned_route_stop[obp.current_stay_num - 1]
previous_unit = obunits[previous_unit_name]
previous_unit.unit.release(previous_request)
previous_unit.num_exits += 1
previous_unit.tot_occ_time += \
self.env.now - obp.entry_ts[obp.current_stay_num - 1]
obp.exit_ts[obp.current_stay_num - 1] = self.env.now
if self.debug:
print("{} entering {} at {:.4f}".format(obp.name, self.name,
self.env.now))
self.num_entries += 1
if self.debug:
if self.env.now > bed_request_ts:
waittime = self.env.now - bed_request_ts
print("{} waited {:.4f} time units for {} bed".format(obp.name,
waittime,
self.name))
# Determine los and then yield for the stay
los = obp.planned_los[obp.current_stay_num]
yield self.env.timeout(los)
# Go to next destination (which could be an exitflow)
next_unit_name = obp.planned_route_stop[obp.current_stay_num + 1]
self.out = obunits[next_unit_name]
if obp.current_stay_num == obp.route_length:
# For ExitFlow object, no process needed
self.out.put(obp)
else:
# Process for putting patient into next bed
self.env.process(self.out.put(obp))
def basic_stats_msg(self):
""" Compute entries, exits, avg los and create summary message.
Returns
-------
str
Message with basic stats
"""
if self.num_exits > 0:
alos = self.tot_occ_time / self.num_exits
else:
alos = 0
msg = "{:6}:\t Entries={}, Exits={}, Occ={}, ALOS={:4.2f}".\
format(self.name, self.num_entries, self.num_exits,
self.unit.count, alos)
return msgThe ExitFlow class
This is a “sink” and is used as a place to trigger any post-processing needs after the patient has completed their entire route through the OB system. For now we are just doing some basic statistical and conservation of flow summaries. More importantly, this class handles releasing the last bed in the patient’s route.
class ExitFlow(object):
""" Patients routed here when ready to exit.
Patient objects put into a Store. Can be accessed later for stats
and logs. A little worried about how big the Store will get.
Parameters
----------
env : simpy.Environment
the simulation environment
debug : boolean
if true then patient details printed on arrival
"""
def __init__(self, env, name, store_obp=True, debug=False):
self.store = simpy.Store(env)
self.env = env
self.name = name
self.store_obp = store_obp
self.debug = debug
self.num_exits = 0
self.last_exit = 0.0
def put(self, obp):
if obp.bed_requests[obp.current_stay_num] is not None:
previous_request = obp.bed_requests[obp.current_stay_num]
previous_unit_name = obp.planned_route_stop[obp.current_stay_num]
previous_unit = obunits[previous_unit_name]
previous_unit.unit.release(previous_request)
previous_unit.num_exits += 1
previous_unit.tot_occ_time += self.env.now - obp.entry_ts[
obp.current_stay_num]
obp.exit_ts[obp.current_stay_num - 1] = self.env.now
self.last_exit = self.env.now
self.num_exits += 1
if self.debug:
print(obp)
# Store patient
if self.store_obp:
self.store.put(obp)
def basic_stats_msg(self):
""" Create summary message with basic stats on exits.
Returns
-------
str
Message with basic stats
"""
msg = "{:6}:\t Exits={}, Last Exit={:10.2f}".format(self.name,
self.num_exits,
self.last_exit)
return msgMain driver script
This looks much like the code in the previous models detailed in Part 1 of this series. One important change in this model is that the simulation environment variable is never used as a “global” variable. Any class needing access to the simulation environment explicitly requires it in its __init__ method.
# Initialize a simulation environment
simenv = simpy.Environment()
# Compute and display traffic intensities
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
print("rho_obs: {:6.3f}\nrho_ldr: {:6.3f}\nrho_pp: {:6.3f}".format(rho_obs,
rho_ldr,
rho_pp))
# Create nursing units
obs_unit = OBunit(simenv, 'OBS', CAPACITY_OBS, debug=False)
ldr_unit = OBunit(simenv, 'LDR', CAPACITY_LDR, debug=False)
pp_unit = OBunit(simenv, 'PP', CAPACITY_PP, debug=False)
# Define system exit
exitflow = ExitFlow(simenv, 'EXIT', store_obp=False)
# Create dictionary of units keyed by name. This object can be passed along
# to other objects so that the units are accessible as patients "flow".
obunits = {'OBS': obs_unit, 'LDR': ldr_unit, 'PP': pp_unit, 'EXIT': exitflow}
# Create a patient generator
obpat_gen = OBPatientGenerator(simenv, ARR_RATE, debug=False)
# Routing logic
# Currently routing logic is hacked into the OBPatientGenerator
# and OBPatient objects
# Run the simulation for a while
runtime = 1000
simenv.run(until=runtime)
# Patient generator stats
print("\nNum patients generated: {}\n".format(obpat_gen.num_patients_created))
# Unit stats
print(obs_unit.basic_stats_msg())
print(ldr_unit.basic_stats_msg())
print(pp_unit.basic_stats_msg())
# System exit stats
print("\nNum patients exiting system: {}\n".format(exitflow.num_exits))
print("Last exit at: {:.2f}\n".format(exitflow.last_exit))rho_obs: 0.600
rho_ldr: 0.800
rho_pp: 0.800
Num patients generated: 453
OBS : Entries=431, Exits=429, Occ=2, ALOS=4.27
LDR : Entries=429, Exits=425, Occ=4, ALOS=11.68
PP : Entries=425, Exits=402, Occ=23, ALOS=48.40
Num patients exiting system: 402
Last exit at: 998.69
Next Steps
Now that I have a feel for how SimPy works, I’ve got two parallel directions to pursue.
Task 1: Research quality SimPy model
I’ll start creating a more full featured version that replicates the functionality of the Simio based model I used in a recent research project. Some of the goals for this new version will be:
- volume, length of stay, routing information and other key inputs should be read in from structured data input files (probably YAML or JSON format),
- create routing classes to handle both random and scheduled arrivals,
- capability to generate time dependent Poisson arrivals,
- detailed patient flow logging to files,
- more statistical output analysis,
- comprehensive documentation,
- a test suite including validation via queueing models (for tractable cases).
Task 2: Rebuild this model using the R package, simmer
I’m a regular user of R for both research and teaching. So, I was pleasantly surprised to learn about a DES framework for R (!). The simmer package looks relatively new and seems to be positioning itself similarly to SimPy. Here’s the description from their main web page:
simmer is a process-oriented and trajectory-based Discrete-Event Simulation (DES) package for R. Designed to be a generic framework like SimPy or SimJulia, it leverages the power of Rcpp to boost the performance and turning DES in R feasible. As a noteworthy characteristic, simmer exploits the concept of trajectory: a common path in the simulation model for entities of the same type. It is pretty flexible and simple to use, and leverages the chaining/piping workflow introduced by the magrittr package.
Interesting. While I love chaining/piping via magrittr, especially when using the mighty dplyr package, I’m wondering if this will tie simmer to a process oriented simulation world view. I can see how the notion of a trajectory might lend itself to chaining and piping. Let’s find out. I’ll post the simmer version as soon as I get one working.
Reuse
Citation
BibTeX citation:
@online{isken2017,
author = {Mark Isken},
title = {An Object Oriented {SimPy} Patient Flow Simulation Model},
date = {2017-07-11},
langid = {en}
}
For attribution, please cite this work as:
Mark Isken. 2017. “An Object Oriented SimPy Patient Flow
Simulation Model.” July 11, 2017.