Monte Carlo Simulation : Underlying distribution parameter changes during simulation run based on conditions

I am trying to simulate running of two manufacturing lines (Line A and Line B).Their time to failure follows Weibull distribution(shape = 0.7, scale = 12 (mins)). Line A and Line B both produces at rate of 100 products/min. From real world scenario, if any of lines fails, I want to increase rate of other line (say – 120 products/min) till the time failed line has not repaired.

Challenge: With increase in rate, chances of failure increases, hence, scale parameter changes (e.g for rate of 120 products/min, scale parameter changes from 12 mins to 10 mins). I want to model this change in distribution parameter in simulation.

Example :

• Line A and Line B starts running at time 0.
• Generate random number using weibull distribution for their time to failure (TTF) – Line A : 15 mins and Line B : 20 mins.
• Line A failed at 15 mins and got repaired at 20 mins. For this 5 mins, Line B ran at increased speed.
• At t= 15 min, generate a random number using weibull (with changed parameter) to get next TTF. The process repeats till the simulation time.

Currently I am using Python Simpy to code the logic but can’t find a way to model this. Any help or reference would be very helpful. Here is my try but i am definitely missing something.

import simpy
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import weibull_min


class Machine(object):

    def __init__(self, env, name,scale_parameter,shape_parameter, mean_repair_time,std_repair_time,increased_rate):
        self.env = env
        self.name = name
        self.scale_parameter = scale_parameter
        self.shape_parameter = shape_parameter
        self.mean_repair_time = mean_repair_time
        self.std_repair_time = std_repair_time
        self.increased_rate = increased_rate
        self.broken = False
        self.processing_list = [0,]
        self.machine_status = 1
        

        self.process = env.process(self.run())
        # Start the failure process
        env.process(self.check_machines_status())
    def run(self):
        """
        Run as long as the simulation runs.
        """
        while True:
            try:
               
                yield self.env.timeout(self.mean_time_to_failure())
                self.processing_list.append(self.env.now)
                print(f'{self.env.now:.2f} {self.name} is in failure.')
                trial_resource.get(1)
                yield self.env.timeout(self.mean_time_to_repair())
                print(f'{self.env.now:.2f} {self.name} is repaired.')
                self.processing_list.append(env.now)
                trial_resource.put(1)
                

            
            except simpy.Interrupt:
                self.machine_status = 0              
                yield self.env.timeout(self.updated_mean_time_to_failure())
                print(f'{self.env.now:.2f} {self.name} is in updated failure.')
                #trial_resource.get(1)
                self.broken = True
                yield self.env.timeout(self.mean_time_to_repair())
                print(f'{self.env.now:.2f} {self.name} is in updated repaired.')
                trial_resource.put(1)
                self.machine_status =1
                
    
    def check_machines_status(self):
        """Periodically check the status of running machines. If any machine fails
        interrupt the process"""
        while True:
            print(self.env.now,trial_resource.level)
            print(self.name)
            
            if trial_resource.level < trial_resource.capacity and self.broken == False and self.machine_status == 1:
                # Only break the machine if it is currently working.
                self.process.interrupt()
                print('Machine running process interrupted %d' % env.now)

            yield env.timeout(1)  


    def mean_time_to_failure(self):
        x = int(weibull_min.rvs(self.shape_parameter, loc=0, scale= self.scale_parameter, size=1).tolist()[0])
        if x == 0:
            x = 1

        return x

    def updated_mean_time_to_failure(self):
        correction_factor = (1-self.increased_rate)/100
        x = int(weibull_min.rvs(self.shape_parameter*correction_factor, loc=0, scale= self.scale_parameter, size=1).tolist()[0])
        if x == 0:
            x = 1

        return x


    def mean_time_to_repair(self):
        x = int(np.random.lognormal(self.mean_repair_time,self.std_repair_time))
        if x ==0:
            x =1
        return x


env = simpy.Environment()
trial_resource = simpy.Container(env,init=3,capacity=3)
machine_1 = Machine(env, 'M1', 12, 0.65, 0.51,1,10)
machine_2 = Machine(env, 'M2', 14, 0.65, 0.51,1,10)
machine_3 = Machine(env, 'M3', 8, 0.65, 0.51,1,10)

env.run(until = 12)
print(machine_1.processing_list)
print(machine_2.processing_list)
print(machine_3.processing_list)

Answer

