I've been working on a problem where I must simulate patients at a treatment centre in a Simpy 4.0.1 program, that has multiple categories of patients (say "A", "B" and "C"). The patients "arrive" based on timeouts generated from numpy.random.exponential() and "occupy a bed" (Simpy resources) at the centre based on numpy.random.lognormal(). Each patient type has a different set of values to control this sampling. This is a simple FIFO process and has no priority involved, the patients simply arrive and stay at different rates.
The first way I tried was to create 3 processes for each category of patient, and let Simpy handle their generation:
import simpy
from numpy.random import exponential, lognormal
counter = 0
total_beds = 5
run_length = 10
mean_iat_a, mean_iat_b, mean_iat_c = 1.0, 2.0, 3.0
normal_mean_a, normal_mean_b, normal_mean_c = 4.0, 5.0, 6.0
normal_stdev_a, normal_stdev_b, normal_stdev_c = 7.0, 8.0, 9.0
class Patient:
def __init__(self, pid, ptype):
self.pid = pid
self.ptype = ptype
self.time_arrived = None
self.time_admitted = None
self.time_discharged = None
def generate_arrivals(env, ptype):
while True:
global counter
counter += 1
patient = Patient(counter, ptype)
env.process(perform_treatment(env, patient))
# sample arrivals
if patient.ptype == 'A':
interarival_time = exponential(scale=mean_iat_a)
elif patient.ptype == 'B':
interarival_time = exponential(scale=mean_iat_b)
else:
interarival_time = exponential(scale=mean_iat_c)
yield env.timeout(interarival_time)
def perform_treatment(env, patient):
patient.time_arrived = env.now
print(f'Patient {patient.pid} - {patient.ptype} - arrived at {patient.time_arrived}')
with beds.request() as req:
yield req
patient.time_admitted = env.now
print(f'Patient {patient.pid} - {patient.ptype} - admitted at {patient.time_admitted}')
# sample duration of stay
if patient.ptype == 'A':
stay_duration = lognormal(mean=normal_mean_a, sigma=normal_stdev_a)
elif patient.ptype == 'B':
stay_duration = lognormal(mean=normal_mean_b, sigma=normal_stdev_b)
else:
stay_duration = lognormal(mean=normal_mean_c, sigma=normal_stdev_c)
yield env.timeout(stay_duration)
patient.time_discharged = env.now
print(f'Patient {patient.pid} - {patient.ptype} - discharged at {patient.time_discharged}')
env = simpy.Environment()
beds = simpy.Resource(env, capacity=total_beds)
env.process(generate_arrivals(env, 'A'))
env.process(generate_arrivals(env, 'B'))
env.process(generate_arrivals(env, 'C'))
env.run(until=run_length)
The problem with this approach was that there is a patient of each category created at the beginning of their respective process (it always creates patients 1, 2 & 3 of each category at time = 0), which doesn't reflect a real-world scenario:
Patient 1 - A - arrived at 0
Patient 2 - B - arrived at 0
Patient 3 - C - arrived at 0
Patient 1 - A - admitted at 0
Patient 2 - B - admitted at 0
Patient 3 - C - admitted at 0
Patient 1 - A - discharged at 0.04165029350880402
Patient 4 - A - arrived at 0.6503311494436321
Patient 4 - A - admitted at 0.6503311494436321
Patient 4 - A - discharged at 0.6626671650563922
Patient 5 - B - arrived at 0.6868621026906724
.
.
etc
So far, I've been toying with the idea of using:
numpy.random.multinomial() with equal probabilities.import simpy
from numpy.random import exponential, lognormal, multinomial
counter = 0
total_beds = 5
run_length = 10
ptypes = ['A', 'B', 'C']
probvals = [1.0 / len(ptypes) for _ in ptypes]
class Patient:
def __init__(self, pid, ptype):
self.pid = pid
self.ptype = ptype
self.time_arrived = None
self.time_admitted = None
self.time_discharged = None
class TypeA(Patient):
def __init__(self, pid, ptype):
super().__init__(pid, ptype)
self.mean_iat = 1.0
self.normal_mean = 4.0
self.normal_stdev = 7.0
class TypeB(Patient):
def __init__(self, pid, ptype):
super().__init__(pid, ptype)
self.mean_iat = 2.0
self.normal_mean = 5.0
self.normal_stdev = 8.0
class TypeC(Patient):
def __init__(self, pid, ptype):
super().__init__(pid, ptype)
self.mean_iat = 3.0
self.normal_mean = 6.0
self.normal_stdev = 9.0
def generate_arrivals(env):
while True:
global counter
counter += 1
ptype = ptypes[multinomial(n=1, pvals=probvals).argmax()]
# From here is where I'm stuck
# Can't understand how to drive the correct object instantiation
# to replace the if-elif-else from the previous example
patient = Patient(counter, ptype)
# Not sure of the logic from here on out
env.process(perform_treatment(env, patient))
interarival_time = exponential(scale=patient.mean_iat)
yield env.timeout(interarival_time)
def perform_treatment(env, patient):
patient.time_arrived = env.now
print(f'Patient {patient.pid} - {patient.ptype} - arrived at {patient.time_arrived}')
with beds.request() as req:
yield req
patient.time_admitted = env.now
print(f'Patient {patient.pid} - {patient.ptype} - admitted at {patient.time_admitted}')
# Again, how do I replace the if-elif-else from the previous example
stay_duration = lognormal(mean=patient.normal_mean, sigma=patient.normal_stdev)
yield env.timeout(stay_duration)
patient.time_discharged = env.now
print(f'Patient {patient.pid} - {patient.ptype} - discharged at {patient.time_discharged}')
env = simpy.Environment()
beds = simpy.Resource(env, capacity=total_beds)
env.process(generate_arrivals(env))
env.run(until=run_length)
However, I can't seem to understand how to go about it using Simpy. Any inputs on my approach, or alternative ideas would be greatly appreciated!
Note: I'm aware that using numpy.random.multinomial() results in only a single kind of patient being generated at a time, but that's probably a different question altogether.
in your first code, can you move your gen patient code, can you move your timeout to the top of the loop. Of course this means you would not get any patients at time 0
every process you start before the env.run starts at time 0, and everything that happens before the first yield happens at time 0. calling a env.timeout does not affect the other process. I think of each process and its own thread. So when your program calls env.run. It grabs one of your process and creates a patient and then waits at the yield env. timeout(), it then grabs the second process, still at time 0, and creates a patient and pauses at the yield. ect... ect... in short time does not move in a process until it gets to a yield. Does this help? – Michael 11 hours ago Delete when the first line in the process is a yield env.timeout() nothing happens until the timeout is done and the time has been updated. It has the effect of scheduling when the process should star