Monte Carlo Analysis of a Discrete Event Simulation of a non deterministic manufacturing process with multi resources using SimPy

import random
import simpy
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import schemdraw
from schemdraw.flow import *

def print_stats(res):
    """Simpy native function to allow looking into queues and resource utilization while the simulation runs"""
    print(f'{res.count} of {res.capacity} slots are allocated.')
    print(f'  Users: {res.users}')
    print(f'  Queued events: {res.queue}')
    print(f'\n')
    
def all_processes(env, batch_i, verbose=True):
    """Modeling the whole manufacturing process that is defined in the process_activities dictionary by scheduling, processing 
    and releasing batches from/to resources"""
    
    for p_activity in process_activities:

        if verbose==True:
            print("> New production batch " + str(batch_i) + " scheduled for " + 
                  p_activity.get('Name') + "@ current time: " + str(np.round(env.now, 2)))

        resource = p_activity.get('Resource')
        request_time = env.now
       
        resource_request = resource.request()
        
        batch_activity_df.loc[(batch_i, p_activity.get('Name')), 'Scheduled'] = env.now
        
        yield resource_request
        
        if verbose==True:
            print("> Start processing batch " + str(batch_i) + " scheduled for " + 
                  p_activity.get('Name') + " @current time: " + str(np.round(env.now, 2)))
        
        batch_activity_df.loc[(batch_i, p_activity.get('Name')), 'Processing'] = env.now
        
        duration = activity(gen_type=p_activity.get('gen_type'), 
                                       gauss_mean=p_activity.get('Average_Duration'), 
                                       gauss_std=p_activity.get('Standard_Deviation'),
                                       unif_start=p_activity.get('Start_Range'), 
                                       unif_end=p_activity.get('End_Range')
                           )

        if (env.now > request_time):
            if verbose==True:
                print("batch {} had to wait " + str(np.round(env.now - request_time,2)) + " mins on activity {}".format(batch_i, 
                                                                            p_activity.get('Name')))
            batch_activity_df.loc[(batch_i, p_activity.get('Name')), 'Wait_time'] = env.now - request_time
        try:    
            yield env.timeout(duration)
        except: # negative durations may occur due to erroneous distributions
            print(duration, p_activity.get('Name'))
            
        resource.release(resource_request)
        batch_activity_df.loc[(batch_i, p_activity.get('Name')), 'Releasing'] = env.now
        if verbose==True:
            print("> Finished Batch " + str(batch_i) + " processed on activity " + 
                  p_activity.get('Name') + " in " + str(np.round(duration,2)) 
                  + " minutes => current time: " + str(np.round(env.now, 2)))

def activity(gen_type=None, gauss_mean=None, gauss_std=None, unif_start=None, unif_end=None):
    """Minor process activity property helper function"""
    if gen_type == "gauss":
        duration = random.gauss(gauss_mean, gauss_std)
    elif gen_type == "uniform":
        duration = random.uniform(unif_start, unif_end)
    return duration


def batch_generator(env, verbose=True):
    """Manufacturing Batch generation with batch seeded Poisson distributed arrival times"""
    batch_i = 0
    
    while True:
        
        inter_arrival_time = np.random.default_rng(batch_i).poisson(45, 1)
        
        yield env.timeout(inter_arrival_time)
        if verbose==True:
            print('new batch generated: ', batch_i)
        env.process(all_processes(env, batch_i, verbose=verbose))
        batch_i += 1

#rng = np.random.default_rng()
#plt.hist(rng.poisson(45, 1))

global performance_series
performance_df = pd.DataFrame(columns = ['num_complete_batches', 'capacity_vector'])


global utilization_df
utilization_df = pd.DataFrame(columns = ['Raw Material Goods receipt',
                                         'Pre-Cutting',
                                         'Polishing',
                                         'Adding Component 1',
                                         'Adding Component 2',
                                         'Adding Component 3',
                                         'Final assembly',
                                         'Quality control',
                                         'Packaging', 'capacity_vector'])