This solution can handle more then 2 machines. When a machine breaks it sends a message to all of the other machines notifying them of the break down. A machine also sends a notification when it gets fixed. Each machine tracks how many machines are broken and speeds up production if 1 or more machines is broken. Time till failure (ttf) is also sped up. when rates change, the progress of a current in progress part is factored in to calculating the new finish time. The interrupt is used to both trigger a breakdown, but to also reset when a work in progress will finish. also check out the machine shop example in the docs. I ran out of time to work on this so it may not be perfect

"""
Machine shop example where machines change their
production rate and break down rates in response
to other machines breaking down or getting fixed

If one machine breaks down, the remaining machines 
will speed up until all machines are fixed

speeding up the production rate also speeds up the
time till failure ttf

rates can change in the middle of making a part

Programmmer: Michael R. Gibbs
"""

import simpy
import random

class Machine():
    """
    Machine that makes parts.
    Machine breaks down at ramdom time
    The machine has a normal mode, and a fast mode

    Machine speeds up when it receives a message that
    another machine has broken down and slows back down
    when it receives messages that all broken machines are fixed
    
    env:    simulation environment
    id:     name of the machine
    rate:   number of parts per time unit (normal mode)
    fast_rate:  number of parts per time unit (fast mode)
    ttf_gen:    zero param func to generate ttf (normal mode)
    fast_ttf_gen:   zero param func to generate ttf (fast mode)
    """

    def __init__(self, env, id, rate, fast_rate, ttf_gen, fast_ttf_gen):

        self.id = id
        self.env = env
        self.rate = rate
        self.fast_rate = fast_rate
        self.ttf_gen = ttf_gen
        self.fast_ttf_gen = fast_ttf_gen
        self.broken = False
        self.fast_mode = False
        self.remaining_time = 0
        self.current_rate = self.rate

        self.mach_list = []
        self.broke_cnt = 0

        # start the part making, an the count down till failure
        self.make_parts_proc = self.env.process(self.make_parts())
        self.breakdown_proc = self.env.process(self.gen_breakdown())

    def make_parts(self):
        """
        Main loop to manufacture parts

        interupts are used to trigger updates
        when rates change and when breakdowns occure
        """

        while True:

            if self.remaining_time <= 0:
                # starting a new part
                print(f'{self.env.now:.2f}: mach: {self.id} has started a part')
                
                # need to store times so other methods can upate
                # the procesing state
                self.total_part_time = 1 / self.current_rate
                self.remaining_time = self.total_part_time 

            while self.remaining_time > 0:
                # need to loop incase get inturrupted 
                # while in the middle of making a part

                try:
                    self.part_start_time = self.env.now # used to track progress
                    yield self.env.timeout(self.remaining_time)

                    # got here without being interupted, done with part
                    self.remaining_time = 0

                except simpy.Interrupt:
                    # can get inturpted because processing rate changed
                    # or a breakdown has happend
                    # if rate changed, we are using the inturput to update 
                    # the processing timeout
                                            
                    if self.broken:
                        # processing is interuped, fix machine

                        # update processing progress
                        self.remaining_time -= (self.env.now - self.part_start_time)
                        
                        print(f'{self.env.now:.2f}: mach: {self.id} has broken down')
                        
                        # notify other machines that this machine has broke
                        for m in self.mach_list:
                            m.someone_broke(self)

                        # time out for fixing
                        yield self.env.timeout(5)

                        # notify other machines that this machine is fixed
                        for m in self.mach_list:
                            m.someone_fixed(self)
                        
                        print(f'{self.env.now:.2f}: mach: {self.id} has been fixed')

                        # start a new breakdown count down
                        self.breakdown_proc = self.env.process(self.gen_breakdown())
                        self.broken = False

            print(f'{self.env.now:.2f}: mach: {self.id} has finish a part')

    def gen_breakdown(self):
        """
            counts down to failure and uses
            a interrupt to stop processing
            and start repairs

            using class properties instead of local
            variables so other methods can update
            the countdown state of ttf
        """

        if not self.broken:
            # get breakdown based on current fast mode
            if self.fast_mode:
                self.total_ttf = self.fast_ttf_gen()
            else:
                self.total_ttf = self.ttf_gen()

            self.remaining_ttf = self.total_ttf

            while self.remaining_ttf > 0:
                self.ttf_start = self.env.now
                print(f'{self.env.now:.2f}: mach: {self.id} has {self.remaining_ttf} till breakdown')
                try:
                    yield self.env.timeout(self.remaining_ttf)

                    # failure has orrured
                    self.broken = True
                    self.make_parts_proc.interrupt()
                    self.remaining_ttf = 0

                except simpy.Interrupt:
                    # the state has been updated
                    # the timeout has been interupted
                    # so it can be restarted with new state
                    print(f'{self.env.now:.2f}: mach: {self.id} updating ttf {self.remaining_ttf}')
        print(f'{self.env.now:.2f}: mach: {self.id} ttf gen exit')

    def someone_broke(self, mach):
        """
        Another machine is notifing this machine that it has broken
        and this machine needs to change to fast mode, if not already
        in fast mode
        """

        self.broke_cnt += 1
        
        print(f'{self.env.now:.2f}: mach: {self.id} received mess that mach: {mach.id} has broke, broke cnt: {self.broke_cnt}')
    
        if not self.fast_mode:
            self.fast_mode = True

            # update the ttf based on the fast mode ttf func
            # keep the same progress so if we were 90% of of the 
            # old ttf, then set things so we are still 90% of new ttf

            # update with the last bit of progress
            self.remaining_ttf -= self.env.now - self.ttf_start
            per_ttf_left = (self.remaining_ttf/self.total_ttf)

            # update based on fast ttf
            self.total_ttf = self.fast_ttf_gen() 
            self.remaining_ttf = per_ttf_left * self.total_ttf
            if self.remaining_ttf <= 0:
                # special case when notification happens at same time as breakdown
                self.remaining_ttf = 0.001
                
            # update the part processing state

            # update the last bit of processing progress
            self.remaining_time -= self.env.now - self.part_start_time

            # update state based on new processing fast rate
            # like ttf keep the same perenct of progress
            self.current_rate = self.fast_rate
            old_total = self.total_part_time
            self.total_part_time = 1 / self.current_rate
            per_left = self.remaining_time / old_total
            self.remaining_time = self.total_part_time * per_left

            if not self.broken:
                # if broken nothing to interrupt
                # new states will used when machine
                # is fixed and processing starts up again
                self.breakdown_proc.interrupt()
                self.make_parts_proc.interrupt()


                
    def someone_fixed(self, mach):
        """
        Another machine is notifing this machine that it has been fixed
        and this machine needs to change to normal mode, if there are
        no more broken machines
        """

        self.broke_cnt -= 1

        print(f'{self.env.now:.2f}: mach: {self.id} received mess that mach: {mach.id} is fixed, broke cnt: {self.broke_cnt}')
        
        # only change if all machines are fixed and broke cnt is 0
        if self.broke_cnt <= 0:
            self.broke_cnt = 0

            if self.fast_mode:
                self.fast_mode = False

            # update the ttf based on the normal mode ttf func
            # keep the same progress so if we were 90% of of the 
            # old ttf, then set things so we are still 90% of new ttf

            # update with the last bit of progress
            self.remaining_ttf -= self.env.now - self.ttf_start
            per_ttf_left = (self.remaining_ttf/self.total_ttf)

            self.total_ttf = self.ttf_gen()
            self.remaining_ttf = per_ttf_left * self.total_ttf

            if self.remaining_ttf <= 0:
                # special case when notifcation happens at breakdown time
                self.remaining_ttf = 0.001
                
            # update state based on new processing normal rate
            # like ttf keep the same perenct of progress
            self.remaining_time -= self.env.now - self.part_start_time
            self.current_rate = self.rate
            old_total = self.total_part_time
            self.total_part_time = 1 / self.current_rate
            per_left = self.remaining_time / old_total
            self.remaining_time = self.total_part_time * per_left

            if not self.broken:
                # if broken nothing to interrupt
                # new states will be used when machine
                # is fixed and processing starts up again
               self.breakdown_proc.interrupt()
               self.make_parts_proc.interrupt()



    def set_mach_list(self, mach_list):
        """
        set the list of machines to be notified if this machine
        breaks down, or is fixed
        """
        self.mach_list = mach_list 

# ttf generator function
#
# by wrapping a dist in a lambda I 
# created create a 0 param function
# that I can pass to the Machine class
# To change the dist I just need to 
# update the lambda, no hard codeing 
# of dist func parameters in the Machine
# class code
ttf = lambda : random.randint(8,10)
fast_ttf = lambda: random.randint(5,7)

# create sim
env = simpy.Environment()

mach_list = []
machines_cnt = 2 # can be more then 2 
for i in range(1, machines_cnt + 1):
    m = Machine(env, i, 5, 8, ttf, fast_ttf)
    mach_list.append(m)

# build list of machines to notify
# when a machine breaks or gets fixed
for m in mach_list:
    # filter out the current machine
    # don't want to send to self
    other_m = mach_list.copy()
    other_m.remove(m)

    m.set_mach_list(other_m)

env.run(until = 50)

print("end of simulation")