Emulate Simpy's Machine Shop example

[galenseilis]

Introduction

I am doing a compare and contrast of the features available in Simpy and Ciw. The Simpy documentation has this machine shop example. I am having some difficulty emulating it using Ciw, and I would appreciate troubleshooting assistance.

Description

Some of my approach is supposed to work by implicit representation of certain aspects of the machine shop example.

Operational vs Broken Represented as Nodes

I have attempted to represent the state of a machine being broken by being in either Node 1 (operational) or Node 2 (broken). Going from operational to broken is achieved by a machine at Node 1 being completely served, sending it off to Node 2.

Infinite Servers to Remove Waiting in Line to Wait to be Broken

Machines do not wait in line to be broken, so I figured that the number of servers at Node 1 being infinite would mean that any machine at Node 1 begins service as soon as that part of the state is updated. Thus being served at Node 1 is being used to represent the waiting of a machine to break.

Initial Conditions & Arrivals

All of the machines begin in an operational state, so they should begin at Node 1. I wrote a time-dependent distribution InitialBatch which is supposed add the 10 machines to Node 1 at the start of the simulation (t=0).

There are no arrivals except for the initial batch of machines, so I have set the arrival distributions to None for both nodes.

Calculating Parts Produced

I have not implemented it yet, but I believe that the number of parts can be calculated after the simulation using the records. Specifically, the machine shops produce an average of one part every 10 minutes with a standard deviation of two minutes. The original Simpy implementation actually breaks if you remove their seed and do repeated trials because they use a normal distribution which can (and does occasionally) sample negative values. But you cannot produce a part in negative time, and it actually leads to an error. I raised this issue on the Simpy GiLab.

Another option might be to look at the state trackers, which I just have not read enough about yet.

Repair Process

When a machine breaks it may have to wait its turn to be repaired by the repairman, so I have left service discipline to the default (which at the time of writing I believe is FIFO). Every machine repair takes 30 minutes, so I assigned a deterministic distribution with value=30 to the service time.

Repairman's Other Jobs

The repairman's other jobs are immaterial to the behaviour of the machines and the number of parts produced. Machines prempt the repairman's other jobs. But if we wanted to compute the amount of jobs the repairman completed, we could again compute this from the records about where the machines were. Whenever none of the machines are broken then the repairman isn't repairing and is instead doing other jobs.

Routing

Machines at Node 1 always travel to Node 2, and machines at Node 2 always travel to Node 1. I tried to deal with this a couple of different ways. The first way was to have the stochastic matrix [[0,1],[1,0]].

The second approach, implemented below, was to use a custom routing function. Since all machines should start at Node 1, I only pass the routing function in the first position and use the ciw.no_routing placeholder as recommended in the documentation.

Attempt

I have attempted an implementation, which I have provided below.

import ciw

def routing_function(ind):
    return [2,1]

class InitialBatch(ciw.dists.Distribution):
    '''Initial batch of arrivals.'''

    def __init__(self, initial_arrivals):
        '''
        initial_arrivals (int): Number of initial arrivals.
        '''
        super().__init__()

        if isinstance(initial_arrivals, int):
            self.initial_arrivals = initial_arrivals
        else:
            raise TypeError('An integer was not given for n.')

    def sample(self, t, ind=None):
        '''Sample the initial arrivals.'''

        if t == 0:
            return self.initial_arrivals

        return 0

N = ciw.create_network(
    arrival_distributions=[None, None],
    service_distributions=[ciw.dists.Exponential(rate=1/300),
                           ciw.dists.Deterministic(value=30.0)],
    routing=[routing_function, ciw.no_routing],
    batching_distributions=[InitialBatch(10), None],
    number_of_servers=[float('inf'), 1]
)

sim = ciw.Simulation(N)
sim.simulate_until_max_time(4 * 7 * 24 * 60)

recs = sim.get_all_records()
print(recs)

Symptoms

• sim.current_time is Inf at the end of simulation.
• recs are an empty list.
• When I use a progress bar it indicates that the simulation does in fact run.
• At the end of the simulation there are no individuals in any of the nodes.

[galenseilis]

I think I misunderstood how batches work, especially time-dependent batches. I thought it was setting arrivals of some integer number of individuals at a time determined in time-dependent batch arrival function. Rather, I think I've learned that it is only controlling information about the number of arrivals, not the when they arrive.

