diff --git a/.gitignore b/.gitignore index 0342d56..f0ab546 100644 --- a/.gitignore +++ b/.gitignore @@ -1,31 +1,20 @@ **/.idea/ **/__pycache__/ **/.ipynb_checkpoints -*.pyc - -*.vscode/ - -mow_demo2.ipynb - -simpy_test.ipynb +validation/ +logs/ +simenv/ -maintsim/test.ipynb +setup.py -*.csv *.ipynb -!demo.ipynb +*.csv +*.html +*.pkl +*.pyc +*.txt +*.vscode/ +*.xlsx -maintsim/.ipynb_checkpoints/test-checkpoint.ipynb -maintsim.py -simpy_test.ipynb -maintsim_example.html -maintsim_example.ipynb -mow_demo.ipynb -mow_demo2.ipynb -pm_reliability.ipynb -pmow_demo.ipynb -maintsim/concurrent downtime demo.ipynb -maintsim/degradation.py -maintsim/Simpy demo.ipynb -maintsim/system viz.ipynb -maintsim/demo.html \ No newline at end of file +!requirements.txt +!LICENSE.txt \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt index 944c0a3..6b59a08 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,6 +1,10 @@ MIT License +<<<<<<< HEAD Copyright (c) 2019 Michael Hoffman +======= +Copyright (c) 2020 Michael Hoffman +>>>>>>> dev Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index 2005445..afe4349 100644 --- a/README.md +++ b/README.md @@ -1,159 +1,197 @@ -# maintsim +# Simantha +<<<<<<< HEAD `maintsim` can be used to model a discrete manufacturing system where components degrade over time and receive maintenance. Users can define the configuration and parameters of the system, as well as the maintenance policy to be carried out. It is built on the `SimPy` discrete-event simulation package. ## Installing maintsim `pip install maintsim` +======= +Simantha uses discrete event simulation to model the behavior of discrete manufacturing systems, in particular the production and maintenance functions of a system. +>>>>>>> dev ## Using this package ### Requirements -maintsim relies on the following packages in addition to Python 3.6+: - -- [SimPy](https://simpy.readthedocs.io/en/latest/) version 3.0.11 -- [pandas](https://pandas.pydata.org/) version >= 0.23.4 -- [SciPy](https://www.scipy.org) version >= 1.1.0 (if specifying random repair times) +Simantha requires Python ≥ 3.6 and [SciPy](https://www.scipy.org/) ≥ 1.5.2 for running tests. ### Setting up a manufacturing system -The workflow begins by creating a `System` object that is defined by the following parameters: +A system consists of various assets that are linked together by a specified routing. The available assets are described below. -#### Configuration parameters -- `process_times` - a list of process times for each machine in a serial line. -- `interarrival_time` - the time between part arrivals to the first machine in the system. The default is 1, which ensures the first machine is never starved. -- `buffer_sizes` - a list of buffer sizes for each machine, or an integer value for the buffer size of all machines. Default is 1. For n machines there are n-1 buffers. +All asset classes take an optional `name` parameter which is useful for identifying unique objects in the simulation trace. -#### Failure parameters -- `failure_mode` - currently either `'degradation'` or `None`. Each machine is subject to the same mode of degradation. By default machines do not degrade. - - `'degradation'` - machine degradation occurs according to a discrete-state Markovian process. Currently only Markovian degradation is supported. Requires specification of `failure_params`, described in the following subsection. -- `planned_failures` - a list of planned failures to occur during the simulation time. Each of the form `(location, time, duration)`. Planned failures do not adhere to maintenance capacity constraints and have not been thoroughly tested in tandem with random failures. +#### Source -##### Markovian degradation parameters +`simantha.Source` -There are several options for specifying the mode of degradation using the `failure_params` argument which should be passed as a dictionary. -- Constant degradation rate - passed as the value of the `degradation rate` key either as a single float or list of floats for each machine. The value is the probability of degrading by one unit at each time step and so should be between 0 and 1. Creates an upper bidiangonal degradation transition matrix. -- Failed state - if the degradation rate is specified, the maximum health (failed) state for each machine can be passed as the value of the `failed state` key. The default failed state is 10. -- Complete degradation transition matrix - can be specified as a single `numpy` array (in which case each machine will be subject to the same degradation profile) or a list of arrays for each machine as the value of the `degradation transition` key. +A source introduces parts to the system. Each source should be directly upstream of a machine. -If `failure_mode = 'degradation'` is passed to a `System` object, either the `degradation rate` or `degradation transition` must be defined. Degradation profiles beyond a constant degradation rate have not yet been thoroughly tested. +Parameters +- `interarrival_time` - time between arrivals generated by the source. If an arrival occurs but there is no available space for the part, the new arrival will be seized as soon as space becomes available. By default, the interarrival time is zero, meaning machines downstream of the source will never be starved. -#### Maintenance parameters -- `maintenance_policy` - currently either `'CM'` or `'CBM'`. - - `'CM'` - "corrective maintenance", the default policy, machines will only be repaired upon complete failure, as determined by the mode of degradation. - - `'CBM'` - "condition-based maintenance", preventive maintenance is performed once a machine's condition reached a prescribed threshold. -- `maintenance_params` - the parameters that define the specified maintenance policy. For `CBM`, a list of thresholds at which to schedule maintenance. - - Currently each machine has 11 health states, with 0 being perfect health and 10 being the failed state. The maintenance threshold should be in this range. -- `repair_params` - a dictionary of `scipy.stats` frozen discrete distributions of time to repair based on repair type. - - For example, `repair_params = {'CM': stats.randint(10, 20), 'CBM': stats.randint(20, 40)}`. -- `maintenance_capacity` - the maximum number of maintenance jobs that can be executed simultaneously. Currently if the number of simultaneously scheduled maintenance jobs exceeds the capacity they will be handled in a first in, first out (FIFO) manner. -- `maintenance_costs` - a dictionary of the cost of each type of maintenance job by type. +#### Machine -#### System state parameters +`simantha.Machine` -These parameters can be set to initialize the system state before simulation. By default the system begins with empty buffers and machines in perfect health. -- `initial_remaining_process` - a list of remaining processing times for each machine. By default this is equal to the machine's total processing time when it does not have a part. -- `initial_buffer` - a list of initial levels for each buffer. -- `initial_health` - a list of initial health states for each machine. -Valid settings for the initial system state are not currently verified automatically. +While a machine operates it will continuously take parts from an upstream container (if the container is not empty), process the part, and place the part in a downstream container (if the container is not full). Machines can be subject to periodic degradation and failure, at which point they will require maintenance action before they can be restored. +Currently, machines may only follow a Markovian degradation process[^1]. Under this degradation mode, a degradation transition matrix is specified as a parameter of a machine. State transitions of this Markov process represent changes in the degradation level of the machine. In general, it is assumed that there is one absorbing state that represents machine failure. -### Creating a custom maintenance scheduler +[^1]: Chan, G. K., & Asgarpoor, S. (2006). Optimum maintenance policy with Markov processes. Electric power systems research, 76(6-7), 452-456. https://doi.org/10.1016/j.epsr.2005.09.010 -The `System` object can take an additional `scheduler` object that will determine how maintenance jobs are scheduled when the number of machines due for maintenance exceeds the maintenance capacity. By default a system will use a FIFO scheduler as defined by the `Scheduler` class. A custom scheduler can also be created that inherits from the `maintsim.Scheduler` class. This new class should include a `choose_next` method that accepts the current maintenance queue as an argument. This method should then return a list of machine objects that are to be assigned maintenance. The `choose_next` method is executed each time a maintenance resource is release from a job. +Parameters +- `cycle_time` - the duration of time a machine must process a part before it can be placed at the next station. +- `degradation_matrix` - the degradation transition matrix for Markovian degradation. For an `n`x`n` transition matrix, degradation levels are integer-labeled from `0`, indicating perfect health, to `n-1`, indicating the machine is failed. +- `cbm_threshold` - the condition-based maintenance threshold for creating a request for a maintenance resource. For an `n`x`n` degradation transition matrix, this threshold should be in the interval `[0, n-1]`. +- `pm_distribution` - random distribution of the duration of any preventive maintenance job performed on a machine. Any maintenenance that is performed on a machine that is not in the failed state is considered preventive. +- `cm_distribution` - random distribution of the duration of any corrective maintenance job performed on a machine. A machine in a failed state can only be serviced by corrective maintenance. -An example of a custom scheduler that uses Monte Carlo tree search (as implemented by the [MCTS](https://github.com/pbsinclair42/MCTS) package) is shown below: +Currently, only a uniform random distrubtion is supported and should take the form `{'uniform': [a, b]}` when passed as an argument to a `Machine`. -```python -import maintsim -import mcts - -class MCTSScheduler(maintsim.Scheduler): - ''' - Resolves maintenance scheduling conflicts using Monte Carlo tree search. - ''' - def __init__(self, time_limit=None, iteration_limit=None, **kwds): - ''' - Must specify either a time limit (in seconds) or iteration limit for the - MCTS. - ''' - super().__init__(**kwds) - self.limit = {} - if time_limit and iteration_limit: - print('Error: cannot specify time and iteration limit.') - elif time_limit: - self.limit['timeLimit'] = time_limit * 1000 - else: - self.limit['iterationLimit'] = iteration_limit - - def choose_next(self, queue): - # formulate and solve MCTS - mcts_schedule = mcts.mcts(**self.limit) - best_action = mcts_schedule.search(initialState=MaintenanceState(self.system)) - next_machine = [self.system.machines[best_action-1]] - - return next_machine -``` +#### Buffer + +`simantha.Buffer` + +Buffers serve as passive containers for parts that are not being processed by a machine. Generally, each machine should retrieve parts from one buffer (or a source) and place parts in one buffer (or a sink). If machines are assigned more than one immediate upstream or downstream buffer unexpected transfer behavior can occur. + +Parameters +- `capacity` - the maximum number of parts that may occupy a buffer at one time. + +#### Sink + +`simantha.Sink` + +A sink collects finished parts that exit the system. Each part that is placed in a sink is counted as a part that has been successfully produced by the system. The total level of all sinks is considered to be the overall system production. + +Parameters +- `initial_level` - the number of finished parts in a sink when the system is initialized. Sinks are intialized as empty by default. + +#### Maintainer + +`simantha.Maintainer` + +A maintainer is responsible for repairing machines that have either failed or requested maintenance. If the number of machines requesting maintenance exceeds the maintainer's available capacity, the default behavior is to maintain the machine that requested maintenance the earilest (a first-come, first-serve policy). This behvaior can be modified by overriding the `choose_maintenance_action` method of the `simantha.Maintainer` class. This method takes the current `queue` in the form of a list of `simantha.Machine` instances requesting maintenance as an argument and should return the instance that is assigned maintenance. + +Parameters +- `capacity` - the maximum number of machines that can be maintained simultaneously. + +#### Defining object routing + +Once the objects are created the routing(s) through the system need to be defined using each asset's `define_routing` method. This method takes two arguments, `upstream` and `downstream`, as lists of other objects in the system directly adjacent to the object calling the method. -### Simulating the system -When the system is instantiated, it will initialize by creating the necessary objects including the SimPy `Environment`, the maintenance resource, machines, and buffers. The simulation can be run by calling the `simulate` method of the `System` object with the following parameters: +#### System -- `title` - the tile of the simulation, used for naming any files that are saved. -- `warmup_time` - the time that the simulation will run before collecting data. Useful for ensuring the system is in steady state before observation. -- `sim_time` - the duration of time that simulation will run once the warmup is completed. Metrics will be determined based on the system performance during this time. -- `seed` - random seed for the simulation. A given seed should always produce the same results. -- `verbose` - boolean. `True` will print out a summary of the simulation run. `False` will suppress all printed output, which may be preferred if many replications are being run. +`simantha.System` +All of the objects are passed to a `System` object that can then be simulated. -#### Simulation data +Parameters +- `objects` - a list of objects including sources, machines, buffers, and sinks that make up the system. +- `maintainer` - an instance of `Maintainer` if one has been created. If one has not been created and machines are subject to degradation, a maintainer with the default first-in, first-out policy will be used. -Several data frames are created to record data of a simulation run and stored as attributes of the `System` object. +### Simulating a system -- `state_data` - remaining processing time of each machine and buffer levels at each time step. -- `production_data` - production volume (in units) and throughput (in units/time) of each machine at every time step. -- `machine_data` - status of each machine at every time step, including if machines are function and blocked or starved. -- `queue_data` - the number of machines waiting for maintenance at each time step. -- `maintenance_data` - the log of maintenance activities including the time at which the activity occurred, the type of activity (corrective, preventive, etc.), what the activity was (failure or repair), and the duration or time to failure. +A `simantha.System` object has two methods for simulating its behavior: +- `simulate` - simulates system behavior for the specified duration of time. Arguments of this method include + - `warm_up_time` - the duration of time the system operation is simulated before performance statistics are gathered. Since objects in the system are empty by default, specifying a warm up period will allow the system to reach steady state performance which is often of interest when conducting any analysis. + - `simulation_time` - the duration of time to simulate system behavior. Performance statistics will be gathered during this period. In total, the duration of simulated time is `warm_up_time` + `simulation_time`. + - `verbose` - `True` or `False`, indicating whether a summary of the simulation run should be displayed. `True` by default. + - `collect_data` - `True` or `False`, indicating whether or not data is collected for indiviudial objects in the system. If many simulation runs are conducted, setting this to `False` may improve performance. +- `iterate_simulation` - conduct multiple simulation runs of a system. Useful for estimating the average performance of a particular system whose behavior is random. This method uses Python's [multiprocessing](https://docs.python.org/3.8/library/multiprocessing.html) to call the `simulate` method in parallel. Arguments to this method are + - `replications` - the number of simulation runs to conduct. + - `warm_up_time` - used the same as in the `simulate` method and applied to each replication. + - `simulation_time` - used the same as in the `simulate` method. + - `verbose` - `True` or `False`, indicating whether or not a summary of all replications is displayed at completion. `True` by default. + - `jobs` - the number of worker processes that will be created by the multiprocessing module. By default, one worker is used which is the equivalent of running all replications in series. -#### Iterating the simulation +### Example usage -A system can be simulated several times using the `System.iterate_simulation` method. The arguments for this method are: +Below is a simple example demonstrating the production of a single machine. -- `replications` - number of times to run the simulation. -- `warmup_time` - the time that the simulation will run before collecting data. Performance statistics are only collected after this time has elapsed. -- `sim_time` - the duration of time each replication will be simulated. -- `objective` - the objective values that will be returned once all replications are complete. Options include: - - `production` - the production volume in units of the system. - - `ppl` - the permanent production loss of the system over the specified simulation time. - - `availability` - the overall availability of machines in the system as a percentage. -- `verbose` - `True` or `False`, determines whether or not summary statistics will be displayed once all replications are completed. +```python +>>> from simantha import Source, Machine, Sink, System +>>> +>>> source = Source() +>>> M1 = Machine('M1', cycle_time=1) +>>> sink = Sink() +>>> +>>> source.define_routing(downstream=[M1]) +>>> M1.define_routing(upstream=[source], downstream=[sink]) +>>> sink.define_routing(upstream=[M1]) +>>> +>>> system = System(objects=[source, M1, sink]) +>>> system.simulate(simulation_time=100) +Simulation finished in 0.00s +Parts produced: 100 +``` -## A simple example +As a slightly more complex example, the code below constructs a two-machine one-buffer line where each machine is subject to degradation. The degradation transition matrix will represent Bernoulli reliability with p = 0.01.[^2] Additionally, the maintainer capacity is set to 1, indicating that only one machine may be repaired at a time. + +[^2]: Li, J., & Meerkov, S. M. (2008). Production systems engineering. Springer Science & Business Media. https://doi.org/10.1007/978-0-387-75579-3 + +```python +>>> from simantha import Source, Machine, Buffer, Sink, System, Maintainer +>>> +>>> degradation_matrix = [[0.99, 0.01], [0., 1.]] +>>> cm_distribution = {'uniform': [10, 20]} +>>> +>>> source = Source() +>>> M1 = Machine( +... 'M1', +... cycle_time=1, +... degradation_matrix=degradation_matrix, +... cm_distribution=cm_distribution, +... ) +>>> B1 = Buffer(capacity=5) +>>> M2 = Machine( +... 'M2', +... cycle_time=1, +... degradation_matrix=degradation_matrix, +... cm_distribution=cm_distribution +... ) +>>> sink = Sink() +>>> +>>> source.define_routing(downstream=[M1]) +>>> M1.define_routing(upstream=[source], downstream=[B1]) +>>> B1.define_routing(upstream=[M1], downstream=[M2]) +>>> M2.define_routing(upstream=[M2], downstream=[sink]) +>>> sink.define_routing(upstream=[M2]) +>>> +>>> maintainer = Maintaner(capacity=1) +>>> +>>> system = System(objects=[source, M1, B1, M2, sink], maintainer=maintainer) +>>> system.simulate(simulation_time=1000) +Simulation finished in 0.05s +Parts produced: 748 +``` -Here is a minimum example for implmenting a CBM policy: +Lastly, the following example will show a system consisting of two machines in parallel. Both `M1` and `M2` retrieve parts from the same source and place finished parts in the same sink. ```python ->>> import maintsim ->>> from scipy import stats +>>> from simantha import Source, Machine, Sink, System +>>> +>>> source = Source() +>>> M1 = Machine('M1', cycle_time=1) +>>> M2 = Machine('M2', cycle_time=1) +>>> sink = Sink() >>> ->>> system = maintsim.System(process_times=[3, 5, 4], -... buffer_sizes=5, -... failure_mode='degradation', -... failure_params={'degradation rate':[0.25, 0.1, 0.2]}, -... maintenance_policy='CBM', -... maintenance_params={'CBM threshold': [8, 6, 7]}, -... repair_params={'CM': stats.randint(20,30), -... 'CBM': stats.randint(10,20)}, -... maintenance_capacity=1) ->>> system.simulate(warmup_time=100, sim_time=500) -Simulation complete in 0.89 s - - Units produced: 31 - System availability: 68.93% +>>> source.define_routing(downstream=[M1, M2]) +>>> for machine in [M1, M2]: +... machine.define_routing(upstream=[source], downstream=[sink]) +... +>>> sink.define_routing(upstream=[M1, M2]) +>>> +>>> system = System(objects=[source, M1, M2, sink]) +>>> system.simulate(simulation_time=100) +Simulation finished in 0.01s +Parts produced: 200 ``` +The elements of these examples can be used to more complex configurations of arbitrary structure. Additionally, the `Machine` class may represent other operations in a manufacturing system, such as an inspection station or a material handling process. + ## Planned features Key planned features include diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..ec47264 --- /dev/null +++ b/requirements.txt @@ -0,0 +1 @@ +scipy==1.5.2 \ No newline at end of file diff --git a/simantha/Asset.py b/simantha/Asset.py new file mode 100644 index 0000000..f1d531b --- /dev/null +++ b/simantha/Asset.py @@ -0,0 +1,31 @@ +import random + +class Asset: + """ + Parent class for assets in the system. All objects should extend this class, and + there should be no instances. + """ + def __init__(self, name, selection_priority=1): + self.name = name + # assets with higher priority will be selected over those with lower priority + # when competing for resources or space. + self.selection_priority = selection_priority + + self.upstream = [] + self.downstream = [] + + def initialize(self): + pass + + def get_candidate_givers(self, blocked=False): + """ + Returns a list of assets that can give a part to this asset from among the + upstream assets. + + If 'blocked=True' then only return assets that are currently blocked. + """ + candidates = [] + for candidate in self.upstream: + if blocked: + if candidate.blocked: + candidates.append(candidate) diff --git a/simantha/Buffer.py b/simantha/Buffer.py new file mode 100644 index 0000000..78cf6a4 --- /dev/null +++ b/simantha/Buffer.py @@ -0,0 +1,93 @@ +class Buffer: + def __init__(self, name='Buffer', capacity=float('inf'), initial_level=0): + self.name = name + self.capacity = capacity + self.initial_level = initial_level + self.level = initial_level + + self.env = None + + self.level_data = {'time': [0], 'level': [initial_level]} + + self.pending_requests = [] + + def initialize(self): + self.level = self.initial_level + + self.reserved_content = 0 + self.reserved_vacancy = 0 + + if self.env.collect_data: + self.level_data = {'time': [0], 'level': [self.initial_level]} + + def can_get_part(self): + return self.level + self.reserved_vacancy < self.capacity + + def reserve_content(self, quantity=1): + self.reserved_content += 1 + + def get(self, quantity=1): + if not self.is_empty(): + self.level -= quantity + self.reserved_content -= quantity + + if self.env.collect_data: + self.level_data['time'].append(self.env.now) + self.level_data['level'].append(self.level) + + else: + raise RuntimeError('Attempting to take more parts than available.') + + def reserve_vacancy(self, quantity=1): + self.reserved_vacancy += 1 + + def put(self, quantity=1): + if not self.is_full(): + self.level += quantity + self.reserved_vacancy -= 1 + + if self.env.collect_data: + self.level_data['time'].append(self.env.now) + self.level_data['level'].append(self.level) + + else: + raise RuntimeError('Attempting to put part in full buffer.') + + def is_empty(self): + #return self.level - self.reserved_content == 0 + return self.level == 0 + + def is_full(self): + #return self.level + self.reserved_vacancy == self.capacity + return self.level == self.capacity + + def can_give(self): + return self.level - self.reserved_content > 0 + + def can_receive(self): + return self.level + self.reserved_vacancy < self.capacity + + # def get_available_space(self): + # return self.capacity - self.level + + # def get_available_parts(self): + # return self.level + + def define_routing(self, upstream=[], downstream=[]): + self.upstream = upstream + self.downstream = downstream + + def get_candidate_givers(self, only_free=False, blocked=False): + if blocked: + # get only candidate givers that can give a part + return [obj for obj in self.get_candidate_givers() if obj.blocked] + else: + return [obj for obj in self.upstream if obj.can_give()] + + def get_candidate_receivers(self, only_free=False, starved=False): + if starved: + return [obj for obj in self.get_candidate_receivers() if obj.starved] + else: + # get only candidate receivers that can accept a part + return [obj for obj in self.downstream if obj.can_receive()] + \ No newline at end of file diff --git a/simantha/Machine.py b/simantha/Machine.py index 96e25db..e6d78ea 100644 --- a/simantha/Machine.py +++ b/simantha/Machine.py @@ -1,359 +1,420 @@ +import random import time +import warnings -import numpy as np -from scipy import stats -import simpy - - -class Machine: - degradation_matrix = np.array([[0.9, 0.1, 0., 0., 0., 0., 0., 0., 0., 0., 0.], - [0., 0.9, 0.1, 0., 0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0.9, 0.1, 0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0.9, 0.1, 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0.9, 0.1, 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0.9, 0.1, 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0.9, 0.1, 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0., 0.9, 0.1, 0., 0.], - [0., 0., 0., 0., 0., 0., 0., 0., 0.9, 0.1, 0.], - [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.9, 0.1], - [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]]) +from .Asset import Asset +from .simulation import * +class Machine(Asset): def __init__( - self, - env, - system, - index, - cycle_time, - in_buffer=None, - out_buffer=None, - degradation_matrix=degradation_matrix, - maintenance_threshold=999, - repairman=None, - initial_health=0, - initial_remaining_process=None, - allow_new_failures=True - ): + self, + name=None, + cycle_time=1, + selection_priority=1, - self.env = env - self.system = system + degradation_matrix=[[1,0],[0,1]], # By default, never degrade + cbm_threshold=None, + planned_failure=None, # Optional planned failure, in the form of (time, duration) + + pm_distribution=5, + cm_distribution=10, - self.index = index + # Initial machine state + initial_health=0, + initial_remaining_process=None + ): + # User-specified parameters + self.name = name - self.cycle_time = cycle_time - self.in_buffer = in_buffer - self.out_buffer = out_buffer - self.degradation_matrix = degradation_matrix - self.failed_state = degradation_matrix.shape[0] - 1 - self.allow_new_failures = allow_new_failures - - # maintenance parameters - self.maintenance_threshold = min( - [self.failed_state, maintenance_threshold] - ) + self.cycle_time = Distribution(cycle_time) - # production state - self.parts_made = 0 - if initial_remaining_process: - self.remaining_processing_time = initial_remaining_process - self.has_part = True + if initial_remaining_process is not None: + self.initial_remaining_process = initial_remaining_process else: - self.remaining_processing_time = cycle_time - self.has_part = False + self.initial_remaining_process = self.get_cycle_time() + self.remaining_process_time = self.initial_remaining_process + self.selection_priority = selection_priority - # maintenance state + # Initial machine state + self.has_finished_part = False + self.initial_health = initial_health self.health = initial_health - self.health_data = np.array([[0, 0]]) - self.maintenance_state = 'healthy' # healthy, unhealthy, failed - - if self.health >= self.maintenance_threshold: - self.in_queue = True - self.allow_new_failures = True + self.degradation_matrix = degradation_matrix + self.failed_health = len(degradation_matrix) - 1 + self.cbm_threshold = cbm_threshold or self.failed_health # if not specified, CM is used + if self.health == self.failed_health: + self.failed = True else: - self.in_queue = False - self.allow_new_failures = allow_new_failures + self.failed = False + self.assigned_maintenance = False + self.maintainer = None + self.has_reserved_part = False - self.failed = False - self.under_repair = False - self.repair_type = None - self.time_entered_queue = np.inf - self.repairman = repairman - self.repair_event = self.env.event() + self.pm_distribution = Distribution(pm_distribution) + self.cm_distribution = Distribution(cm_distribution) + + self.planned_failure = planned_failure + + # check if planned failures and degradation are specified (may cause errors) + if planned_failure is not None and degradation_matrix[0][0] != 1: + warnings.warn( + 'Specifying planned failures along with degradtion is untested and may cause errors.' + ) - # machine data - self.production_data = [] - self.maintenance_data = [] + # Routing + self.upstream = [] + self.downstream = [] - self.production_process = self.env.process(self.production()) - self.degradation_process = self.env.process(self.degradation()) - - - def production(self): - # if failed, should sleep until woken up - while True: - try: - if self.in_buffer and (not self.has_part): - yield self.in_buffer.get(1) - self.update_in_buffer_data(self.env.now) - self.has_part = True - - # TEMP DEBUG - #if (self.index == self.system.n - 1) and (not self.system.mcts_system): - # print(f'M{self.index} got part at t={self.env.now}') - - while self.remaining_processing_time: - yield self.env.timeout(1) - self.remaining_processing_time -= 1 - self.remaining_processing_time = self.cycle_time - - if self.out_buffer: - yield self.out_buffer.put(1) - self.update_out_buffer_data(self.env.now) - self.has_part = False - - # TEMP DEBUG - #if (self.index == self.system.n - 1) and (not self.system.mcts_system): - # print(f'M{self.index} put part at t={self.env.now}') - - if ( - (self.env.now >= self.system.warm_up_time) - or (self.system.mcts_system) - ): - self.parts_made += 1 - self.production_data.append( - [self.env.now-self.system.warm_up_time, self.parts_made] - ) - - # TEMP DEBUG - #if (self.index == self.system.n - 1) and (not self.system.mcts_system): - # if self.system.mcts_system: - # print('MCTS: ', end='') - # print(f'M{self.index} finished part at t={self.env.now}, total production: {self.parts_made}') - - except simpy.Interrupt: - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} stopped production at t={self.env.now}') - - yield self.repair_event # wait for repair to finish - - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} repaired at t={self.env.now}') - - self.repair_event = self.env.event() - self.degradation_process = self.env.process(self.degradation()) - - self.system.repairman.schedule_maintenance() - - - def degradation(self): - try: - if self.maintenance_threshold < self.failed_state: - # time until enter maintenance queue - ttq = self.get_time_to_queue() - yield self.env.timeout(ttq) - - self.repair_type = 'preventive' - self.time_entered_queue = self.env.now - self.in_queue = True - - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} entered queue at t={self.env.now} with health {self.get_health(self.env.now)}') - - self.system.repairman.schedule_maintenance() - - # time until complete failure - ttf = self.get_time_to_failure() - yield self.env.timeout(ttf) - self.repair_type = 'corrective' - self.failed = True - if not self.in_queue: - self.in_queue = True - self.time_entered_queue = self.env.now - self.system.repairman.schedule_maintenance() - - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} failed at t={self.env.now}') + # Machine status + self.has_part = False + self.under_repair = False + self.in_queue = False + self.remaining_ttr = None + + # Machine statistics + self.parts_made = 0 + self.downtime = 0 - self.production_process.interrupt() + # Simulation data + self.production_data = {'time': [0], 'production': [0]} + self.health_data = {'time': [0], 'health': [self.health]} + self.maintenance_data = {'time': [], 'event': []} + + self.env = None - except simpy.Interrupt: - # degradation interrupted if maintenance is scheduled before failure - yield self.repair_event + def initialize(self): + self.remaining_process_time = self.initial_remaining_process + self.health = self.initial_health + if self.health == self.failed_health: + self.failed = True + else: + self.failed = False + self.time_entered_queue = -1 - - def repair(self): + self.has_part = False + self.under_repair = False self.in_queue = False - self.under_repair = True + self.remaining_ttr = None - self.system.repairman.update_queue_data() + self.target_giver = None + self.target_receiver = None - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} starting repair at t={self.env.now}') + self.reserved_content = 0 + self.reserved_vacancy = 0 - #self.maintenance_request = self.system.repairman.request(priority=1) - if self.repair_type == 'preventive': - self.degradation_process.interrupt() - self.production_process.interrupt() + self.blocked = False + self.starved = True + + # Initialize statistics + self.parts_made = 0 + self.downtime = 0 - # clear health data beyond this point - self.health_data = self.health_data[self.health_data[:,0] < self.env.now] + # Schedule planned failures + if self.planned_failure is not None: + self.env.schedule_event( + self.planned_failure[0], self, self.maintain_planned_failure + ) - ttr = self.get_time_to_repair() - self.maintenance_data.append([self.env.now-self.system.warm_up_time, - self.repair_type, ttr]) - yield self.env.timeout(ttr) + # Initialize data + if self.env.collect_data: + self.production_data = {'time': [0], 'production': [0]} + self.health_data = {'time': [0], 'health': [self.health]} + self.maintenance_data = {'time': [], 'event': []} - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} repair completed at t={self.env.now}') + # Schedule initial events + time_to_degrade = self.get_time_to_degrade() + self.env.schedule_event( + time_to_degrade, self, self.degrade, f'{self.name}.initialize' + ) - self.health = 0 # machine completely restored - self.update_health_data(at=self.env.now) + self.initialize_addon_processes() - self.failed = False - self.under_repair = False - self.remaining_processing_time = self.cycle_time + def initialize_addon_processes(self): + pass + + def reserve_vacancy(self, quantity=1): + self.reserved_vacancy += 1 - if not self.system.allow_new_failures: - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} disallowing additional failures at t={self.env.now}') + def get_part(self): + # Choose a random upstream container from which to take a part. + assert self.target_giver is not None, f'No giver identified for {self.name}' + self.target_giver.get(1) - self.allow_new_failures = False - self.in_queue = False + self.has_part = True - self.repair_event.succeed() - self.system.repairman.utilization -= 1 # release repairman + self.env.schedule_event( + self.env.now+self.get_cycle_time(), + self, + self.request_space, + f'{self.name}.get_part at {self.env.now}' + ) - #self.system.repairman.update_queue_data() + # check if this event unblocked another machine + for asset in self.target_giver.upstream: + if asset.can_give() and self.target_giver.can_receive(): + source = f'{self.name}.get_part at {self.env.now}' + self.env.schedule_event( + self.env.now, asset, asset.request_space, source + ) - def get_time_to_repair(self): - if self.get_health(self.env.now) < self.failed_state: - self.repair_type = 'preventive' + self.target_giver = None + + def request_space(self): + #request_space_start = time.time() + self.has_finished_part = True + candidate_receivers = [obj for obj in self.downstream if obj.can_receive()] + if len(candidate_receivers) > 0: + self.target_receiver = random.choice(candidate_receivers) + self.target_receiver.reserve_vacancy(1) + source = f'{self.name}.request_space at {self.env.now}' + self.env.schedule_event(self.env.now, self, self.put_part, source) else: - self.repair_type = 'corrective' + self.blocked = True + + def put_part(self): + assert self.target_receiver is not None, f'No receiver identified for {self.name}' - if self.repair_type == 'preventive': - ttr = self.system.pm_distribution.rvs() - elif self.repair_type == 'corrective': - ttr = self.system.cm_distribution.rvs() + self.target_receiver.put(1) + + if self.env.now > self.env.warm_up_time: + self.parts_made += 1 + self.has_finished_part = False + self.has_part = False + + if self.env.now > self.env.warm_up_time and self.env.collect_data: + self.production_data['time'].append(self.env.now) + self.production_data['production'].append(self.parts_made) + + source = f'{self.name}.put_part at {self.env.now}' + self.env.schedule_event(self.env.now, self, self.request_part, source) + + # check if this event fed another machine + for asset in self.target_receiver.downstream: + if self.target_receiver.can_give() and asset.can_receive() and not asset.has_content_request(): + source = f'{self.name}.put_part at {self.env.now}' + self.env.schedule_event(self.env.now, asset, asset.request_part, source) + + self.target_receiver = None + + def request_part(self): + candidate_givers = [obj for obj in self.upstream if obj.can_give()] + if len(candidate_givers) > 0: + self.starved = False + self.target_giver = random.choice(candidate_givers) + self.target_giver.reserve_content(1) + source = f'{self.name}.request_part at {self.env.now}' + self.env.schedule_event(self.env.now, self, self.get_part, source) else: - raise ValueError('Invalid repair type') - - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} {self.repair_type} TTR:', ttr) - - return ttr - - - def get_time_to_queue(self): - # time to enter the maintenance queue - #current_health = self.health - if ( - (1 in np.diagonal(self.degradation_matrix)[:-1]) - or (not self.allow_new_failures) - ): - #if self.system.debug: - # print(f'M{self.index} TTQ is inf') - ttq = np.inf + self.starved = True + + def degrade(self): + source = f'{self.name}.degrade at {self.env.now}' + self.health += 1 + + if self.env.collect_data: + self.health_data['time'].append(self.env.now) + self.health_data['health'].append(self.health) + + time_to_degrade = self.get_time_to_degrade() + if self.health == self.failed_health: + self.env.schedule_event(self.env.now, self, self.fail, source) + elif self.health == self.cbm_threshold: + self.env.schedule_event(self.env.now, self, self.enter_queue, source) + self.env.schedule_event( + self.env.now+time_to_degrade, self, self.degrade, source + ) else: - ttq = 0 - while self.health < self.maintenance_threshold: # TODO: verify '<' here - previous_health = self.health - all_states = np.arange(self.failed_state + 1) - self.health = np.random.choice( - all_states, p=self.degradation_matrix[self.health] - ) - if self.health != previous_health: - self.update_health_data(at=self.env.now+ttq) - ttq += 1 + self.env.schedule_event( + self.env.now+time_to_degrade, self, self.degrade, source + ) - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} TTQ: {ttq} at t={self.env.now}') + def enter_queue(self): + if not self.in_queue: + if self.env.collect_data: + self.maintenance_data['time'].append(self.env.now) + self.maintenance_data['event'].append('enter queue') - return ttq - - - def get_time_to_failure(self): - if ( - (1 in np.diagonal(self.degradation_matrix)[:-1]) - or ((not self.allow_new_failures) and (not self.in_queue)) - ): - #if self.system.debug: - # print(f'M{self.index} TTF is inf') - ttf = np.inf - else: - ttf = 0 - while self.health != self.failed_state: - previous_health = self.health - all_states = np.arange(self.failed_state + 1) - self.health = np.random.choice( - all_states, p=self.degradation_matrix[self.health] - ) - if self.health != previous_health: - self.update_health_data(at=self.env.now+ttf) - ttf +=1 + self.time_entered_queue = self.env.now + self.in_queue = True + + if not self.failed and self.maintainer.is_available(): + source = f'{self.name}.enter_queue at {self.env.now}' + self.env.schedule_event( + self.env.now, self.maintainer, self.maintainer.inspect, source + ) + + def fail(self): + self.failed = True + self.downtime_start = self.env.now + + if not self.in_queue: + self.enter_queue() - if self.system.debug: - if self.system.mcts_system: - print('MCTS: ', end='') - print(f'M{self.index} TTF: {ttf} at t={self.env.now}') + if self.env.collect_data: + self.maintenance_data['time'].append(self.env.now) + self.maintenance_data['event'].append('failure') - return ttf + self.cancel_all_events() + if self.maintainer.is_available(): + source = f'{self.name}.fail at {self.env.now}' + self.env.schedule_event( + self.env.now, self.maintainer, self.maintainer.inspect, source + ) + + def get_cycle_time(self): + return self.cycle_time.sample() + + def get_time_to_degrade(self): + if 1 in self.degradation_matrix[self.health]: + return float('inf') + + ttd = 0 + next_health = self.health + while next_health == self.health: + ttd += 1 + next_health = random.choices( + population=range(self.failed_health+1), + weights=self.degradation_matrix[self.health], + k=1 + )[0] + return ttd - def update_health_data(self, at=None): - self.health_data = np.append( - self.health_data, [[at, self.health]], axis=0 + def maintain(self): + if not self.failed: + self.downtime_start = self.env.now + self.has_part = False + self.has_finished_part = False + self.under_repair = True + + if self.env.collect_data: + self.maintenance_data['time'].append(self.env.now) + self.maintenance_data['event'].append('begin maintenance') + + self.in_queue = False + time_to_repair = self.get_time_to_repair() + + self.cancel_all_events() + + source = f'{self.name}.maintain at {self.env.now}' + self.env.schedule_event(self.env.now+time_to_repair, self, self.restore, source) + + def maintain_planned_failure(self): + self.failed = True + self.downtime_start = self.env.now + self.under_repair = True + + if self.env.collect_data: + self.maintenance_data['time'].append(self.env.now) + self.maintenance_data['event'].append('planned failure') + + self.cancel_all_events() + + time_to_repair = self.planned_failure[1] + source = f'{self.name}.maintain_planned_failure at {self.env.now}' + self.env.schedule_event( + self.env.now+time_to_repair, self, self.restore, source ) + def restore(self): + self.health = 0 + self.under_repair = False + self.failed = False + + self.maintainer.utilization -= 1 - def update_in_buffer_data(self, at=None): - if at in self.in_buffer.buffer_data[:,0]: - i = np.where(self.in_buffer.buffer_data[:,0] == at) - self.in_buffer.buffer_data[i] = [at, self.in_buffer.level] - else: - self.in_buffer.buffer_data = np.append( - self.in_buffer.buffer_data, - [[at, self.in_buffer.level]], - axis=0 - ) + self.downtime += (self.env.now - self.downtime_start) + if self.env.collect_data: + self.maintenance_data['time'].append(self.env.now) + self.maintenance_data['event'].append('repaired') - def update_out_buffer_data(self, at=None): - if at in self.out_buffer.buffer_data[:,0]: - i = np.where(self.out_buffer.buffer_data[:,0] == at) - self.out_buffer.buffer_data[i] = [at, self.out_buffer.level] - else: - self.out_buffer.buffer_data = np.append( - self.out_buffer.buffer_data, - [[at, self.out_buffer.level]], - axis=0 - ) + self.health_data['time'].append(self.env.now) + self.health_data['health'].append(self.health) + + source = f'{self.name}.restore at {self.env.now}' + self.env.schedule_event(self.env.now, self, self.request_part, source) + time_to_degrade = self.get_time_to_degrade() + self.env.schedule_event( + self.env.now+time_to_degrade, self, self.degrade, source + ) + + # Repairman to scan queue once released + self.env.schedule_event( + self.env.now, self.maintainer, self.maintainer.inspect, source + ) + self.repair_addon_processes() - def get_health(self, time): - # returns machine health at a specific time - time_steps = self.health_data[:,0] + def repair_addon_processes(self): + pass + + def requesting_maintenance(self): + return ( + (not self.under_repair) + and ((self.failed) or (self.health >= self.cbm_threshold)) + and (not self.assigned_maintenance) + ) - if time in time_steps: - index = np.where(self.health_data[:,0] == time)[0][0] + def get_time_to_repair(self): + if self.failed: + return self.cm_distribution.sample() else: - index = np.searchsorted(time_steps, time) - 1 + return self.pm_distribution.sample() - return self.health_data[index, 1] + def define_routing(self, upstream=[], downstream=[]): + self.upstream = upstream + self.downstream = downstream + + def can_receive(self): + return ( + (not self.under_repair) + and (not self.failed) + and (not self.has_part) + ) + + def can_give(self): + return ( + (self.has_finished_part) + and (not self.under_repair) + and (not self.failed) + ) or ( + (self.has_finished_part) + and (self.downtime_start == self.env.now) + ) + + def has_content_request(self): + # check if a machine has an existing request for a part + for event in self.env.events: + if ( + ((event.location is self) and (event.action.__name__ == 'request_part')) + or ((event.location is self) and (event.action.__name__ == 'get_part')) + ): + return True + return False + + def has_vacancy_request(self): + for event in self.env.events: + if (event.location is self) and (event.action.__name__ == 'request_space'): + return True + return False + + def cancel_all_events(self): + # cancel all events scheduled on this machine + for event in self.env.events: + if event.location == self: + event.canceled = True + + def get_candidate_givers(self, only_free=False, blocked=False): + if blocked: + # get only candidate givers that can give a part + return [obj for obj in self.get_candidate_givers() if obj.blocked] + else: + return [obj for obj in self.upstream if obj.can_give()] + + def get_candidate_receivers(self, only_free=False, starved=False): + if starved: + return [obj for obj in self.get_candidate_receivers() if obj.starved] + else: + # get only candidate receivers that can accept a part + return [obj for obj in self.downstream if obj.can_receive()] diff --git a/simantha/Maintainer.py b/simantha/Maintainer.py new file mode 100644 index 0000000..4b7aa38 --- /dev/null +++ b/simantha/Maintainer.py @@ -0,0 +1,43 @@ +import random + +from .simulation import * + +class Maintainer: + # Basic maintainer, follows FIFO by default + def __init__(self, name='repairman', capacity=float('inf')): + self.name = name + self.capacity = capacity + + self.utilization = 0 + + self.env = None + self.system = None + + def initialize(self): + self.utilization = 0 + + def is_available(self): + return self.utilization < self.capacity + + def inspect(self): + # if available, check for machines requesting repair + current_queue = self.get_queue() + if (not self.is_available()) or (len(current_queue) == 0): + # No available capacity and/or empty queue + return + else: + machine = self.choose_maintenance_action(current_queue) + self.utilization += 1 + machine.in_queue = False + machine.under_repair = True + source = f'{self.name}.inspect at {self.env.now}' + self.env.schedule_event(self.env.now, machine, machine.maintain, source) + + def choose_maintenance_action(self, queue): + # default fifo policy, break ties randomly + earliest_request = min(m.time_entered_queue for m in queue) + candidates = [m for m in queue if m.time_entered_queue == earliest_request] + return random.choice(candidates) + + def get_queue(self): + return [machine for machine in self.system.machines if machine.in_queue] diff --git a/simantha/Sink.py b/simantha/Sink.py new file mode 100644 index 0000000..1125645 --- /dev/null +++ b/simantha/Sink.py @@ -0,0 +1,35 @@ +class Sink: + def __init__(self, name='Sink', initial_level=0): + self.name = name + self.capacity = float('inf') + self.initial_level = initial_level + self.level = initial_level + + self.selection_priority = 1 + + self.env = None + + # self.level_data = {'time': [0], 'level': [initial_level]} + + def initialize(self): + self.level = self.initial_level + + def reserve_vacancy(self, quantity=1): + return + + def put(self, quantity=1): + if self.env.now > self.env.warm_up_time: + self.level += quantity + + # self.level_data['time'].append(self.env.now) + # self.level_data['level'].append(self.level) + + def define_routing(self, upstream=[], downstream=[]): + self.upstream = upstream + self.downstream = downstream + + def can_give(self): + return False + + def can_receive(self): + return True diff --git a/simantha/Source.py b/simantha/Source.py new file mode 100644 index 0000000..6d64759 --- /dev/null +++ b/simantha/Source.py @@ -0,0 +1,78 @@ +import random + +class Source: + def __init__( + self, + name='Source', + interarrival_time=None + ): + self.name = name + self.interarrival_time = interarrival_time + self.last_arrival = 0 + + if self.interarrival_time is None: + self.level = float('inf') + else: + self.level = 0 + + self.define_routing() + + self.env = None + + def initialize(self): + # schedule inital get events for each downstream machine + # currently it's assumed that no machine pulling from a source is starved + self.reserved_content = 0 + + for receiver in self.downstream: + if receiver.can_receive(): + receiver.starved = False + self.env.schedule_event( + self.env.now, + receiver, + receiver.request_part, + f'{self.name}.arrival at {self.env.now}' + ) + + def generate_arrival(self): + if self.interarrival_time is None: + return + + self.last_arrival += 1 + if (self.last_arrival >= self.interarrival_time) and self.is_empty(): + self.level += 1 + self.last_arrival = 0 + + def get(self, quantity=1): + self.reserved_content -= quantity + if self.interarrival_time is None: + return + + if not self.is_empty(): + self.level -= quantity + else: + raise RuntimeError('Attempting to take part from source before arrival.') + + def reserve_content(self, quantity=1): + self.reserved_content += quantity + + def is_empty(self): + return self.level == 0 + + def define_routing(self, upstream=[], downstream=[]): + self.upstream = upstream + self.downstream = downstream + + def can_give(self): + # TODO: this assumes receivers are never starved + return True + + def get_candidate_givers(self, only_free=False, blocked=False): + return self.upstream + + def get_candidate_receivers(self, only_free=False, starved=False): + if only_free: + # get only candidate receivers that can accept a part + return [obj for obj in self.get_candidate_receivers() if obj.can_receive()] + else: + return [obj for obj in self.downstream if obj.can_receive()] \ No newline at end of file diff --git a/simantha/System.py b/simantha/System.py index d824c6f..5f4d038 100644 --- a/simantha/System.py +++ b/simantha/System.py @@ -1,4 +1,5 @@ import copy +<<<<<<< HEAD import time from joblib import Parallel, delayed @@ -248,11 +249,110 @@ def continue_simulation(self, simulation_time, debug=False): self.env.run(until=self.env.now+simulation_time) +======= +import multiprocessing +import random +import time +import warnings + +from .simulation import Environment +from .Source import Source +from .Sink import Sink +from .Machine import Machine +from .Buffer import Buffer +from .Maintainer import Maintainer + +class System: + def __init__( + self, + objects=[], + maintainer=Maintainer() + ): + self.objects = objects + self.sources = [] + self.machines = [] + self.buffers = [] + self.sinks = [] + for obj in objects: + if type(obj) == Source: + self.sources.append(obj) + #elif type(obj) == Machine: + elif isinstance(obj, Machine): + self.machines.append(obj) + elif type(obj) == Buffer: + self.buffers.append(obj) + elif type(obj) == Sink: + self.sinks.append(obj) + self.maintainer = maintainer + + # put machines at the front as they should be initialized first + self.objects.sort(key=lambda obj: not isinstance(obj, Machine)) + + def initialize(self): + for machine in self.machines: + machine.remaining_process_time = machine.initial_remaining_process + machine.parts_made = 0 + machine.health = machine.initial_health + + if len(machine.upstream) > 1 or len(machine.downstream) > 1: + warnings.warn( + 'System configuration includes machines with more than one asset ' + + 'upstream or downstream, which may result in unexepcted behavior. ' + + 'It is recommended that the system be rearranged such that each ' + + 'machine gives and takes from only one buffer.' + ) + + for buffer in self.buffers: + buffer.level = buffer.initial_level + + for sink in self.sinks: + sink.level = sink.initial_level + + self.maintainer.utilization = 0 + + def simulate( + self, + warm_up_time=0, + simulation_time=0, + verbose=True, + trace=False, + collect_data=True + ): + start = time.time() + for machine in self.machines: + machine.maintainer = self.maintainer + + self.env = Environment(trace=trace, collect_data=collect_data) + for obj in self.objects: + # should initialize machines first + obj.env = self.env + obj.initialize() + + self.maintainer.env = self.env + self.maintainer.system = self + self.maintainer.utilization = 0 + + self.warm_up_time = warm_up_time + self.simulation_time = simulation_time + + self.env.run(warm_up_time, simulation_time) + + # clean up data here + for machine in self.machines: + if machine.under_repair or machine.failed: + machine.downtime += (self.env.now - machine.downtime_start) + + stop = time.time() + if verbose: + print(f'Simulation finished in {stop-start:.2f}s') + print(f'Parts produced: {sum([sink.level for sink in self.sinks])}') +>>>>>>> dev def iterate_simulation( self, replications, warm_up_time=0, +<<<<<<< HEAD simulation_time=100, objective='production', verbose=True @@ -323,3 +423,72 @@ def main(): if __name__ == '__main__': main() +======= + simulation_time=0, + store_system_state=False, + verbose=True, + jobs=1, + seedseed=0 + ): + """Replicate multiple simulation runs for a specified system. Statistics for + each run will gathered after the "warm_up_time" has elapsed. Currently the + following statistics are gathered: + - Machine + - Production (units) + - Availability (proportion of time not failed or under maintenance) + - Sink + - Level (units): completed parts that have exited the system + + A nested dictionary is returned with "replications" samples of each statistic. + """ + start = time.time() + with multiprocessing.Pool(jobs) as p: + args = [ + (seed, warm_up_time, simulation_time, store_system_state) + for seed in range(seedseed, seedseed+replications) + ] + samples = p.starmap(self.simulate_in_parallel, args) + stop = time.time() + + if verbose: + print(f'Finished {replications} replications in {stop-start:.2f}s') + + return samples + + def simulate_in_parallel( + self, + seed, + warm_up_time, + simulation_time, + store_system_state=False + ): + random.seed() + + self.simulate( + warm_up_time, + simulation_time, + verbose=False, + collect_data=store_system_state + ) + + availability = [ + (1 - machine.downtime/(warm_up_time+simulation_time)) + for machine in self.machines + ] + + machine_production = [machine.parts_made for machine in self.machines] + + system_production = sum([sink.level for sink in self.sinks]) + + if store_system_state: + system_state = copy.deepcopy(self) + else: + system_state = None + + return ( + system_production, + machine_production, + availability, + system_state + ) +>>>>>>> dev diff --git a/simantha/__init__.py b/simantha/__init__.py index f722b4b..6a1b4a8 100644 --- a/simantha/__init__.py +++ b/simantha/__init__.py @@ -1,4 +1,22 @@ +<<<<<<< HEAD from .Evaluator import * from .Machine import * from .Repairman import * -from .System import * \ No newline at end of file +from .System import * +======= +from .System import System +from .Source import * +from .Machine import * +from .Buffer import * +from .Sink import * +from .Maintainer import * +from .simulation import * +from .utils import * + +#__name__ = 'simantha' + +#__all__ = ['System', 'Maintainer', 'Machine'] + +#__title__ = 'simantha' +#__author__ = 'Michael Hoffman ' +>>>>>>> dev diff --git a/simantha/examples/ConditionBasedMaintenance.py b/simantha/examples/ConditionBasedMaintenance.py new file mode 100644 index 0000000..23e0c55 --- /dev/null +++ b/simantha/examples/ConditionBasedMaintenance.py @@ -0,0 +1,58 @@ +""" +Example of a condition-based maintenance policy. The CBM threshold determines the health +index level at which a machine requests maintenance. +""" + +import random + +from simantha import Source, Machine, Buffer, Sink, Maintainer, System, utils + +def main(): + degradation_matrix = [ + [0.9, 0.1, 0., 0., 0. ], + [0., 0.9, 0.1, 0., 0. ], + [0., 0., 0.9, 0.1, 0. ], + [0., 0., 0., 0.9, 0.1], + [0., 0., 0., 0., 1. ] + ] + cbm_threshold = 3 + pm_distribution = {'geometric': 0.25} + cm_distribution = {'geometric': 0.10} + + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cbm_threshold=cbm_threshold, + pm_distribution=pm_distribution, + cm_distribution=cm_distribution + ) + B1 = Buffer(capacity=10) + M2 = Machine( + name='M2', + cycle_time=2, + degradation_matrix=degradation_matrix, + cbm_threshold=cbm_threshold, + pm_distribution=pm_distribution, + cm_distribution=cm_distribution + ) + sink = Sink() + + objects = [source, M1, B1, M2, sink] + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + maintainer = Maintainer(capacity=1) + + system = System(objects=objects, maintainer=maintainer) + + random.seed(1) + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/MachineExtensibility.py b/simantha/examples/MachineExtensibility.py new file mode 100644 index 0000000..396bf17 --- /dev/null +++ b/simantha/examples/MachineExtensibility.py @@ -0,0 +1,95 @@ +""" +An example of extensible machine behavior. Random noise is added to machine health index +readings and monitored periodically instead of continuously. +""" + +import random + +import numpy as np + +from simantha import Source, Machine, Sink, Maintainer, System, simulation, utils + +class SensingEvent(simulation.Event): + # A new event type is used to assign the correct priority to the sensing event. In + # this case, sensing events are scheduled just after machine degradation events. + def get_action_priority(self): + return 5.5 + +class ConditionMonitoredMachine(Machine): + def __init__(self, sensing_interval=1, sensor_noise=0, **kwargs): + self.sensing_interval = sensing_interval + self.sensor_noise = sensor_noise + self.sensor_data = {'time': [], 'reading': []} + + super().__init__(**kwargs) + + def initialize_addon_processes(self): + self.env.schedule_event( + time=self.env.now, + location=self, + action=self.sense, + source=f'{self.name} initial addon process', + event_type=SensingEvent + ) + + def repair_addon_processes(self): + self.env.schedule_event( + time=self.env.now, + location=self, + action=self.sense, + source=f'{self.name} repair addon process at {self.env.now}', + event_type=SensingEvent + ) + + def sense(self): + self.sensor_reading = self.health + np.random.normal(0, self.sensor_noise) + + self.sensor_data['time'].append(self.env.now) + self.sensor_data['reading'].append(self.sensor_reading) + + self.env.schedule_event( + time=self.env.now+self.sensing_interval, + location=self, + action=self.sense, + source=f'{self.name} sensing at {self.env.now}', + event_type=SensingEvent + ) + +def main(): + degradation_matrix = utils.generate_degradation_matrix(h_max=10, p=0.1) + cm_distribution = {'geometric': 0.1} + + source = Source() + M1 = ConditionMonitoredMachine( + name='M1', + cycle_time=2, + degradation_matrix=degradation_matrix, + cm_distribution=cm_distribution, + sensing_interval=2, + sensor_noise=1 + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1]) + + system = System(objects=[source, M1, sink]) + + random.seed(1) + system.simulate(simulation_time=6*60) + + # Print true health and corresponding sensor reading + rows = 12 + print('\ntime health sensor reading') + for time, reading in zip( + M1.sensor_data['time'][:rows], M1.sensor_data['reading'][:rows] + ): + timestamp = max([t for t in M1.health_data['time'] if t <= time]) + idx = M1.health_data['time'].index(timestamp) + health = M1.health_data['health'][idx] + print(f'{time:<4} {health:<3} {reading:>8.4f}') + + +if __name__ == '__main__': + main() diff --git a/simantha/examples/MaintainerExtensibility.py b/simantha/examples/MaintainerExtensibility.py new file mode 100644 index 0000000..c5e2022 --- /dev/null +++ b/simantha/examples/MaintainerExtensibility.py @@ -0,0 +1,70 @@ +""" +An example of the extensibility of the Maintainer class. A "least processing time first" +(LPT) policy is implemented. +""" + +import random + +from simantha import Source, Machine, Buffer, Sink, Maintainer, System, utils + +class LptMaintainer(Maintainer): + """Chooses the maintenance action with the longest expected duration first.""" + def choose_maintenance_action(self, queue): + def expected_repair_time(machine): + if machine.failed: # Machine requires corrective maintenance + return machine.cm_distribution.mean + else: # Machine requires preventive maintenance + return machine.pm_distribution.mean + + return max(queue, key=expected_repair_time) + +def main(): + degradation_matrix = [ + [0.9, 0.1, 0., 0., 0. ], + [0., 0.9, 0.1, 0., 0. ], + [0., 0., 0.9, 0.1, 0. ], + [0., 0., 0., 0.9, 0.1], + [0., 0., 0., 0., 1. ] + ] + + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.1} + ) + B1 = Buffer(capacity=10) + M2 = Machine( + name='M2', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.075} + ) + B2 = Buffer(capacity=10) + M3 = Machine( + name='M3', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.05} + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[B2]) + B2.define_routing(upstream=[M2], downstream=[M3]) + M3.define_routing(upstream=[B2], downstream=[sink]) + sink.define_routing(upstream=[M3]) + + objects = [source, M1, B1, M2, B2, M3, sink] + maintainer = LptMaintainer() + + system = System(objects, maintainer) + + random.seed(1) + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/ParallelStations.py b/simantha/examples/ParallelStations.py new file mode 100644 index 0000000..46db343 --- /dev/null +++ b/simantha/examples/ParallelStations.py @@ -0,0 +1,42 @@ +from simantha import Source, Machine, Buffer, Sink, System, utils + +def main(): + source = Source() + + M1 = Machine('M1', cycle_time=2) + M2 = Machine('M2', cycle_time=2) + station1 = [M1, M2] + + B1 = Buffer(capacity=5) + + M3 = Machine('M3', cycle_time=3) + M4 = Machine('M4', cycle_time=3) + M5 = Machine('M5', cycle_time=3) + station2 = [M3, M4, M5] + + B2 = Buffer(capacity=5) + + M6 = Machine('M6', cycle_time=2) + M7 = Machine('M7', cycle_time=2) + station3 = [M6, M7] + + sink = Sink() + + source.define_routing(downstream=station1) + for machine in station1: + machine.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=station1, downstream=station2) + for machine in station2: + machine.define_routing(upstream=[B1], downstream=[B2]) + B2.define_routing(upstream=station2, downstream=station3) + for machine in station3: + machine.define_routing(upstream=[B2], downstream=[sink]) + sink.define_routing(upstream=station3) + + objects = [source] + station1 + [B1] + station2 + [B2] + station3 + [sink] + system = System(objects) + + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/SimulationReplication.py b/simantha/examples/SimulationReplication.py new file mode 100644 index 0000000..66baa17 --- /dev/null +++ b/simantha/examples/SimulationReplication.py @@ -0,0 +1,114 @@ +""" +This example compares the production of two systems, one under a corrective maintenance +policy and the other under a condition-based maintenance policy. The simulation +replication functionality is used to estimate the average production of both systems. +""" + +import random + +from simantha import Source, Machine, Buffer, Sink, Maintainer, System, utils + +def main(): + # Parameters for both systems + degradation_matrix = [ + [0.9, 0.1, 0., 0., 0. ], + [0., 0.9, 0.1, 0., 0. ], + [0., 0., 0.9, 0.1, 0. ], + [0., 0., 0., 0.9, 0.1], + [0., 0., 0., 0., 1. ] + ] + pm_distribution = {'geometric': 0.25} # Average PM duration: 4 minutes + cm_distribution = {'geometric': 0.10} # Average CM duration: 10 minutes + + # Define the corrective maintenance system + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution=cm_distribution + ) + B1 = Buffer(capacity=10) + M2 = Machine( + name='M2', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution=cm_distribution + ) + sink = Sink() + + objects = [source, M1, B1, M2, sink] + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + maintainer = Maintainer(capacity=1) + + corrective_maintenance_system = System(objects=objects, maintainer=maintainer) + + # Define the condition-based maintenance system + cbm_threshold = 3 + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cbm_threshold=cbm_threshold, + pm_distribution=pm_distribution, + cm_distribution=cm_distribution + ) + B1 = Buffer(capacity=10) + M2 = Machine( + name='M2', + cycle_time=1, + degradation_matrix=degradation_matrix, + cbm_threshold=cbm_threshold, + pm_distribution=pm_distribution, + cm_distribution=cm_distribution + ) + sink = Sink() + + objects = [source, M1, B1, M2, sink] + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + maintainer = Maintainer(capacity=1) + + condition_based_maintenance_system = System(objects=objects, maintainer=maintainer) + + # Simulate both systems + replications = 50 + random.seed(1) + + cm_results = corrective_maintenance_system.iterate_simulation( + replications=replications, + warm_up_time=utils.DAY, + simulation_time=utils.WEEK, + jobs=10, # Simulate in parallel using 10 cores + verbose=False + ) + cm_production = [r[0] for r in cm_results] + cm_average = sum(cm_production) / replications + + cbm_results = condition_based_maintenance_system.iterate_simulation( + replications=replications, + warm_up_time=utils.DAY, + simulation_time=utils.WEEK, + jobs=10, + verbose=False + ) + cbm_production = [r[0] for r in cbm_results] + cbm_average = sum(cbm_production) / replications + + print(f'Average corrective maintenance production: {cm_average:.2f} units') + print(f'Average condition-based maintenance production: {cbm_average:.2f} units') + +if __name__ == '__main__': + main() diff --git a/simantha/examples/SingleMachine.py b/simantha/examples/SingleMachine.py new file mode 100644 index 0000000..bde7512 --- /dev/null +++ b/simantha/examples/SingleMachine.py @@ -0,0 +1,18 @@ +from simantha import Source, Machine, Sink, System + +def main(): + source = Source() + M1 = Machine(name='M1', cycle_time=1) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1]) + + system = System(objects=[source, M1, sink]) + + system.simulate(simulation_time=100) + +if __name__ == '__main__': + main() + \ No newline at end of file diff --git a/simantha/examples/SingleMachineDegradation.py b/simantha/examples/SingleMachineDegradation.py new file mode 100644 index 0000000..77dd59a --- /dev/null +++ b/simantha/examples/SingleMachineDegradation.py @@ -0,0 +1,39 @@ +""" +An example of Markovian degradation of a single machine. The degradation matrix +specifies the health index transition probabilities. Once the machine reaches the +maximum health index it requires a corrective maintenance action to be restored. +""" + +import random + +from simantha import Source, Machine, Sink, System, utils + +def main(): + degradation_matrix = [ + [0.9, 0.1, 0., 0., 0. ], + [0., 0.9, 0.1, 0., 0. ], + [0., 0., 0.9, 0.1, 0. ], + [0., 0., 0., 0.9, 0.1], + [0., 0., 0., 0., 1. ] + ] + + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.1} + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1]) + + system = System([source, M1, sink]) + + random.seed(1) + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/TwoMachineOneBuffer.py b/simantha/examples/TwoMachineOneBuffer.py new file mode 100644 index 0000000..16713f5 --- /dev/null +++ b/simantha/examples/TwoMachineOneBuffer.py @@ -0,0 +1,21 @@ +from simantha import Source, Machine, Buffer, Sink, System, utils + +def main(): + source = Source() + M1 = Machine(name='M1', cycle_time=1) + B1 = Buffer(capacity=5) + M2 = Machine(name='M2', cycle_time=1) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + system = System([source, M1, B1, M2, sink]) + + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/TwoMachineOneBufferDegradation.py b/simantha/examples/TwoMachineOneBufferDegradation.py new file mode 100644 index 0000000..aee011a --- /dev/null +++ b/simantha/examples/TwoMachineOneBufferDegradation.py @@ -0,0 +1,43 @@ +import random + +from simantha import Source, Machine, Buffer, Sink, Maintainer, System, utils + +def main(): + degradation_matrix = [ + [0.9, 0.1, 0., 0., 0. ], + [0., 0.9, 0.1, 0., 0. ], + [0., 0., 0.9, 0.1, 0. ], + [0., 0., 0., 0.9, 0.1], + [0., 0., 0., 0., 1. ] + ] + + source = Source() + M1 = Machine( + name='M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.1} + ) + B1 = Buffer(capacity=5) + M2 = Machine( + name='M2', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'geometric': 0.1} + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + maintainer = Maintainer(capacity=1) + system = System([source, M1, B1, M2, sink], maintainer=maintainer) + + random.seed(1) + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/examples/TwoMachineParallel.py b/simantha/examples/TwoMachineParallel.py new file mode 100644 index 0000000..55c7f5d --- /dev/null +++ b/simantha/examples/TwoMachineParallel.py @@ -0,0 +1,19 @@ +from simantha import Source, Machine, Sink, System, utils + +def main(): + source = Source() + M1 = Machine(name='M1', cycle_time=1) + M2 = Machine(name='M2', cycle_time=1) + sink = Sink() + + source.define_routing(downstream=[M1, M2]) + M1.define_routing(upstream=[source], downstream=[sink]) + M2.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1, M2]) + + system = System([source, M1, M2, sink]) + + system.simulate(simulation_time=utils.WEEK) + +if __name__ == '__main__': + main() diff --git a/simantha/simulation.py b/simantha/simulation.py new file mode 100644 index 0000000..fae6a69 --- /dev/null +++ b/simantha/simulation.py @@ -0,0 +1,243 @@ +import bisect +import pickle +import random +import sys +import time +import warnings + +class Event: + action_priority = [ + # Events at the end of the last time step + 'generate_arrival', # Priority: 0 (highest priority) + 'request_space', # 1 + 'put_part', # 2 + 'restore', # 3 + + # Events at the start of the current time step + 'maintain_planned_failure', # 4 + 'degrade', # 5 + 'enter_queue', # 6 + 'fail', # 7 + 'inspect', # 8 + 'maintain', # 9 + 'request_part', # 10 + 'get_part', # 11 + + # Simulation runtime events + 'terminate' # 12 (lowest priority) + ] + + action_priority = { + action: priority for priority, action in enumerate(action_priority) + } + + def __init__(self, time, location, action, source='', priority=0, status=''): + self.time = time + self.location = location + self.action = action + self.source = source + self.priority = priority + self.status = status + + self.tiebreak = random.random() + + self.canceled = False + self.executed = False + + def get_action_priority(self): + if self.action.__name__ in self.action_priority.keys(): + return self.action_priority[self.action.__name__] + else: + return float('inf') + + def execute(self): + if not self.canceled: + self.action() + else: + self.status = 'canceled' + self.executed = True + + def __lt__(self, other): + return ( + self.time, + #self.action_priority[self.action.__name__], + self.get_action_priority(), + self.priority, + self.tiebreak + ) < ( + other.time, + #other.event_priority[other.action.__name__], + other.get_action_priority(), + other.priority, + other.tiebreak + ) + +class Environment: + """The main simulation environment for Simantha. This is designed to be an + enviroment specifically for use with Simantha objects and is not intended to be a + general simulation engine. In general, users of Simantha should not need to + instantiate an Environment object. + """ + def __init__(self, name='environment', trace=False, collect_data=True): + self.events = [] + self.name = name + self.now = 0 + + self.terminated = False + + self.trace = trace + if self.trace: + self.event_trace = { + 'time': [], + 'location': [], + 'action': [], + 'source': [], + 'priority': [], + 'status': [], + 'index': [] + } + + self.collect_data = collect_data + + def run(self, warm_up_time=0, simulation_time=0): + """Simulate the system for the specified run time or until no simulation events + remain. + """ + self.now = 0 + self.warm_up_time = warm_up_time + self.simulation_time = simulation_time + self.terminated = False + self.events.append(Event(warm_up_time+simulation_time, self, self.terminate)) + self.event_index = 0 + + self.events.sort() + + while self.events and not self.terminated: + self.step() + self.event_index += 1 + + if self.trace: + self.export_trace() + + def step(self): + """Find and execute the next earliest simulation event. Simultaneous events are + executed in order according to their event type priority, then their + user-assigned priority. If these values are equal then ties are broken randomly. + """ + next_event = self.events.pop(0) + + self.now = next_event.time + + try: + if self.trace: + self.trace_event(next_event) + next_event.execute() + except: + self.export_trace() + print('Failed event:') + print(f' time: {next_event.time}') + print(f' location: {next_event.location.name}') + print(f' action: {next_event.action.__name__}') + print(f' priority: {next_event.priority}') + sys.exit() + + def schedule_event( + self, time, location, action, source='', priority=0, event_type=Event + ): + new_event = Event(time, location, action, source, priority) + bisect.insort(self.events, new_event) + + def terminate(self): + self.terminated = True + + def trace_event(self, event): + if self.trace: + self.event_trace['time'].append(self.now) + self.event_trace['location'].append(event.location.name) + self.event_trace['action'].append(event.action.__name__) + self.event_trace['source'].append(event.source) + self.event_trace['priority'].append(event.priority) + if event.canceled: + self.event_trace['status'].append('canceled') + else: + self.event_trace['status'].append(event.status) + self.event_trace['index'].append(self.event_index) + + def export_trace(self): + if self.trace: + trace_file = open(f'{self.name}_trace.pkl', 'wb') + pickle.dump(self.event_trace, trace_file) + trace_file.close() + + +def rng(dist): + if 'constant' in dist.keys(): + # return deterministic value + return dist['constant'] + elif 'uniform' in dist.keys(): + # uniform distribution, specifed as + # {'unifom':[a, b]} + # returns a number between a and b, inclusive + a, b = dist['uniform'] + return random.randrange(a, b+1) + else: + raise NotImplementedError(f'Invalid distribution specified: {dist}') + +class Distribution: + """ + A class for random number generation in Simantha. Several built-in distributions are + available, but any class that returns a single integer value via a `sample` method + can be used. The built-in distributions are discrete uniform, specified by passing + {'uniform': [a, b]} to the distribution object, and geometric, specified via + {'geometric': p}. Constant integer values are also permitted. + """ + def __init__(self, distribution): + if type(distribution) == int: + self.distribution_type = 'constant' + self.distribution_parameters = distribution + elif type(distribution) != dict: + raise ValueError( + f'Invalid distribution {distribution}. Distribution should be a dictionary' + ) + elif len(distribution) > 1: + raise ValueError( + f'Invalid distribution {distribution}. Too many dictionary members' + ) + else: + for distribution_type, distribution_parameters in distribution.items(): + self.distribution_type = distribution_type + self.distribution_parameters = distribution_parameters + + if self.distribution_type == 'constant': + self.mean = self.distribution_parameters + elif self.distribution_type == 'uniform': + self.mean = sum(self.distribution_parameters) / 2 + elif self.distribution_type == 'geometric': + self.mean = 1 / self.distribution_parameters + else: + self.mean = None + + def sample(self): + """Returns a single sample from the specified distribution.""" + if self.distribution_type == 'constant': + return self.distribution_parameters + + elif self.distribution_type == 'uniform': + a, b = self.distribution_parameters + return random.randint(a, b) + + elif self.distribution_type == 'geometric': + # Returns the number of trials needed to achieve a single success, where the + # probability of success for each trial is p. + p = self.distribution_parameters + s = 1 + while random.random() > p: + s += 1 + return s + +class ContinuousDistribution: + def __init__(self, distribution): + warnings.warn('Continuous distributions are not thoroughly tested.') + for distribution_type, distribution_parameters in distribution.items(): + self.distribution_type = distribution_type + self.distribution_parameters = distribution_parameters diff --git a/simantha/utils.py b/simantha/utils.py index b6af7a6..3ecf9b6 100644 --- a/simantha/utils.py +++ b/simantha/utils.py @@ -1,3 +1,4 @@ +<<<<<<< HEAD import numpy as np def generate_degradation_matrix(q, dim=10): @@ -6,4 +7,23 @@ def generate_degradation_matrix(q, dim=10): degradation_matrix[i, i] = 1 - q degradation_matrix[i, i+1] = q - return degradation_matrix \ No newline at end of file + return degradation_matrix +======= +def generate_degradation_matrix(p, h_max): + # Returns an upper bidiagonal degradation matrix with probability p of degrading at + # each time step. + degradation_matrix = [] + for h in range(h_max): + transitions = [0] * (h_max + 1) + transitions[h] = 1 - p + transitions[h+1] = p + degradation_matrix.append(transitions) + degradation_matrix.append([0]*h_max + [1]) + return degradation_matrix + +# Time constants (in minutes) +DAY = 24 * 60 +WEEK = 7 * DAY +MONTH = 30 * DAY +YEAR = 365 * DAY +>>>>>>> dev diff --git a/test.py b/test.py new file mode 100644 index 0000000..f86861b --- /dev/null +++ b/test.py @@ -0,0 +1,221 @@ +import random +import unittest + +import scipy.stats + +from simantha import Source, Machine, Buffer, Sink, System +import simantha.simulation +import simantha.utils + +# Degradation transition matrix used for all tests where applicable +degradation_matrix = [ + [0.9, 0.1, 0, 0, 0, 0 ], + [0, 0.9, 0.1, 0, 0, 0 ], + [0, 0, 0.9, 0.1, 0, 0 ], + [0, 0, 0, 0.9, 0.1, 0 ], + [0, 0, 0, 0, 0.9, 0.1], + [0, 0, 0, 0, 0, 1 ] + ] + +class UtilsTests(unittest.TestCase): + """Test Simantha utility functons. + """ + def test_degradation_matrix(self): + degradation_matrix = simantha.utils.generate_degradation_matrix(p=0.1, h_max=5) + target_matrix = [ + [0.9, 0.1, 0, 0, 0, 0 ], + [0, 0.9, 0.1, 0, 0, 0 ], + [0, 0, 0.9, 0.1, 0, 0 ], + [0, 0, 0, 0.9, 0.1, 0 ], + [0, 0, 0, 0, 0.9, 0.1], + [0, 0, 0, 0, 0, 1 ] + ] + #self.assertEqual(degradation_matrix, target_matrix) + self.assertEqual( + simantha.utils.generate_degradation_matrix(p=0.1, h_max=5), + degradation_matrix + ) + +class SimulationTests(unittest.TestCase): + """Tests for the underlying simulation engine. + """ + def test_event_scheduling(self): + # Test scheduling a new event in an empty environment + def dummy_action(): + pass + + env = simantha.simulation.Environment() + new_event = {'time': 0, 'location': None, 'action': dummy_action} + env.schedule_event(**new_event) + + self.assertEqual(len(env.events), 1) + + def test_event_insertion(self): + # Test the ordering of newly scheduled events + def dummy_action(): + pass + + env = simantha.simulation.Environment() + env.schedule_event( + time=10, location=None, action=dummy_action, source='last event' + ) + env.schedule_event( + time=0, location=None, action=dummy_action, source='first event' + ) + env.schedule_event( + time=5, location=None, action=dummy_action, source='middle event' + ) + + event_order = [ev.source for ev in env.events] + + self.assertEqual(event_order, ['first event', 'middle event', 'last event']) + + def test_constant_distribution(self): + # Test sampling of a constant value + constant = simantha.simulation.Distribution({'constant': 42}) + + self.assertTrue(constant.sample(), 42) + + def test_uniform_distribution(self): + # Test sampling of a discrete uniform distribution + random.seed(1) + low, high = [10, 60] + uniform = simantha.simulation.Distribution({'uniform': [low, high]}) + rvs = [uniform.sample() for _ in range(1000)] + + # H_0: The sample is drawn from the specified distribution + _, p_value = scipy.stats.kstest(rvs, scipy.stats.randint(low, high+1).cdf) + + self.assertGreater(p_value, 0.05) + + def test_geometric_distribution(self): + # Test sampling of a geometric distribution + success = 1/100 + geometric = simantha.simulation.Distribution({'geometric': success}) + rvs = [geometric.sample() for _ in range(1000)] + + # H_0: The sample is drawn from the specified distribution + _, p_value = scipy.stats.kstest(rvs, scipy.stats.geom(success).cdf) + + self.assertGreater(p_value, 0.05) + + +class SingleMachineDeterministicTests(unittest.TestCase): + """Basic simulation behavior of a single machine.""" + def build_system(self): + source = Source() + M1 = Machine('M1', cycle_time=1) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1]) + + return System(objects=[source, M1, sink]) + + def test_production(self): + """Test production of a single machine. With cycle time 1, the production count + should be equal to the simulation time.""" + system = self.build_system() + system.simulate(simulation_time=1000, verbose=False) + self.assertEqual(system.machines[0].parts_made, 1000) + + def test_production_warm_up(self): + """Test the warm up period for a single machine. No statistics should be + gathered during the warm up, so production should be equal to simulation + time, regardless of warm up duration.""" + system = self.build_system() + system.simulate(warm_up_time=500, simulation_time=500, verbose=False) + self.assertEqual(system.machines[0].parts_made, 500) + + +class SingleMachineStochsticTests(unittest.TestCase): + """Testing simulation randomness for one machine.""" + def build_system(self): + source = Source() + M1 = Machine( + 'M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution={'constant': 10} + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[sink]) + sink.define_routing(upstream=[M1]) + + return System(objects=[source, M1, sink]) + + def test_production(self): + random.seed(1) + system = self.build_system() + + system.simulate(simulation_time=1000, verbose=False) + + # Assert that the number of parts made is at most the simulation time + self.assertLessEqual(system.machines[0].parts_made, 1000) + + +class TwoMachineDeterministicTests(unittest.TestCase): + """Tests for a two-machine one-buffer line. + """ + def build_system(self): + source = Source() + M1 = Machine('M1', cycle_time=1) + B1 = Buffer('B1', capacity=5) + M2 = Machine('M2', cycle_time=1) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + return System(objects=[source, M1, B1, M2, sink]) + + def test_production(self): + system = self.build_system() + system.simulate(simulation_time=1000, verbose=False) + + self.assertEqual(sum([s.level for s in system.sinks]), 999) + + +class TwoMachineStochasticTests(unittest.TestCase): + def build_system(self): + cm = {'constant': 10} + source = Source() + M1 = Machine( + 'M1', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution=cm + ) + B1 = Buffer('B1', capacity=5) + M2 = Machine( + 'M2', + cycle_time=1, + degradation_matrix=degradation_matrix, + cm_distribution=cm + ) + sink = Sink() + + source.define_routing(downstream=[M1]) + M1.define_routing(upstream=[source], downstream=[B1]) + B1.define_routing(upstream=[M1], downstream=[M2]) + M2.define_routing(upstream=[B1], downstream=[sink]) + sink.define_routing(upstream=[M2]) + + return System(objects=[source, M1, B1, M2, sink]) + + def test_production(self): + system = self.build_system() + system.simulate(simulation_time=1000, verbose=False) + + self.assertLessEqual(system.machines[-1].parts_made, 1000) + + +if __name__ == '__main__': + random.seed(1) + unittest.main()