<jupyter_start><jupyter_text>12.5.2 Continuous monitoring and finite discrete state spaceScripts for Monte Carlo simulations and numerical computations corresponding to Markov graphes on Figure 12.10 are presented. The notebook is structured as follows:+ Generic functions for Monte Carlo simulations and numerical computation for cases 1-5 on Figure 12.10+ Case 1: use of the generic functions for getting state probabilities in case 1 - No preventive maintenance+ Case 2: use of the generic functions for getting state probabilities in case 2 - Perfect preventive maintenance+ Case 3: use of the generic functions for getting state probabilities in case 3 - Unperfect preventive maintenance+ Case 4: use of the generic functions for getting state probabilities in case 4 - Two levels preventive maintenance+ Case 5: use of the generic functions for getting state probabilities in case 5 - Two steps preventive maintenance+ Compare cases     Generic functions<jupyter_code>import numpy as np
import pandas as pd
from scipy.linalg import expm

%matplotlib notebook
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Functions for Monte Carlo simulation<jupyter_code>def RandExp(rate):
    '''Draw random samples from an exponential distribution.
    If rate is null, output is set to infinity.'''
    if rate==0:
        return np.inf
    else:
        return np.random.exponential(scale=1/rate)
def UpdateVar(horizonTime, currentTime, currentState, dT, nextState, 
              timeState, transMatrix):
    '''Update system state and time spent in the different states and achieved transitions.'''
    currentTime = currentTime+dT
    flagSim = True
    if currentTime>=horizonTime:
        flagSim = False
        dT = dT-(currentTime-horizonTime)
        currentTime = currentTime+dT
    timeState[currentState] = timeState[currentState]+dT
    transMatrix[currentState,nextState] = transMatrix[currentState,nextState]+1
    currentState = nextState
    return (flagSim, currentTime, currentState, timeState, transMatrix)
def OneHistory(para):
    '''Simulate one history of the system until simulation time horizon is reached.
    '''
    # Intial conditions
    currentState = 3
    currentTime = 0
    flagSim = True
    # Counters initialization
    nState = 4
    timeState = np.zeros(nState)
    transMatrix = np.zeros((nState, nState))
    while flagSim:
        if (currentState==3):
            # Evaluate concurrent transitions and determine the actual one
            dT = RandExp(rate=para['lambda32'])
            nextState = 2
        elif (currentState==2):
            # Evaluate concurrent transitions and determine the actual one
            dT23 = RandExp(rate=para['mu23'])
            dT21 = RandExp(rate=para['lambda21'])
            dT = np.min([dT23, dT21])
            if (dT==dT23):
                nextState = 3
            elif (dT==dT21):
                nextState = 1
            else:
                raise ValueError('Approximation error')
        elif (currentState==1):
            # Evaluate concurrent transitions and determine the actual one
            dT12 = RandExp(rate=para['mu12'])
            dT13 = RandExp(rate=para['mu13'])
            dT10 = RandExp(rate=para['lambda10'])
            dT = np.min([dT12, dT13, dT10])
            if (dT==dT12):
                nextState = 2
            elif (dT==dT13):
                nextState = 3
            elif (dT==dT10):
                nextState = 0
            else:
                raise ValueError('Approximation error')
        elif (currentState==0):
            # Evaluate concurrent transitions and determine the actual one
            dT = RandExp(rate=para['mu03'])
            nextState = 3
        else:
            raise ValueError('Unknown state')
        # Update variables
        (flagSim, currentTime, currentState, timeState, transMatrix) = UpdateVar(
            para['horizon'], currentTime, currentState, dT,  nextState, 
            timeState, transMatrix)
    return (timeState, transMatrix)
def NHistories(nbN, para):
    '''Simulate several histories and aggregate time spent in the different states and achieved
    transitions (Monté-carlo approach).
    '''
    # Simulate a first history to assess output shapes
    (timeStateOne, transMatrixOne) = OneHistory(para)
    timeState = np.empty((nbN, timeStateOne.shape[0]))
    transMatrix = np.empty((nbN, transMatrixOne.shape[0], transMatrixOne.shape[1]))
    # Main loop
    for id in range(nbN):
        (timeState[id,:], transMatrix[id,:,:]) = OneHistory(para)
    return (timeState, transMatrix)<jupyter_output><empty_output><jupyter_text>Functions for numerical computation<jupyter_code>def NumComp(nbSamp, para):
    mA = np.array([
        [np.nan, 0, 0, para['mu03']],
        [para['lambda10'], np.nan, para['mu12'], para['mu13']],
        [0, para['lambda21'], np.nan, para['mu23']],
        [0, 0, para['lambda32'], np.nan]])
    np.fill_diagonal(mA, -np.nansum(mA, axis=1))
    vP0 = np.array([0, 0, 0, 1])
    vPmean = np.zeros((4, ))
    for t in np.linspace(0, para['horizon'], nbSamp):
        vPend = vP0@expm(mA*t)
        vPmean = vPmean+vPend
    vPmean = vPmean/nbSamp
    return vPmean, vPend