I have attempted a workaround to put individuals on a queue, which I have posted here. Let me know if the workaround actually works. If it does, please feel to close this issue. If it doesn't, please explain how the attempted workaround does not work.

Thanks! 😊

[geraintpalmer]

Hi @galenseilis, thanks for the interest in the library.

I think you might have overcomplicated the initial batch. Here is my solution:

import ciw

N = ciw.create_network(
    arrival_distributions=[ciw.dists.Sequential([0, float('inf')]), None],
    service_distributions=[ciw.dists.Exponential(rate=1/300),
                           ciw.dists.Deterministic(value=30.0)],
    routing=[[0, 1], [1, 0]],
    batching_distributions=[ciw.dists.Deterministic(10), None],
    number_of_servers=[float('inf'), 1]
)

sim = ciw.Simulation(N)
max_sim_time = 4 * 7 * 24 * 60
sim.simulate_until_max_time(max_sim_time)

recs = sim.get_all_records()

Here I say that the arrivals at Node 1 follow a Sequential distribution, which just lists all inter-arrival times sequentially: 0 is immediately after the simulation begins, then infinity fo no arrivals ever after. Batching at Node 1 is always 10 (but there will only ever be 1 arrival time, so it results that this will only ever happen once).

Note that a the routing_function you defined means that arriving customers are first routed to Node 2, then back to Node 1, and then will leave the system. But using the transition matrix [[0, 1], [1, 0]] means they will continue to carry on breaking down and never leave the system.

Then to get the total uptime for each machine:

def get_uptime_machine_i(i, recs):
    # All uptimes before the currently interrupted session
    completed_uptimes = recs[(recs['id_number'] == i) & recs['node']==1]['service_time'].sum()

    # Get the current uptime interrupted by ending the simulation
    inds = sim.get_all_individuals()
    ind_i = [ind for ind in inds if ind.id_number == i][0]
    if ind_i.node == 1:
        current_uptime = max_sim_time - ind_i.arrival_date
    else:
        current_uptime = 0

    # Return the sum
    return completed_uptimes + current_uptime

and to get the total time the repairman was free (not repairing machines):

def get_repairman_freetime(recs, max_sim_time):
    # Get all completed repairs
    completed_services = recs[recs['node'] == 2]['service_time'].sum()

    # Get the times of the repairs interrupted by ending the simulation
    current_services = sum([max_sim_time - ind.arrival_date for ind in sim.get_all_individuals() if ind.node == 2])

    # Get the total service time
    total_services = completed_services + current_services

    # Get the free time of the repairman
    return max_sim_time - total_services

Let me know if this is the desired behaviour.

[geraintpalmer]

I have attempted a workaround to put individuals on a queue, which I have posted here. Let me know if the workaround actually works. If it does, please feel to close this issue. If it doesn't, please explain how the attempted workaround does not work.

Thanks for this!

I think we need more than just attaching the simulation object to the new individuals, as we need to set a few different attributes to get the desired effect.

This was my attempt in a previous project:

import math
import ciw

def begin_service_if_possible_accept(node, next_individual):
    """
    Begins the service of the next individual (at acceptance point):
        - Sets the arrival date as the current time
        - If there is a free server or there are infinite servers:
            - Attach the server to the individual (only when servers are not infinite)
        - Get service start time, service time, service end time
        - Update the server's end date (only when servers are not infinite)
    """
    free_server = node.find_free_server(next_individual)
    if free_server is None and math.isinf(node.c) is False:
        node.decide_preempt(next_individual)
    if free_server is not None or math.isinf(node.c):
        if math.isinf(node.c) is False:
            node.attach_server(free_server, next_individual)
        next_individual.service_start_date = 0
        next_individual.service_time = node.get_service_time(next_individual)
        next_individual.service_end_date = next_individual.service_time
        if not math.isinf(node.c):
            free_server.next_end_service_date = next_individual.service_end_date

# from the ciw source code; needs to be adapted to accommodate the pre-existing wait list
def accept(individual_id, individual_class, arrival_date, node, simulation):
    """
    Accepts a new customer to the queue:
      - record all other information at arrival point
      - update state tracker
    """
    simulation.current_time = arrival_date
    next_individual = ciw.Individual(id_number=individual_id, customer_class=individual_class, priority_class=simulation.network.priority_class_mapping[individual_class], simulation=simulation)
    next_individual.node = node.id_number
    next_individual.queue_size_at_arrival = 'Unknown'
    node.individuals[next_individual.priority_class].append(next_individual)
    node.number_of_individuals += 1
    node.simulation.statetracker.change_state_accept(node, next_individual)
    next_individual.arrival_date = arrival_date
    begin_service_if_possible_accept(node, next_individual)

    simulation.nodes[0].number_of_individuals += 1
    simulation.nodes[0].number_of_individuals_per_class[next_individual.customer_class] += 1
    simulation.nodes[0].number_accepted_individuals += 1
    simulation.nodes[0].number_accepted_individuals_per_class[next_individual.customer_class] += 1

