Collecting results from a single run

A tool like simpy allows you to collect your data flexibly using an approach that makes sense to you! Some options are:

1. Code an auditor / observer process. This process will periodically observe the state of the system. We can use this to collect information on current state at time t. For example, how many patients are queuing and how many have a call in progress between by time of day.

2. Store process metrics during a run and perform calculations at the end of a run. For example, if you want to calculate mean patient waiting time then store each patient waiting time in a list and calculate the mean at the end of the run.

3. Conduct and audit or calculate running statistics as the simulation executes an event. For example, as a patient completes a call we can calculate a running mean of waiting times and a running total of the operators are taking calls. The latter measure can then be used to calculate server utilisation. You could also use this approach to audit queue length where the queue length is recorded each time request for a resource is made (and/or when a resource is released).

1. Imports

import simpy
import numpy as np
import itertools

2. Calculating mean waiting time

The second strategy to results collection is to store either a reference to a quantitative value (e.g. waiting time) during the run. Once the run is complete you will need to include a procedure for computing the metric of interest. An advantage of this strategy is that it is very simple, captures all data, and has minimal computational overhead during a model run! A potential disadvantage is that for complex simulation you may end up storing a large amount of data in memory. In these circumstances, it may be worth exploring event driven strategies to reduce memory requirements.

In our example, we will store each patient’s waiting time in a python list. At the end of the run we will loop through these references and calculate mean waiting time and operator utilisation.

To do this I’m going to declare a list with notebook level scope. The service function will then append a waiting_time for a caller to the list each time the caller enters service.

2.1 Notebook level variables for results collection.

We will create a python dictionary called results to store result collection variables. This means that it is simple to add new variables in at a later date.

The dictionary has notebook level scope. This means that any functions or class in the notebook can access and/or append to the list access via the key waiting_times.

results = {}
results['waiting_times'] = []

2.2 A helper function

We will create a helper function called trace that wraps print. We can set a variable called TRACE that switches printing patient level results on and off.

def trace(msg):
    '''
    Turing printing of events on and off.
    
    Params:
    -------
    msg: str
        string to print to screen.
    '''
    if TRACE:
        print(msg)

2.3 Service and arrival functions

The only modification we need to make is to the service function. We will add in a line of code to record the waiting_time of the caller as they enter service.

def service(identifier, operators, env):
    '''
    simulates the service process for a call operator

    1. request and wait for a call operator
    2. phone triage (triangular)
    3. exit system
    
    Params:
    ------
    
    identifier: int 
        A unique identifer for this caller
        
    operators: simpy.Resource
        The pool of call operators that answer calls
        These are shared across resources.
        
    env: simpy.Environment
        The current environent the simulation is running in
        We use this to pause and restart the process after a delay.
    
    '''
    # record the time that call entered the queue
    start_wait = env.now

    # request an operator
    with operators.request() as req:
        yield req

        # record the waiting time for call to be answered
        waiting_time = env.now - start_wait
        
        # MODIFICATION - store the waiting time.
        results['waiting_times'].append(waiting_time)

        trace(f'operator answered call {identifier} at ' \
              + f'{env.now:.3f}')

        # sample call duration.
        call_duration = np.random.triangular(left=5.0, mode=7.0,
                                             right=10.0)
        
        # schedule process to begin again after call_duration
        yield env.timeout(call_duration)

        # print out information for patient.
        trace(f'call {identifier} ended {env.now:.3f}; ' \
              + f'waiting time was {waiting_time:.3f}')

def arrivals_generator(env, operators):
    '''
    IAT is exponentially distributed

    Parameters:
    ------
    env: simpy.Environment
        The simpy environment for the simulation

    operators: simpy.Resource
        the pool of call operators.
    '''

    # use itertools as it provides an infinite loop 
    # with a counter variable that we can use for unique Ids
    for caller_count in itertools.count(start=1):

        # 100 calls per hour (units = hours). 
        # Time between calls is 1/100
        inter_arrival_time = np.random.exponential(60/100)
        yield env.timeout(inter_arrival_time)

        trace(f'call arrives at: {env.now:.3f}')

        # create a new simpy process for this caller.
        # we pass in the caller id, the operator resources, and env.
        env.process(service(caller_count, operators, env))

# model parameters
RUN_LENGTH = 1000
N_OPERATORS = 13

# MODIFICATION - turn off caller level results.
TRACE = False

# create simpy environment and operator resources
env = simpy.Environment()
operators = simpy.Resource(env, capacity=N_OPERATORS)

env.process(arrivals_generator(env, operators))
env.run(until=RUN_LENGTH)
print(f'end of run. simulation clock time = {env.now}')

# MODIFICATION calculate results on notebook level variables.
mean_wt = np.mean(results['waiting_times'])
print(f'Mean waiting time was {mean_wt:.2f}')

end of run. simulation clock time = 1000
Mean waiting time was 4.73