def GetResults(case, timeState, transMatrix, probStateMean, probStateEnd):
    '''Get and display simulations results'''
    # Probability of each states
    res = {'Case': case}
    for id in range(4):
        res['State{:d}'.format(id)] = [np.mean(timeState[:, id]/np.sum(timeState[:, :], axis=1)), 
                                       probStateMean[id], probStateEnd[id]]
        print('State {:d} -> Probability: {:.5f} (Sim.) {:.5f} (Cal1.) {:.5f} (Cal2.)'.format(
            id, res['State{:d}'.format(id)][0], res['State{:d}'.format(id)][1], 
            res['State{:d}'.format(id)][2], axis=1))
    # Average simulated time
    res['AveSimTime'] = np.mean(np.sum(timeState[:, :], axis=1))
    print('Average simulated time (h): {:.2f}'.format(res['AveSimTime']))
    # Average operating time
    res['AveOpeTime'] = np.mean(np.sum(timeState[:, 1:], axis=1))
    print('Average operating time (h): {:.2f}'.format(res['AveOpeTime']))
    # Average repairing time
    res['aveRepTime'] = np.mean(timeState[:, 0])
    print('Average repairing time (h): {:.2f}'.format(res['aveRepTime']))
    # Average number of corrective maintenance
    res['aveCorMaint'] = np.mean(transMatrix[:, 0, 3])
    print('Average number of corrective maintenance: {:.2f}'.format(res['aveCorMaint']))
    # Average number of preventive maintenance
    res['avePreMaint'] = np.mean(transMatrix[:, 1, 2])+np.mean(transMatrix[:, 1, 3]+np.mean(transMatrix[:, 2, 3])) 
    print('Average number of preventive maintenance: {:.2f}'.format(res['avePreMaint']))    
    return res<jupyter_output><empty_output><jupyter_text>CASE 1: NO PREVENTIVE MAINTENANCE<jupyter_code># Parameters
nbHist = int(1e5)
nbSamp = 1000
para = {
    'horizon' : 10*365*24, # Horizon time at 10 years, unit is hour
    'lambda32': 3e-5,
    'lambda21': 3e-5,
    'lambda10': 3e-5,
    'mu03': 1e-3,
    'mu12': 0,
    'mu13': 0,
    'mu23': 0}
# Simulation
(timeState, transMatrix) = NHistories(nbHist, para)
# Numerical computation
(probStateMean, probStateEnd) = NumComp(nbSamp, para)
# Get simulations results
print('Case 1: No preventive maintenance\n')
res1 = GetResults(1, timeState, transMatrix, probStateMean, probStateEnd)<jupyter_output>Case 1: No preventive maintenance

State 0 -> Probability: 0.00604 (Sim.) 0.00604 (Cal1.) 0.00978 (Cal2.)
State 1 -> Probability: 0.20549 (Sim.) 0.20501 (Cal1.) 0.32647 (Cal2.)
State 2 -> Probability: 0.32916 (Sim.) 0.32913 (Cal1.) 0.34252 (Cal2.)
State 3 -> Probability: 0.45932 (Sim.) 0.45982 (Cal1.) 0.32123 (Cal2.)
Average simulated time (h): 87600.00
Average operating time (h): 87070.92
Average repairing time (h): 529.08
Average number of corrective maintenance: 0.54
Average number of preventive maintenance: 0.00<jupyter_text>CASE 2: PERFECT PREVENTIVE MAINTENANCE<jupyter_code># Parameters
nbHist = int(1e5)
nbSamp = 1000
para = {
    'horizon' : 10*365*24,
    'lambda32': 3e-5,
    'lambda21': 3e-5,
    'lambda10': 3e-5,
    'mu03': 1e-3,
    'mu12': 0,
    'mu13': 1e-3,
    'mu23': 0}
# Simulation
(timeState, transMatrix) = NHistories(nbHist, para)
# Numerical computation
(probStateMean, probStateEnd) = NumComp(nbSamp, para)
# Get simulations results
print('Case 2: Perfect preventive maintenance\n')
res2 = GetResults(2, timeState, transMatrix, probStateMean, probStateEnd)<jupyter_output>Case 2: Perfect preventive maintenance

State 0 -> Probability: 0.00034 (Sim.) 0.00034 (Cal1.) 0.00043 (Cal2.)
State 1 -> Probability: 0.01151 (Sim.) 0.01151 (Cal1.) 0.01427 (Cal2.)
State 2 -> Probability: 0.40171 (Sim.) 0.40064 (Cal1.) 0.49025 (Cal2.)
State 3 -> Probability: 0.58644 (Sim.) 0.58751 (Cal1.) 0.49505 (Cal2.)
Average simulated time (h): 87600.00
Average operating time (h): 87570.38
Average repairing time (h): 29.62
Average number of corrective maintenance: 0.03
Average number of preventive maintenance: 1.02<jupyter_text>CASE 3: UNPERFECT PREVENTIVE MAINTENANCE<jupyter_code># Parameters
nbHist = int(1e5)
nbSamp = 1000
para = {
    'horizon' : 10*365*24,
    'lambda32': 3e-5,
    'lambda21': 3e-5,
    'lambda10': 3e-5,
    'mu03': 1e-3,
    'mu12': 1e-3,
    'mu13': 0,
    'mu23': 0}