global error_check
error_check = pd.DataFrame(columns = ['Raw Material Goods receipt',
                                         'Pre-Cutting',
                                         'Polishing',
                                         'Adding Component 1',
                                         'Adding Component 2',
                                         'Adding Component 3',
                                         'Final assembly',
                                         'Quality control',
                                         'Packaging'])

global batch_cycle_time_df
batch_cycle_time_df = pd.DataFrame(columns = ['run_id', 'batch_id', 'process_cycle_time', 'capacity_vector'])

global wait_time_df
wait_time_df = pd.DataFrame(columns = ['activity', 'mean', 'run_id', 'capacity_vector'])

def run_monte_carlo(capacity_vector):
    """Conducting Monte Carlo simulation runs with a given capacity vector of the defined process in the 
    process_activities dict implemented with Simpy Simulation Framework"""    
    run_i = 0
    runs = 50
   
    while run_i < runs:
    #for run_i in range(runs):
       
        random.seed(run_i)
                
        env = simpy.Environment()
        
        global process_activities
        process_activities = [{'Name': 'Raw Material Goods receipt',
                               'gen_type': "gauss",
          'Average_Duration': 50,
          'Standard_Deviation': 5,
          'Resource': simpy.Resource(env, capacity=capacity_vector[0])},

         {'Name': 'Pre-Cutting', 'gen_type': "gauss",
          'Average_Duration': 20,
          'Standard_Deviation': 3,
          'Resource': simpy.Resource(env, capacity=capacity_vector[1])},

         {'Name': 'Polishing', 'gen_type': "gauss",
          'Average_Duration': 20,
          'Standard_Deviation': 5,
          'Resource': simpy.Resource(env, capacity=capacity_vector[2])},

         {'Name': 'Adding Component 1', 'gen_type': "gauss",
          'Average_Duration': 5,
          'Standard_Deviation': 1,
          'Resource': simpy.Resource(env, capacity=capacity_vector[3])},

         {'Name': 'Adding Component 2', 'gen_type': "gauss",
          'Average_Duration': 10,
          'Standard_Deviation': 2,
          'Resource': simpy.Resource(env, capacity=capacity_vector[4])},

         {'Name': 'Adding Component 3', 'gen_type': "gauss",
          'Average_Duration': 8,
          'Standard_Deviation': 1,
          'Resource': simpy.Resource(env, capacity=capacity_vector[5])},

         {'Name': 'Final assembly', 'gen_type': "gauss",
          'Average_Duration': 90,
          'Standard_Deviation': 10,
          'Resource': simpy.Resource(env, capacity=capacity_vector[6])},

         {'Name': 'Quality control', 'gen_type': "gauss",
          'Average_Duration': 10,
          'Standard_Deviation': 1,
          'Resource': simpy.Resource(env, capacity=capacity_vector[7])},

         {'Name': 'Packaging', 'gen_type': "uniform",
          'Start_Range': 5,
          'End_Range': 10,
          'Resource': simpy.Resource(env, capacity=capacity_vector[8])}
        ]

        global batch_activity_df
        batch_activity_df = pd.DataFrame(columns = ['batchid', 'activity', 'Scheduled', 
                                            'Processing', 'Releasing', 'Wait_time']).set_index(['batchid', 'activity'])
        
        #six day week
        #duration in mins
        sim_duration = 6*24*60

        b = env.process(batch_generator(env, verbose=False))
        c = env.run(until=sim_duration)
        
        process_stats(batch_activity_df, capacity_vector, sim_duration, run_i)
        
        # wait times stats across activities
        #wait_time_stats_list.append(pd.pivot_table(data=batch_activity_df, index=['activity'], 
        #       aggfunc=['min', 'max', 'mean', 'median', 'sum'], values='Wait_time',
        #       fill_value=0
        #      ))
        
        print('monte carlo run: ' + str(run_i))
        run_i += 1

        #del env
        #del b
        #del c