def create_existing_customers_from_list(backlog, Q):
    customer_count = 1
    for row in backlog:
        customer_id = row[0]
        customer_class = row[1]
        customer_arrival_date = row[2]
        customer_node = row[3]
        accept(customer_id, customer_class, customer_arrival_date, Q.nodes[customer_node], Q)
        customer_count += 1

and then to manually insert customers like so:

N = ciw.create_network(
    arrival_distributions=[ciw.dists.Exponential(5), ciw.dists.Exponential(5)],
    service_distributions=[ciw.dists.Exponential(6), ciw.dists.Exponential(6)],
    number_of_servers=[2, 2],
    routing=[[0, 0.5], [0, 0]]
)

Q = ciw.Simulation(N)

backlog = [
    [1, 'Customer', -5.5, 1],
    [2, 'Customer', -4.3, 2],
    [3, 'Customer', -2.2, 1],
    [4, 'Customer', -0.7, 1]
]

create_existing_customers_from_list(backlog, Q)

But I think having some Ciw functions to do this for us would be a great idea.

[galenseilis]

I have attempted a workaround to put individuals on a queue, which I have posted here. Let me know if the workaround actually works. If it does, please feel to close this issue. If it doesn't, please explain how the attempted workaround does not work.

Thanks for this!

I think we need more than just attaching the simulation object to the new individuals, as we need to set a few different attributes to get the desired effect.

This was my attempt in a previous project:

import math
import ciw

def begin_service_if_possible_accept(node, next_individual):
    """
    Begins the service of the next individual (at acceptance point):
        - Sets the arrival date as the current time
        - If there is a free server or there are infinite servers:
            - Attach the server to the individual (only when servers are not infinite)
        - Get service start time, service time, service end time
        - Update the server's end date (only when servers are not infinite)
    """
    free_server = node.find_free_server(next_individual)
    if free_server is None and math.isinf(node.c) is False:
        node.decide_preempt(next_individual)
    if free_server is not None or math.isinf(node.c):
        if math.isinf(node.c) is False:
            node.attach_server(free_server, next_individual)
        next_individual.service_start_date = 0
        next_individual.service_time = node.get_service_time(next_individual)
        next_individual.service_end_date = next_individual.service_time
        if not math.isinf(node.c):
            free_server.next_end_service_date = next_individual.service_end_date

# from the ciw source code; needs to be adapted to accommodate the pre-existing wait list
def accept(individual_id, individual_class, arrival_date, node, simulation):
    """
    Accepts a new customer to the queue:
      - record all other information at arrival point
      - update state tracker
    """
    simulation.current_time = arrival_date
    next_individual = ciw.Individual(id_number=individual_id, customer_class=individual_class, priority_class=simulation.network.priority_class_mapping[individual_class], simulation=simulation)
    next_individual.node = node.id_number
    next_individual.queue_size_at_arrival = 'Unknown'
    node.individuals[next_individual.priority_class].append(next_individual)
    node.number_of_individuals += 1
    node.simulation.statetracker.change_state_accept(node, next_individual)
    next_individual.arrival_date = arrival_date
    begin_service_if_possible_accept(node, next_individual)

    simulation.nodes[0].number_of_individuals += 1
    simulation.nodes[0].number_of_individuals_per_class[next_individual.customer_class] += 1
    simulation.nodes[0].number_accepted_individuals += 1
    simulation.nodes[0].number_accepted_individuals_per_class[next_individual.customer_class] += 1

def create_existing_customers_from_list(backlog, Q):
    customer_count = 1
    for row in backlog:
        customer_id = row[0]
        customer_class = row[1]
        customer_arrival_date = row[2]
        customer_node = row[3]
        accept(customer_id, customer_class, customer_arrival_date, Q.nodes[customer_node], Q)
        customer_count += 1

and then to manually insert customers like so:

N = ciw.create_network(
    arrival_distributions=[ciw.dists.Exponential(5), ciw.dists.Exponential(5)],
    service_distributions=[ciw.dists.Exponential(6), ciw.dists.Exponential(6)],
    number_of_servers=[2, 2],
    routing=[[0, 0.5], [0, 0]]
)

Q = ciw.Simulation(N)

backlog = [
    [1, 'Customer', -5.5, 1],
    [2, 'Customer', -4.3, 2],
    [3, 'Customer', -2.2, 1],
    [4, 'Customer', -0.7, 1]
]

create_existing_customers_from_list(backlog, Q)

But I think having some Ciw functions to do this for us would be a great idea.

Cool! Thank you for attempting this. I'll give it a try.

I agree, having some supported way to start with a desired initial state would be quite beneficial. It can be especially useful for planning to compare different scenarios assuming that the model starts in something similar to the actual current state of a waitlist.

[galenseilis]

I have attempted a workaround to put individuals on a queue, which I have posted here. Let me know if the workaround actually works. If it does, please feel to close this issue. If it doesn't, please explain how the attempted workaround does not work.

Thanks for this!

I think we need more than just attaching the simulation object to the new individuals, as we need to set a few different attributes to get the desired effect.

This was my attempt in a previous project:

import math
import ciw

def begin_service_if_possible_accept(node, next_individual):
    """
    Begins the service of the next individual (at acceptance point):
        - Sets the arrival date as the current time
        - If there is a free server or there are infinite servers:
            - Attach the server to the individual (only when servers are not infinite)
        - Get service start time, service time, service end time
        - Update the server's end date (only when servers are not infinite)
    """
    free_server = node.find_free_server(next_individual)
    if free_server is None and math.isinf(node.c) is False:
        node.decide_preempt(next_individual)
    if free_server is not None or math.isinf(node.c):
        if math.isinf(node.c) is False:
            node.attach_server(free_server, next_individual)
        next_individual.service_start_date = 0
        next_individual.service_time = node.get_service_time(next_individual)
        next_individual.service_end_date = next_individual.service_time
        if not math.isinf(node.c):
            free_server.next_end_service_date = next_individual.service_end_date

# from the ciw source code; needs to be adapted to accommodate the pre-existing wait list
def accept(individual_id, individual_class, arrival_date, node, simulation):
    """
    Accepts a new customer to the queue:
      - record all other information at arrival point
      - update state tracker
    """
    simulation.current_time = arrival_date
    next_individual = ciw.Individual(id_number=individual_id, customer_class=individual_class, priority_class=simulation.network.priority_class_mapping[individual_class], simulation=simulation)
    next_individual.node = node.id_number
    next_individual.queue_size_at_arrival = 'Unknown'
    node.individuals[next_individual.priority_class].append(next_individual)
    node.number_of_individuals += 1
    node.simulation.statetracker.change_state_accept(node, next_individual)
    next_individual.arrival_date = arrival_date
    begin_service_if_possible_accept(node, next_individual)

    simulation.nodes[0].number_of_individuals += 1
    simulation.nodes[0].number_of_individuals_per_class[next_individual.customer_class] += 1
    simulation.nodes[0].number_accepted_individuals += 1
    simulation.nodes[0].number_accepted_individuals_per_class[next_individual.customer_class] += 1

def create_existing_customers_from_list(backlog, Q):
    customer_count = 1
    for row in backlog:
        customer_id = row[0]
        customer_class = row[1]
        customer_arrival_date = row[2]
        customer_node = row[3]
        accept(customer_id, customer_class, customer_arrival_date, Q.nodes[customer_node], Q)
        customer_count += 1

and then to manually insert customers like so:

N = ciw.create_network(
    arrival_distributions=[ciw.dists.Exponential(5), ciw.dists.Exponential(5)],
    service_distributions=[ciw.dists.Exponential(6), ciw.dists.Exponential(6)],
    number_of_servers=[2, 2],
    routing=[[0, 0.5], [0, 0]]
)

Q = ciw.Simulation(N)

backlog = [
    [1, 'Customer', -5.5, 1],
    [2, 'Customer', -4.3, 2],
    [3, 'Customer', -2.2, 1],
    [4, 'Customer', -0.7, 1]
]

create_existing_customers_from_list(backlog, Q)

But I think having some Ciw functions to do this for us would be a great idea.

Yes, I ran this example and it seems to work! I agree that being able to do this with Ciw would be desirable.