# Simulation
(timeState, transMatrix) = NHistories(nbHist, para)
# Numerical computation
(probStateMean, probStateEnd) = NumComp(nbSamp, para)
# Get simulations results
print('Case 3: Unperfect preventive maintenance\n')
res3 = GetResults(3, timeState, transMatrix, probStateMean, probStateEnd)<jupyter_output>Case 3: Unperfect preventive maintenance

State 0 -> Probability: 0.00053 (Sim.) 0.00052 (Cal1.) 0.00077 (Cal2.)
State 1 -> Probability: 0.01780 (Sim.) 0.01774 (Cal1.) 0.02561 (Cal2.)
State 2 -> Probability: 0.61867 (Sim.) 0.61875 (Cal1.) 0.88122 (Cal2.)
State 3 -> Probability: 0.36300 (Sim.) 0.36299 (Cal1.) 0.09241 (Cal2.)
Average simulated time (h): 87600.00
Average operating time (h): 87553.46
Average repairing time (h): 46.54
Average number of corrective maintenance: 0.05
Average number of preventive maintenance: 1.58<jupyter_text>CASE 4: TWO LEVELS PREVENTIVE MAINTENANCE<jupyter_code># Parameters
nbHist = int(1e5)
nbSamp = 1000
para = {
    'horizon' : 10*365*24,
    'lambda32': 3e-5,
    'lambda21': 3e-5,
    'lambda10': 3e-5,
    'mu03': 1e-3,
    'mu12': 0,
    'mu13': 1e-3,
    'mu23': 1e-3}
# Simulation
(timeState, transMatrix) = NHistories(nbHist, para)
# Numerical computation
(probStateMean, probStateEnd) = NumComp(nbSamp, para)
# Get simulations results
print('Case 4: Two levels preventive maintenance\n')
res4 = GetResults(4, timeState, transMatrix, probStateMean, probStateEnd)<jupyter_output>Case 4: Two levels preventive maintenance

State 0 -> Probability: 0.00002 (Sim.) 0.00002 (Cal1.) 0.00002 (Cal2.)
State 1 -> Probability: 0.00076 (Sim.) 0.00081 (Cal1.) 0.00082 (Cal2.)
State 2 -> Probability: 0.02792 (Sim.) 0.02796 (Cal1.) 0.02828 (Cal2.)
State 3 -> Probability: 0.97130 (Sim.) 0.97121 (Cal1.) 0.97087 (Cal2.)
Average simulated time (h): 87600.00
Average operating time (h): 87597.84
Average repairing time (h): 2.16
Average number of corrective maintenance: 0.00
Average number of preventive maintenance: 2.55<jupyter_text>CASE 5: TWO STEPS PREVENTIVE MAINTENANCE<jupyter_code># Parameters
nbHist = int(1e5)
nbSamp = 1000
para = {
    'horizon' : 10*365*24,
    'lambda32': 3e-5,
    'lambda21': 3e-5,
    'lambda10': 3e-5,
    'mu03': 1e-3,
    'mu12': 1e-3,
    'mu13': 0,
    'mu23': 1e-3}
# Simulation
(timeState, transMatrix) = NHistories(nbHist, para)
# Numerical computation
(probStateMean, probStateEnd) = NumComp(nbSamp, para)
# Get simulations results
print('Case 5: Two steps preventive maintenance\n')
res5 = GetResults(5, timeState, transMatrix, probStateMean, probStateEnd)<jupyter_output>Case 5: Two steps preventive maintenance

State 0 -> Probability: 0.00002 (Sim.) 0.00002 (Cal1.) 0.00003 (Cal2.)
State 1 -> Probability: 0.00081 (Sim.) 0.00083 (Cal1.) 0.00085 (Cal2.)
State 2 -> Probability: 0.02893 (Sim.) 0.02873 (Cal1.) 0.02908 (Cal2.)
State 3 -> Probability: 0.97024 (Sim.) 0.97042 (Cal1.) 0.97005 (Cal2.)
Average simulated time (h): 87600.00
Average operating time (h): 87597.90
Average repairing time (h): 2.10
Average number of corrective maintenance: 0.00
Average number of preventive maintenance: 2.63<jupyter_text>Compare cases<jupyter_code>df = pd.DataFrame([res1, res2, res3, res4, res5]).set_index('Case')
df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13093720114759325, 'fraction_numerical': 0.07260121134842205, 'mean_word_length': 3.640902366863905, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '638440', 'license': 'MIT', 'n_tokens_mistral': 4818, 'n_tokens_neox': 4137, 'n_words': 1194, 'path': '12_PreventiveMaintenance/12_05_ConditionBasedMaintenance/12_05_02_ContinuousMonitoringAndFiniteDiscreteStateSpace1.ipynb', 'repo_name': 'RausandBarros/ReliabilityBookScripts', 'size': '12548'}