def process_stats(batch_activity_df, capacity_vector, sim_duration, run_i):
    """Processing stats for each sim run for later analysis"""
    
    ba_df = batch_activity_df.copy(deep=True)
    ba_df['Releasing'].fillna(sim_duration, inplace=True)
    ba_df['Processing_duration'] = ba_df.Releasing - ba_df.Processing
    process_durations_series = ba_df.groupby('activity')['Processing_duration'].sum()
    
    #check
    caps = np.multiply(capacity_vector, sim_duration)
    #print(capacity_vector)
    #print(process_durations_series)
    
    # activity durations
    process_durations_df = pd.DataFrame(data=[process_durations_series, 
                                          pd.Series(np.multiply(capacity_vector, sim_duration),
                                                    name='Capacity',
                                                    index=list((object['Name'] for object in process_activities)))],
                                   ).T
    
    if len(process_durations_df.Capacity) == len(process_durations_series):
        process_durations_df['Utilization_perc'] = np.round(process_durations_df.Processing_duration.astype(float).divide(process_durations_df.Capacity.astype(float))*100,2)
    else:
        print('lengths: ', len(process_durations_df.Capacity), len(process_durations_series))
    
    #error check
    error_check.loc[run_i] = process_durations_series
    
    # resource utilization
    #print(type(process_durations_df['Utilization_perc']))
    util = process_durations_df['Utilization_perc'].copy()
    util['capacity_vector'] = str(capacity_vector)
    #display(util)
    utilization_df.loc[len(utilization_df)] = util
       
    # performance
    # count complete batches
    # must use untouched batch_activity_df!
    # max batch_id + 1 for the count
    max_compl_batch = batch_activity_df.loc[(slice(None), 'Packaging'), ['Scheduled','Processing','Releasing']].sort_values('batchid', ascending=False).dropna(how='any').index.max()[0]
    
    performance_df.loc[len(performance_df)] = [max_compl_batch+1, str(capacity_vector)]
    
    # wait_time
   
    wt_df = batch_activity_df.copy().reset_index()
    compl_batches_df = wt_df.loc[wt_df['batchid'] <= max_compl_batch]
    #compl_batches_df = batch_activity_df.query('0 <= batchid <='+ str(max_compl_batch))
    run_activity_mean_wait_time_df = pd.pivot_table(data=compl_batches_df, index=['activity'], 
                                       aggfunc=['mean'], values='Wait_time',
                                       fill_value=0)
    run_activity_mean_wait_time_df['run_id'] = run_i
    run_activity_mean_wait_time_df['capacity_vector'] = str(capacity_vector)
    run_activity_mean_wait_time_df.columns = run_activity_mean_wait_time_df.columns.droplevel(1)
    run_activity_mean_wait_time_df.reset_index(inplace=True)
    #print(run_activity_mean_wait_time_df)
    global wait_time_df
    wait_time_df = pd.concat([wait_time_df, run_activity_mean_wait_time_df])
    
    # process cycle time
    # must use untouched batch_activity_df!
    for batch_id in range(max_compl_batch+1):
        start = float(batch_activity_df.loc[(batch_id, 'Raw Material Goods receipt'), ['Processing']])
        end = float(batch_activity_df.loc[(batch_id, 'Packaging'), ['Releasing']])
        process_cycle_time = end - start
        batch_cycle_time_df.loc[len(batch_cycle_time_df)] = [run_i, batch_id, float(process_cycle_time), str(capacity_vector)]

#run_monte_carlo([1,1,1,1,1,1,2,1,1])

capacity_vectors = [
                    [1,1,1,1,1,1,2,1,1],
                    [2,1,1,1,1,1,2,1,1],
                    [2,1,1,1,1,1,3,1,1]
                   ]

# Call the run_monte_carlo function across all capacities scenarios in form of capacity vectors
for capacity_vector in capacity_vectors:
    run_monte_carlo(capacity_vector)

performance_df

set1 = ['Adding Component 1', 'Adding Component 2', 'Adding Component 3', 'Quality control', 'Packaging', 'capacity_vector']
set2 = ['Raw Material Goods receipt', 'Pre-Cutting','Polishing']
set3 = ['Final assembly', 'capacity_vector']

utilization_df2 = utilization_df[utilization_df['capacity_vector'] == '[2, 1, 1, 1, 1, 1, 2, 1, 1]'].copy()

fig, axs = plt.subplots(3, figsize=(7, 11))
fig.suptitle('Resource utilization distributions across process activities \n for Scenario 1 ([2, 1, 1, 1, 1, 1, 2, 1, 1])) in Percent')

sns.boxplot(data=utilization_df2[set1], orient='h', ax=axs[0])
sns.stripplot(data=utilization_df2[set1], orient='h',size=3, palette='dark:white', linewidth=0.5, ax=axs[0], #hue='capacity_vector'
             )

sns.boxplot(data=utilization_df2[set2], orient='h', ax=axs[1])
sns.stripplot(data=utilization_df2[set2], orient='h',size=4,  palette='dark:white', linewidth=0.5, ax=axs[1])

sns.boxplot(data=utilization_df2[set3], orient='h')
sns.stripplot(data=utilization_df2[set3], orient='h',size=5,  palette='dark:white', linewidth=0.5, ax=axs[2])

plt.show()

sns.displot(data=utilization_df2[set3], kind="ecdf")
plt.title("Empirical cumulative distribution function of 'Final assembly' activity resource utilization")
plt.xlabel('Resource Utilization in percent')
plt.show()

def pdf(performance_df):
    """Analyse and visualize performance stats ex post across runs"""
    pivot = pd.pivot_table(data=performance_df, index=['capacity_vector'], 
               aggfunc=[ 'count', 'mean', 'sem']
                        ,# values='num_complete_batches',
               #fill_value=0
              ).droplevel(1, axis=1)
    
    df = pivot.copy()
   
    from scipy import stats
    confidence = 0.90
    dof = len(df)-1  
    critical_val_z = np.abs(stats.t.ppf((1-confidence)/2,dof))
    
    df['ci_upper_limit'] = df['mean'] + critical_val_z * df['sem']
    df['ci_lower_limit'] = df['mean'] - critical_val_z * df['sem']
    df['performance_index_percent'] = np.round(df['mean'] / df['mean'].min() *100, 2)
    
    title = 'Performance Comparison for two new capacity configurations vs. baseline (completed batches per work week)'
    print(title)
    display(df)
    
    sns.set(rc={'figure.figsize':(6,9)})
    sns.set_style('dark')
    sns.barplot(data=performance_df, y='num_complete_batches', x='capacity_vector', errorbar=('ci', 90))
    plt.title(title)
    plt.ylabel('Average Number of completed batches')
    plt.xlabel('Capacity configurations as vectors')
    plt.show()

pdf(performance_df)

try:
    wait_time_df['mean'] = wait_time_df['mean'].str[0]
except:
    pass

wait_time_df.groupby('capacity_vector').mean()

pd.pivot_table(data=batch_cycle_time_df, index=['capacity_vector'], 
               aggfunc=['min', 'max', np.mean, 'median', 'sum'], values='process_cycle_time')

batch_cycle_time_df.sort_values('process_cycle_time', ascending=False)

# appendix
# just visualizing the process

fontsize=11
boxwidth=2.5
counter=0
schemdraw.theme(theme='dark')

with schemdraw.Drawing() as d:
    d+= Start().label('Start', fontsize=fontsize)
    d+= Arrow().right(d.unit)
      
    for p_activity in process_activities:
        name = p_activity.get('Name')
        
        d+= Box(w = boxwidth).label(p_activity.get('Name').replace(" ","\n"), fontsize=fontsize)
        if counter < 3:
            d+= Arrow().right(d.unit/2)
        elif counter == 3:
            d+= Arrow().down(d.unit/2)
        elif counter == 4:
            d+= Arrow().left(d.unit)
        else:
            d+= Arrow().left(d.unit/2)
        counter += 1
    
    #End program
    d+= (end := Ellipse().label("End", fontsize=fontsize))
d.save('process_schemadraw.svg')