Stochastic

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
import scipy as sp

def reliability(beta0_b, d_b, L_b, Thoriz):
    # --- Reliability estimation
    # beta_bt: reliability curve over time for bridge b
    # beta0_b: Initial reliability of bridge b at t=0
    # d_b, L_b: when the degradation starts and reaches beta=0 for each bridge
    # Thoriz: maximum time studied
    nb = len(beta0_b)
    beta_bt = np.zeros((nb, Thoriz))
    
    for b in range(nb):
        #horziontal branch
        timesteps = int(np.floor(d_b[b]))
        beta_bt[b, :timesteps] = beta0_b[b]
        
        #degradation till L_b
        servicelife = int(np.floor(L_b[b]))
        t = np.arange(timesteps + 1, servicelife + 1)
        beta_bt[b, timesteps:servicelife] = beta0_b[b] / (L_b[b] - d_b[b]) * (L_b[b] - t)
    
    return beta_bt

def MaintenanceApplication(beta_bt, X_bt, Thoriz):
    # --- Given a maintenance strategy given by X_bt, where 1 means that bridge
    # b undergoes maintenance at time t, this function gives the updated
    # reliability profile
    updatedbeta_bt = beta_bt.copy()
    InterventionTimes = np.where(np.sum(X_bt, axis=0) > 0)[0]
    
    for t in InterventionTimes:
        intervbrid = np.where(X_bt[:, t] == 1)[0]
        timecells = Thoriz - t - 1
        updatedbeta_bt[intervbrid, t+1:Thoriz] = beta_bt[intervbrid, :timecells]
    
    return updatedbeta_bt

def Frequency2Schedule(y_b, nb, nt, Thoriz):
    # Transforms a vector indicating the frequency to a schedule matrix [0/1]
    freq_b = y_b.reshape((nt, nb)).T
    cumFreq_b = cumulateFreq(freq_b)
    X_bt = np.zeros((nb, Thoriz))
    
    for b in range(nb):
        idx = cumFreq_b[b, cumFreq_b[b, :] != 0].astype(int)
        idx = idx[idx <= Thoriz]
        if len(idx) > 0:
            X_bt[b, idx - 1] = 1
    
    return X_bt

def cumulateFreq(freq_b):
    # creates the cumulative frequency
    nb, nt = freq_b.shape
    cumFreq_b = np.zeros((nb, nt))
    
    for b in range(nb):
        f_b = freq_b[b, freq_b[b, :] != 0]
        cumFreq_b[b, :len(f_b)] = np.cumsum(f_b)
    
    return cumFreq_b

def CountNinterv(y_b, nt, nb):
    # from vector y_b, it computes the number of interventions per bridge
    y_b01 = y_b != 0
    freq_b = y_b01.reshape((nt, nb)).T
    ninterv_b = np.sum(freq_b, axis=1)
    return ninterv_b

def ProbFulfilment(beta_bti0, X_bt, Thoriz):
    nb, _, nMCS = beta_bti0.shape

    # --- MCS
    Achivement_bi = np.zeros((nb, nMCS))
    for i in range(nMCS):
        beta_bt = MaintenanceApplication(beta_bti0[:, :, i], X_bt, Thoriz)
        failures_b = np.sum(beta_bt == 0, axis=1)
        Achivement_bi[:, i] = (failures_b == 0)
    
    ProbRel_b = np.sum(Achivement_bi, axis=1) / nMCS
    return ProbRel_b

Thoriz = 25
# for each bridge b:  beta0 mu_d(years) sigma_d(years) mu_L(years) sigma_L(years) C(Meuros)
data = np.array([
    [1, 5, 2 / 1.96, 10, 1 / 1.96, 0.6],
    [1, 6, 1 / 1.96, 11, 2 / 1.96, 0.8],
    [1, 7, 1 / 1.96, 15, 5 / 1.96, 1.0],
    [1, 10, 1 / 1.96, 20, 1 / 1.96, 0.7]
])
beta0_b, mud_b, sigd_b, muL_b, sigL_b, C_b = data[:, 0], data[:, 1], data[:, 2], data[:, 3], data[:, 4], data[:, 5]

# Annual resources (M euro)
R = 1.4

# Sample size of the probabilistic analysis
nMCS=10

# --- Initial reliability estimation
nb = len(beta0_b)  # number of bridges

# nMCS cases of d and L
d_bi, L_bi = np.zeros((nb, nMCS)), np.zeros((nb, nMCS))
for b in range(nb):
    d_bi[b, :] = np.random.normal(mud_b[b], sigd_b[b], nMCS)
    L_bi[b, :] = np.random.normal(muL_b[b], sigL_b[b], nMCS)

beta_bti0 = np.zeros((nb, Thoriz, nMCS))
for i in range(nMCS):
    beta_bti0[:, :, i] = reliability(beta0_b, d_bi[:, i], L_bi[:, i], Thoriz)  # reliability with time

plt.figure()
col = ['r', 'b', 'g', 'm']
for i in range(nMCS):
    for b in range(nb):
        plt.plot(range(1, Thoriz + 1), beta_bti0[b, :, i], f'-{col[b]}', linewidth=0.5)
for b in range(nb):
    minbeta_t0 = reliability([beta0_b[b]], [np.min(d_bi[b, :])], [np.min(L_bi[b, :])], Thoriz)[0,:]  # min reliability
    maxbeta_t0 = reliability([beta0_b[b]], [np.max(d_bi[b, :])], [np.max(L_bi[b, :])], Thoriz)[0,:]  # max reliability
    plt.plot(range(1, Thoriz + 1), minbeta_t0, f'-{col[b]}', linewidth=3)
    plt.plot(range(1, Thoriz + 1), maxbeta_t0, f'-{col[b]}', linewidth=3)
plt.grid()
plt.xlim([1, Thoriz])
plt.ylim([0, 1])
plt.xlabel('Time, t (years)')
plt.ylabel('Reliability index, β')
plt.title('Reliability Index Over Time')
plt.show()

# Deterministic Optimization
# Variables ->  e.g., freq=[10 2 0 0 ; 5, 5 0 0; 0 0 0 0];  

maxInter_b = np.ceil(Thoriz / np.min(L_bi, axis=1)).astype(int) + 1  # max number of interventions expected per bridge given by the lowest service life
nt = max(maxInter_b)  # to dimension the matrix, max interventions experienced by the worst case. 
nvars = nb * nt #number of design variables
integers = np.array([True] * nvars) #all variables are integers

#  Range of definition: [0,max number of years with no interventions==L_b]

lb = np.zeros((nvars, 1))
# upper bound
ub = lb.copy()
pos = 0
for b in range(nb):
    # This guarantees interventions freq. within the service life (if conducted)
    ub[pos:pos + maxInter_b[b]] = np.max(L_bi[b, :], axis=0) - 1  # the largest service life is of relevance here
    # this forces zeros after the max number of interventions expected
    if maxInter_b[b] < nb:
        ub[pos + maxInter_b[b]:pos + nb - 1] = 0
    pos += nt
ub = np.ceil(ub)

def nonlcon(y_b):
    #  From frequency to a schedule matrix [0/1]
    X_bt = Frequency2Schedule(y_b, nb, nt, Thoriz)
    
    # Constraint 1 >>> the last intervention is within Thoriz --> LastInterv<Thoriz --> LastInterv-Thoriz<=1   for all b
    lastInterv_b = np.max(cumulateFreq(y_b.reshape((nt, nb)).T), axis=1)
    c1 = lastInterv_b - Thoriz - 1

    # Constraint 2 >>> C_b*X_bt <= R   for all t
    c2 = np.sum(C_b[:, None] * X_bt, axis=0) - R

    # Constraint 3 >>> Prob(beta_bt>0)>=0.90  --> sum(sum(beta_bt))<=0 for all b
    ProbRel_b = ProbFulfilment(beta_bti0, X_bt, Thoriz)
    c3 = 0.90 - ProbRel_b

    c = np.concatenate([c1, c2, c3])
    return c

cons = sp.optimize.NonlinearConstraint(nonlcon, np.array([-np.inf] * (nb + Thoriz + nb)), np.array([0] * (nb + Thoriz + nb)))

def fun(y_b):
    #  minimize cost trying to go for the largest time between interventions (== max y_b)
    ninterv_b = CountNinterv(y_b, nt, nb)
    return np.dot(C_b, ninterv_b)

bounds = [[lb[i][0], ub[i][0]] for i in range(nvars)]

result = differential_evolution(fun, bounds, constraints=cons, integrality=integers, disp = True, maxiter = 2000)
print(result.x)

differential_evolution step 1: f(x)= 11.7
differential_evolution step 2: f(x)= 11.7
differential_evolution step 3: f(x)= 11.7
differential_evolution step 4: f(x)= 9.3
differential_evolution step 5: f(x)= 9.3
differential_evolution step 6: f(x)= 9.3
differential_evolution step 7: f(x)= 9.3
differential_evolution step 8: f(x)= 9.3
differential_evolution step 9: f(x)= 9.3
differential_evolution step 10: f(x)= 9.3
differential_evolution step 11: f(x)= 9.3
differential_evolution step 12: f(x)= 9.3
differential_evolution step 13: f(x)= 9.3
differential_evolution step 14: f(x)= 9.3
differential_evolution step 15: f(x)= 9.3
differential_evolution step 16: f(x)= 9.3
differential_evolution step 17: f(x)= 9.3
differential_evolution step 18: f(x)= 9.3
differential_evolution step 19: f(x)= 9.3
differential_evolution step 20: f(x)= 9.3
differential_evolution step 21: f(x)= 9.3
differential_evolution step 22: f(x)= 9.3
differential_evolution step 23: f(x)= 9.3
differential_evolution step 24: f(x)= 9.3
differential_evolution step 25: f(x)= 9.3
differential_evolution step 26: f(x)= 9.3
differential_evolution step 27: f(x)= 9.3
differential_evolution step 28: f(x)= 9.3
differential_evolution step 29: f(x)= 9.3
differential_evolution step 30: f(x)= 9.3
differential_evolution step 31: f(x)= 9.3
differential_evolution step 32: f(x)= 9.3
differential_evolution step 33: f(x)= 9.3
differential_evolution step 34: f(x)= 9.3
differential_evolution step 35: f(x)= 9.3
differential_evolution step 36: f(x)= 9.3
differential_evolution step 37: f(x)= 9.3
differential_evolution step 38: f(x)= 9.3
differential_evolution step 39: f(x)= 9.3
differential_evolution step 40: f(x)= 9.3
differential_evolution step 41: f(x)= 9.3
differential_evolution step 42: f(x)= 9.3
differential_evolution step 43: f(x)= 9.3
differential_evolution step 44: f(x)= 9.3
differential_evolution step 45: f(x)= 9.3
differential_evolution step 46: f(x)= 9.3
differential_evolution step 47: f(x)= 9.3
differential_evolution step 48: f(x)= 9.3
differential_evolution step 49: f(x)= 9.3
differential_evolution step 50: f(x)= 9.3
differential_evolution step 51: f(x)= 9.3
differential_evolution step 52: f(x)= 9.3
differential_evolution step 53: f(x)= 9.3
differential_evolution step 54: f(x)= 9.3
differential_evolution step 55: f(x)= 9.3
differential_evolution step 56: f(x)= 9.3
differential_evolution step 57: f(x)= 9.3
differential_evolution step 58: f(x)= 9.3
differential_evolution step 59: f(x)= 9.3
differential_evolution step 60: f(x)= 9.3
differential_evolution step 61: f(x)= 9.3
differential_evolution step 62: f(x)= 9.3
differential_evolution step 63: f(x)= 9.3
differential_evolution step 64: f(x)= 9.3
differential_evolution step 65: f(x)= 9.3
differential_evolution step 66: f(x)= 9.3
differential_evolution step 67: f(x)= 9.3
differential_evolution step 68: f(x)= 9.3
differential_evolution step 69: f(x)= 9.3
differential_evolution step 70: f(x)= 9.3
differential_evolution step 71: f(x)= 9.3
differential_evolution step 72: f(x)= 9.3
differential_evolution step 73: f(x)= 9.3
differential_evolution step 74: f(x)= 9.3
differential_evolution step 75: f(x)= 9.3
differential_evolution step 76: f(x)= 9.3
differential_evolution step 77: f(x)= 9.3
differential_evolution step 78: f(x)= 9.3
differential_evolution step 79: f(x)= 9.3
differential_evolution step 80: f(x)= 9.3
differential_evolution step 81: f(x)= 9.3
differential_evolution step 82: f(x)= 9.3
differential_evolution step 83: f(x)= 9.3
differential_evolution step 84: f(x)= 9.3
differential_evolution step 85: f(x)= 9.3
differential_evolution step 86: f(x)= 9.3
differential_evolution step 87: f(x)= 9.3
differential_evolution step 88: f(x)= 9.3
differential_evolution step 89: f(x)= 9.3
differential_evolution step 90: f(x)= 9.3
differential_evolution step 91: f(x)= 9.3
differential_evolution step 92: f(x)= 9.3
differential_evolution step 93: f(x)= 9.3
differential_evolution step 94: f(x)= 9.3
differential_evolution step 95: f(x)= 9.3
differential_evolution step 96: f(x)= 9.1
differential_evolution step 97: f(x)= 9.1
differential_evolution step 98: f(x)= 9.1
differential_evolution step 99: f(x)= 9.1
differential_evolution step 100: f(x)= 9.1
differential_evolution step 101: f(x)= 9.1
differential_evolution step 102: f(x)= 9.1
differential_evolution step 103: f(x)= 9.1
differential_evolution step 104: f(x)= 9.1
differential_evolution step 105: f(x)= 9.1
differential_evolution step 106: f(x)= 9.1
differential_evolution step 107: f(x)= 9.1
differential_evolution step 108: f(x)= 9.1
differential_evolution step 109: f(x)= 9.1
differential_evolution step 110: f(x)= 9.1
differential_evolution step 111: f(x)= 9.1
differential_evolution step 112: f(x)= 9.1
differential_evolution step 113: f(x)= 9.1
differential_evolution step 114: f(x)= 9.1
differential_evolution step 115: f(x)= 9.1
differential_evolution step 116: f(x)= 9.1
differential_evolution step 117: f(x)= 9.1
differential_evolution step 118: f(x)= 9.1
differential_evolution step 119: f(x)= 9.1
differential_evolution step 120: f(x)= 9.1
differential_evolution step 121: f(x)= 9.1
differential_evolution step 122: f(x)= 8.9
differential_evolution step 123: f(x)= 8.9
differential_evolution step 124: f(x)= 8.9
differential_evolution step 125: f(x)= 8.9
differential_evolution step 126: f(x)= 8.9
differential_evolution step 127: f(x)= 8.9
differential_evolution step 128: f(x)= 8.9
differential_evolution step 129: f(x)= 8.9
differential_evolution step 130: f(x)= 8.9
differential_evolution step 131: f(x)= 8.9
differential_evolution step 132: f(x)= 8.9
differential_evolution step 133: f(x)= 8.9
differential_evolution step 134: f(x)= 8.9
differential_evolution step 135: f(x)= 8.9
differential_evolution step 136: f(x)= 8.9
differential_evolution step 137: f(x)= 8.9
differential_evolution step 138: f(x)= 8.9
differential_evolution step 139: f(x)= 8.9
differential_evolution step 140: f(x)= 8.9
differential_evolution step 141: f(x)= 8.9
differential_evolution step 142: f(x)= 8.9
differential_evolution step 143: f(x)= 8.9
differential_evolution step 144: f(x)= 8.9
differential_evolution step 145: f(x)= 8.9
differential_evolution step 146: f(x)= 8.9
differential_evolution step 147: f(x)= 8.9
differential_evolution step 148: f(x)= 8.9
differential_evolution step 149: f(x)= 8.9
differential_evolution step 150: f(x)= 8.9
differential_evolution step 151: f(x)= 8.9
differential_evolution step 152: f(x)= 8.9
differential_evolution step 153: f(x)= 8.9
differential_evolution step 154: f(x)= 8.9
differential_evolution step 155: f(x)= 8.9
differential_evolution step 156: f(x)= 8.9
differential_evolution step 157: f(x)= 8.9
differential_evolution step 158: f(x)= 8.9
differential_evolution step 159: f(x)= 8.9
differential_evolution step 160: f(x)= 8.9
differential_evolution step 161: f(x)= 8.9
differential_evolution step 162: f(x)= 8.9
differential_evolution step 163: f(x)= 8.9
differential_evolution step 164: f(x)= 8.9
differential_evolution step 165: f(x)= 8.9
differential_evolution step 166: f(x)= 8.9
differential_evolution step 167: f(x)= 8.9
differential_evolution step 168: f(x)= 8.9
differential_evolution step 169: f(x)= 8.9
differential_evolution step 170: f(x)= 8.9
differential_evolution step 171: f(x)= 8.9
differential_evolution step 172: f(x)= 8.9
differential_evolution step 173: f(x)= 8.9
differential_evolution step 174: f(x)= 8.9
differential_evolution step 175: f(x)= 8.9
differential_evolution step 176: f(x)= 8.9
differential_evolution step 177: f(x)= 8.9
differential_evolution step 178: f(x)= 8.9
differential_evolution step 179: f(x)= 8.9
differential_evolution step 180: f(x)= 8.9
differential_evolution step 181: f(x)= 8.9
differential_evolution step 182: f(x)= 8.9
differential_evolution step 183: f(x)= 8.9
differential_evolution step 184: f(x)= 8.9
differential_evolution step 185: f(x)= 8.9
differential_evolution step 186: f(x)= 8.9
differential_evolution step 187: f(x)= 8.9
differential_evolution step 188: f(x)= 8.9
differential_evolution step 189: f(x)= 8.9
differential_evolution step 190: f(x)= 8.9
differential_evolution step 191: f(x)= 8.9
differential_evolution step 192: f(x)= 8.9
differential_evolution step 193: f(x)= 8.9
differential_evolution step 194: f(x)= 8.9
differential_evolution step 195: f(x)= 8.9
differential_evolution step 196: f(x)= 8.9
differential_evolution step 197: f(x)= 8.9
differential_evolution step 198: f(x)= 8.9
differential_evolution step 199: f(x)= 8.9
differential_evolution step 200: f(x)= 8.9
differential_evolution step 201: f(x)= 8.9
differential_evolution step 202: f(x)= 8.9
differential_evolution step 203: f(x)= 8.9
differential_evolution step 204: f(x)= 8.9
differential_evolution step 205: f(x)= 8.9
differential_evolution step 206: f(x)= 8.9
differential_evolution step 207: f(x)= 8.9
differential_evolution step 208: f(x)= 8.9
differential_evolution step 209: f(x)= 8.9
differential_evolution step 210: f(x)= 8.9
differential_evolution step 211: f(x)= 8.9
differential_evolution step 212: f(x)= 8.9
differential_evolution step 213: f(x)= 8.9
differential_evolution step 214: f(x)= 8.9
differential_evolution step 215: f(x)= 8.9
differential_evolution step 216: f(x)= 8.9
differential_evolution step 217: f(x)= 8.9
differential_evolution step 218: f(x)= 8.9
differential_evolution step 219: f(x)= 8.9
differential_evolution step 220: f(x)= 8.9
differential_evolution step 221: f(x)= 8.9
differential_evolution step 222: f(x)= 8.9
differential_evolution step 223: f(x)= 8.9
differential_evolution step 224: f(x)= 8.9
differential_evolution step 225: f(x)= 8.9
differential_evolution step 226: f(x)= 8.9
differential_evolution step 227: f(x)= 8.9
differential_evolution step 228: f(x)= 8.9
differential_evolution step 229: f(x)= 8.9
differential_evolution step 230: f(x)= 8.9
differential_evolution step 231: f(x)= 8.9
differential_evolution step 232: f(x)= 8.9
differential_evolution step 233: f(x)= 8.9
differential_evolution step 234: f(x)= 8.9
differential_evolution step 235: f(x)= 8.9
differential_evolution step 236: f(x)= 8.9
differential_evolution step 237: f(x)= 8.9
differential_evolution step 238: f(x)= 8.9
differential_evolution step 239: f(x)= 8.9
differential_evolution step 240: f(x)= 8.9
differential_evolution step 241: f(x)= 8.9
differential_evolution step 242: f(x)= 8.9
differential_evolution step 243: f(x)= 8.9
differential_evolution step 244: f(x)= 8.9
differential_evolution step 245: f(x)= 8.9
differential_evolution step 246: f(x)= 8.9
differential_evolution step 247: f(x)= 8.9
differential_evolution step 248: f(x)= 8.9
differential_evolution step 249: f(x)= 8.9
differential_evolution step 250: f(x)= 8.9
differential_evolution step 251: f(x)= 8.9
differential_evolution step 252: f(x)= 8.9
differential_evolution step 253: f(x)= 8.9
differential_evolution step 254: f(x)= 8.9
differential_evolution step 255: f(x)= 8.9
differential_evolution step 256: f(x)= 8.9
differential_evolution step 257: f(x)= 8.9
differential_evolution step 258: f(x)= 8.9
differential_evolution step 259: f(x)= 8.9
differential_evolution step 260: f(x)= 8.9
differential_evolution step 261: f(x)= 8.9
differential_evolution step 262: f(x)= 8.9
differential_evolution step 263: f(x)= 8.9
differential_evolution step 264: f(x)= 8.9
differential_evolution step 265: f(x)= 8.9
differential_evolution step 266: f(x)= 8.9
differential_evolution step 267: f(x)= 8.9
differential_evolution step 268: f(x)= 8.9
differential_evolution step 269: f(x)= 8.9
differential_evolution step 270: f(x)= 8.9
differential_evolution step 271: f(x)= 8.9
differential_evolution step 272: f(x)= 8.9
differential_evolution step 273: f(x)= 8.9
differential_evolution step 274: f(x)= 8.9
differential_evolution step 275: f(x)= 8.9
differential_evolution step 276: f(x)= 8.9
differential_evolution step 277: f(x)= 8.9
differential_evolution step 278: f(x)= 8.9
differential_evolution step 279: f(x)= 8.9
differential_evolution step 280: f(x)= 8.4
differential_evolution step 281: f(x)= 8.4
differential_evolution step 282: f(x)= 8.4
differential_evolution step 283: f(x)= 8.4
differential_evolution step 284: f(x)= 8.4
differential_evolution step 285: f(x)= 8.4
differential_evolution step 286: f(x)= 8.4
differential_evolution step 287: f(x)= 8.4
differential_evolution step 288: f(x)= 8.4
differential_evolution step 289: f(x)= 8.4
differential_evolution step 290: f(x)= 8.4
differential_evolution step 291: f(x)= 8.4
differential_evolution step 292: f(x)= 8.4
differential_evolution step 293: f(x)= 8.4
differential_evolution step 294: f(x)= 8.4
differential_evolution step 295: f(x)= 8.4
differential_evolution step 296: f(x)= 8.4
differential_evolution step 297: f(x)= 8.4
differential_evolution step 298: f(x)= 8.4
differential_evolution step 299: f(x)= 8.4
differential_evolution step 300: f(x)= 8.4
differential_evolution step 301: f(x)= 8.4
differential_evolution step 302: f(x)= 8.4
differential_evolution step 303: f(x)= 8.4
differential_evolution step 304: f(x)= 8.4
differential_evolution step 305: f(x)= 8.4
differential_evolution step 306: f(x)= 8.4
differential_evolution step 307: f(x)= 8.4
differential_evolution step 308: f(x)= 8.4
differential_evolution step 309: f(x)= 8.4
differential_evolution step 310: f(x)= 8.4
differential_evolution step 311: f(x)= 8.4
differential_evolution step 312: f(x)= 8.4
differential_evolution step 313: f(x)= 8.4
differential_evolution step 314: f(x)= 8.4
differential_evolution step 315: f(x)= 8.4
differential_evolution step 316: f(x)= 8.4
differential_evolution step 317: f(x)= 8.4
differential_evolution step 318: f(x)= 8.4
differential_evolution step 319: f(x)= 8.4
differential_evolution step 320: f(x)= 8.4
differential_evolution step 321: f(x)= 8.4
differential_evolution step 322: f(x)= 8.4
differential_evolution step 323: f(x)= 8.4
differential_evolution step 324: f(x)= 8.4
differential_evolution step 325: f(x)= 8.4
differential_evolution step 326: f(x)= 8.4
differential_evolution step 327: f(x)= 8.4
differential_evolution step 328: f(x)= 8.200000000000001
differential_evolution step 329: f(x)= 8.200000000000001
differential_evolution step 330: f(x)= 8.200000000000001
differential_evolution step 331: f(x)= 8.200000000000001
differential_evolution step 332: f(x)= 8.200000000000001
differential_evolution step 333: f(x)= 8.200000000000001
differential_evolution step 334: f(x)= 8.200000000000001
differential_evolution step 335: f(x)= 8.200000000000001
differential_evolution step 336: f(x)= 8.200000000000001
differential_evolution step 337: f(x)= 8.200000000000001
differential_evolution step 338: f(x)= 8.200000000000001
differential_evolution step 339: f(x)= 8.200000000000001
differential_evolution step 340: f(x)= 8.200000000000001
differential_evolution step 341: f(x)= 8.200000000000001
differential_evolution step 342: f(x)= 8.200000000000001
differential_evolution step 343: f(x)= 8.200000000000001
differential_evolution step 344: f(x)= 8.200000000000001
differential_evolution step 345: f(x)= 8.200000000000001
differential_evolution step 346: f(x)= 8.200000000000001
differential_evolution step 347: f(x)= 8.200000000000001
differential_evolution step 348: f(x)= 8.200000000000001
differential_evolution step 349: f(x)= 8.200000000000001
differential_evolution step 350: f(x)= 8.200000000000001
differential_evolution step 351: f(x)= 8.200000000000001
differential_evolution step 352: f(x)= 8.200000000000001
differential_evolution step 353: f(x)= 8.200000000000001
differential_evolution step 354: f(x)= 8.200000000000001
differential_evolution step 355: f(x)= 8.200000000000001
differential_evolution step 356: f(x)= 8.200000000000001
differential_evolution step 357: f(x)= 8.200000000000001
differential_evolution step 358: f(x)= 8.200000000000001
differential_evolution step 359: f(x)= 8.200000000000001
differential_evolution step 360: f(x)= 8.200000000000001
differential_evolution step 361: f(x)= 8.200000000000001
differential_evolution step 362: f(x)= 8.200000000000001
differential_evolution step 363: f(x)= 8.200000000000001
differential_evolution step 364: f(x)= 8.200000000000001
differential_evolution step 365: f(x)= 8.200000000000001
differential_evolution step 366: f(x)= 8.200000000000001
differential_evolution step 367: f(x)= 8.200000000000001
differential_evolution step 368: f(x)= 8.200000000000001
differential_evolution step 369: f(x)= 8.200000000000001
differential_evolution step 370: f(x)= 8.200000000000001
differential_evolution step 371: f(x)= 8.200000000000001
differential_evolution step 372: f(x)= 8.200000000000001
differential_evolution step 373: f(x)= 8.200000000000001
differential_evolution step 374: f(x)= 8.200000000000001
differential_evolution step 375: f(x)= 8.200000000000001
differential_evolution step 376: f(x)= 8.200000000000001
differential_evolution step 377: f(x)= 8.200000000000001
differential_evolution step 378: f(x)= 8.200000000000001
differential_evolution step 379: f(x)= 8.200000000000001
differential_evolution step 380: f(x)= 8.200000000000001
differential_evolution step 381: f(x)= 8.200000000000001
differential_evolution step 382: f(x)= 8.200000000000001
differential_evolution step 383: f(x)= 8.200000000000001
differential_evolution step 384: f(x)= 8.200000000000001
differential_evolution step 385: f(x)= 8.200000000000001
differential_evolution step 386: f(x)= 8.200000000000001
differential_evolution step 387: f(x)= 8.200000000000001
differential_evolution step 388: f(x)= 8.200000000000001
differential_evolution step 389: f(x)= 8.200000000000001
differential_evolution step 390: f(x)= 8.200000000000001
differential_evolution step 391: f(x)= 8.200000000000001
differential_evolution step 392: f(x)= 8.200000000000001
differential_evolution step 393: f(x)= 8.200000000000001
differential_evolution step 394: f(x)= 8.200000000000001
differential_evolution step 395: f(x)= 8.200000000000001
differential_evolution step 396: f(x)= 8.200000000000001
differential_evolution step 397: f(x)= 8.200000000000001
differential_evolution step 398: f(x)= 8.200000000000001
differential_evolution step 399: f(x)= 8.200000000000001
differential_evolution step 400: f(x)= 8.200000000000001
differential_evolution step 401: f(x)= 8.200000000000001
differential_evolution step 402: f(x)= 8.200000000000001
differential_evolution step 403: f(x)= 8.200000000000001
differential_evolution step 404: f(x)= 8.200000000000001
differential_evolution step 405: f(x)= 8.200000000000001
differential_evolution step 406: f(x)= 8.200000000000001
differential_evolution step 407: f(x)= 8.200000000000001
differential_evolution step 408: f(x)= 8.200000000000001
differential_evolution step 409: f(x)= 8.200000000000001
differential_evolution step 410: f(x)= 8.200000000000001
differential_evolution step 411: f(x)= 8.200000000000001
differential_evolution step 412: f(x)= 8.200000000000001
differential_evolution step 413: f(x)= 8.200000000000001
differential_evolution step 414: f(x)= 8.200000000000001
differential_evolution step 415: f(x)= 8.200000000000001
differential_evolution step 416: f(x)= 8.200000000000001
differential_evolution step 417: f(x)= 8.200000000000001
differential_evolution step 418: f(x)= 8.200000000000001
differential_evolution step 419: f(x)= 8.200000000000001
differential_evolution step 420: f(x)= 8.200000000000001
differential_evolution step 421: f(x)= 8.200000000000001
differential_evolution step 422: f(x)= 8.200000000000001
differential_evolution step 423: f(x)= 8.200000000000001
differential_evolution step 424: f(x)= 8.200000000000001
differential_evolution step 425: f(x)= 8.200000000000001
differential_evolution step 426: f(x)= 8.200000000000001
differential_evolution step 427: f(x)= 8.200000000000001
differential_evolution step 428: f(x)= 8.200000000000001
differential_evolution step 429: f(x)= 8.200000000000001
differential_evolution step 430: f(x)= 8.200000000000001
differential_evolution step 431: f(x)= 8.200000000000001
differential_evolution step 432: f(x)= 8.200000000000001
differential_evolution step 433: f(x)= 8.200000000000001
differential_evolution step 434: f(x)= 8.200000000000001
differential_evolution step 435: f(x)= 8.200000000000001
differential_evolution step 436: f(x)= 8.200000000000001
differential_evolution step 437: f(x)= 8.200000000000001
differential_evolution step 438: f(x)= 8.200000000000001
differential_evolution step 439: f(x)= 8.200000000000001
differential_evolution step 440: f(x)= 8.200000000000001
differential_evolution step 441: f(x)= 8.200000000000001
differential_evolution step 442: f(x)= 8.200000000000001
differential_evolution step 443: f(x)= 8.200000000000001
differential_evolution step 444: f(x)= 8.200000000000001
differential_evolution step 445: f(x)= 8.200000000000001
differential_evolution step 446: f(x)= 8.200000000000001
differential_evolution step 447: f(x)= 8.200000000000001
differential_evolution step 448: f(x)= 8.200000000000001
differential_evolution step 449: f(x)= 8.200000000000001
differential_evolution step 450: f(x)= 8.200000000000001
differential_evolution step 451: f(x)= 8.200000000000001
differential_evolution step 452: f(x)= 8.200000000000001
differential_evolution step 453: f(x)= 8.200000000000001
differential_evolution step 454: f(x)= 8.200000000000001
differential_evolution step 455: f(x)= 8.200000000000001
differential_evolution step 456: f(x)= 8.200000000000001
differential_evolution step 457: f(x)= 8.200000000000001
differential_evolution step 458: f(x)= 8.200000000000001
differential_evolution step 459: f(x)= 8.200000000000001
differential_evolution step 460: f(x)= 8.200000000000001
differential_evolution step 461: f(x)= 8.200000000000001
differential_evolution step 462: f(x)= 8.200000000000001
differential_evolution step 463: f(x)= 8.200000000000001
differential_evolution step 464: f(x)= 8.200000000000001
differential_evolution step 465: f(x)= 8.200000000000001
differential_evolution step 466: f(x)= 8.200000000000001
differential_evolution step 467: f(x)= 8.200000000000001
differential_evolution step 468: f(x)= 8.200000000000001
differential_evolution step 469: f(x)= 8.200000000000001
differential_evolution step 470: f(x)= 8.200000000000001
differential_evolution step 471: f(x)= 8.200000000000001
differential_evolution step 472: f(x)= 8.200000000000001
differential_evolution step 473: f(x)= 8.200000000000001
differential_evolution step 474: f(x)= 8.200000000000001
differential_evolution step 475: f(x)= 8.200000000000001
differential_evolution step 476: f(x)= 8.200000000000001
differential_evolution step 477: f(x)= 8.200000000000001
differential_evolution step 478: f(x)= 8.200000000000001
differential_evolution step 479: f(x)= 8.200000000000001
differential_evolution step 480: f(x)= 8.200000000000001
differential_evolution step 481: f(x)= 8.200000000000001
differential_evolution step 482: f(x)= 8.200000000000001
differential_evolution step 483: f(x)= 8.200000000000001
differential_evolution step 484: f(x)= 7.500000000000001
differential_evolution step 485: f(x)= 7.500000000000001
differential_evolution step 486: f(x)= 7.500000000000001
differential_evolution step 487: f(x)= 7.500000000000001
differential_evolution step 488: f(x)= 7.500000000000001
differential_evolution step 489: f(x)= 7.500000000000001
differential_evolution step 490: f(x)= 7.500000000000001
differential_evolution step 491: f(x)= 7.500000000000001
differential_evolution step 492: f(x)= 7.500000000000001
differential_evolution step 493: f(x)= 7.500000000000001
differential_evolution step 494: f(x)= 7.500000000000001
differential_evolution step 495: f(x)= 7.500000000000001
differential_evolution step 496: f(x)= 7.500000000000001
differential_evolution step 497: f(x)= 7.500000000000001
differential_evolution step 498: f(x)= 7.500000000000001
differential_evolution step 499: f(x)= 7.500000000000001
differential_evolution step 500: f(x)= 7.500000000000001
differential_evolution step 501: f(x)= 7.500000000000001
differential_evolution step 502: f(x)= 7.500000000000001
differential_evolution step 503: f(x)= 7.500000000000001
differential_evolution step 504: f(x)= 7.500000000000001
differential_evolution step 505: f(x)= 7.500000000000001
differential_evolution step 506: f(x)= 7.500000000000001
differential_evolution step 507: f(x)= 7.500000000000001
differential_evolution step 508: f(x)= 7.500000000000001
differential_evolution step 509: f(x)= 7.500000000000001
differential_evolution step 510: f(x)= 7.500000000000001
differential_evolution step 511: f(x)= 7.500000000000001
differential_evolution step 512: f(x)= 7.500000000000001
differential_evolution step 513: f(x)= 7.500000000000001
differential_evolution step 514: f(x)= 7.500000000000001
differential_evolution step 515: f(x)= 7.500000000000001
differential_evolution step 516: f(x)= 7.500000000000001
differential_evolution step 517: f(x)= 7.500000000000001
differential_evolution step 518: f(x)= 7.500000000000001
differential_evolution step 519: f(x)= 7.500000000000001
differential_evolution step 520: f(x)= 7.500000000000001
differential_evolution step 521: f(x)= 7.500000000000001
differential_evolution step 522: f(x)= 7.500000000000001
differential_evolution step 523: f(x)= 7.500000000000001
differential_evolution step 524: f(x)= 7.500000000000001
differential_evolution step 525: f(x)= 7.500000000000001
differential_evolution step 526: f(x)= 7.500000000000001
differential_evolution step 527: f(x)= 7.500000000000001
differential_evolution step 528: f(x)= 7.500000000000001
differential_evolution step 529: f(x)= 7.500000000000001
differential_evolution step 530: f(x)= 7.500000000000001
differential_evolution step 531: f(x)= 7.500000000000001
differential_evolution step 532: f(x)= 7.500000000000001
differential_evolution step 533: f(x)= 7.500000000000001
differential_evolution step 534: f(x)= 7.500000000000001
differential_evolution step 535: f(x)= 7.500000000000001
differential_evolution step 536: f(x)= 7.500000000000001
differential_evolution step 537: f(x)= 7.500000000000001
differential_evolution step 538: f(x)= 7.500000000000001
differential_evolution step 539: f(x)= 7.500000000000001
differential_evolution step 540: f(x)= 7.500000000000001
differential_evolution step 541: f(x)= 7.500000000000001
differential_evolution step 542: f(x)= 7.500000000000001
differential_evolution step 543: f(x)= 7.500000000000001
differential_evolution step 544: f(x)= 7.500000000000001
differential_evolution step 545: f(x)= 7.500000000000001
differential_evolution step 546: f(x)= 7.500000000000001
differential_evolution step 547: f(x)= 7.500000000000001
differential_evolution step 548: f(x)= 7.500000000000001
differential_evolution step 549: f(x)= 7.500000000000001
differential_evolution step 550: f(x)= 7.500000000000001
differential_evolution step 551: f(x)= 7.500000000000001
differential_evolution step 552: f(x)= 7.500000000000001
differential_evolution step 553: f(x)= 7.500000000000001
differential_evolution step 554: f(x)= 7.500000000000001
differential_evolution step 555: f(x)= 7.500000000000001
differential_evolution step 556: f(x)= 7.500000000000001
differential_evolution step 557: f(x)= 7.500000000000001
differential_evolution step 558: f(x)= 7.500000000000001
differential_evolution step 559: f(x)= 7.500000000000001
differential_evolution step 560: f(x)= 7.500000000000001
differential_evolution step 561: f(x)= 7.500000000000001
differential_evolution step 562: f(x)= 7.500000000000001
differential_evolution step 563: f(x)= 7.500000000000001
differential_evolution step 564: f(x)= 7.500000000000001
differential_evolution step 565: f(x)= 7.500000000000001
differential_evolution step 566: f(x)= 7.500000000000001
differential_evolution step 567: f(x)= 7.500000000000001
differential_evolution step 568: f(x)= 7.500000000000001
differential_evolution step 569: f(x)= 7.500000000000001
differential_evolution step 570: f(x)= 7.500000000000001
differential_evolution step 571: f(x)= 7.500000000000001
differential_evolution step 572: f(x)= 7.500000000000001
differential_evolution step 573: f(x)= 7.500000000000001
differential_evolution step 574: f(x)= 7.500000000000001
differential_evolution step 575: f(x)= 7.500000000000001
differential_evolution step 576: f(x)= 7.500000000000001
differential_evolution step 577: f(x)= 7.500000000000001
differential_evolution step 578: f(x)= 7.500000000000001
differential_evolution step 579: f(x)= 7.500000000000001
differential_evolution step 580: f(x)= 7.500000000000001
differential_evolution step 581: f(x)= 7.500000000000001
differential_evolution step 582: f(x)= 7.500000000000001
differential_evolution step 583: f(x)= 7.500000000000001
differential_evolution step 584: f(x)= 7.500000000000001
differential_evolution step 585: f(x)= 7.500000000000001
differential_evolution step 586: f(x)= 7.500000000000001
differential_evolution step 587: f(x)= 7.500000000000001
differential_evolution step 588: f(x)= 7.500000000000001
differential_evolution step 589: f(x)= 7.500000000000001
differential_evolution step 590: f(x)= 7.500000000000001
differential_evolution step 591: f(x)= 6.9
differential_evolution step 592: f(x)= 6.9
differential_evolution step 593: f(x)= 6.9
differential_evolution step 594: f(x)= 6.9
differential_evolution step 595: f(x)= 6.9
differential_evolution step 596: f(x)= 6.9
differential_evolution step 597: f(x)= 6.9
differential_evolution step 598: f(x)= 6.9
differential_evolution step 599: f(x)= 6.9
differential_evolution step 600: f(x)= 6.9
differential_evolution step 601: f(x)= 6.9
differential_evolution step 602: f(x)= 6.9
differential_evolution step 603: f(x)= 6.9
differential_evolution step 604: f(x)= 6.9
differential_evolution step 605: f(x)= 6.9
differential_evolution step 606: f(x)= 6.9
differential_evolution step 607: f(x)= 6.9
differential_evolution step 608: f(x)= 6.9
differential_evolution step 609: f(x)= 6.9
differential_evolution step 610: f(x)= 6.9
differential_evolution step 611: f(x)= 6.9
differential_evolution step 612: f(x)= 6.9
differential_evolution step 613: f(x)= 6.9
differential_evolution step 614: f(x)= 6.9
differential_evolution step 615: f(x)= 6.9
differential_evolution step 616: f(x)= 6.9
differential_evolution step 617: f(x)= 6.9
differential_evolution step 618: f(x)= 6.9
differential_evolution step 619: f(x)= 6.9
differential_evolution step 620: f(x)= 6.9
differential_evolution step 621: f(x)= 6.9
differential_evolution step 622: f(x)= 6.9
differential_evolution step 623: f(x)= 6.9
differential_evolution step 624: f(x)= 6.9
differential_evolution step 625: f(x)= 6.9
differential_evolution step 626: f(x)= 6.9
differential_evolution step 627: f(x)= 6.9
differential_evolution step 628: f(x)= 6.9
differential_evolution step 629: f(x)= 6.9
differential_evolution step 630: f(x)= 6.9
differential_evolution step 631: f(x)= 6.9
differential_evolution step 632: f(x)= 6.9
differential_evolution step 633: f(x)= 6.9
differential_evolution step 634: f(x)= 6.9
differential_evolution step 635: f(x)= 6.9
differential_evolution step 636: f(x)= 6.9
differential_evolution step 637: f(x)= 6.9
differential_evolution step 638: f(x)= 6.9
differential_evolution step 639: f(x)= 6.9
differential_evolution step 640: f(x)= 6.9
differential_evolution step 641: f(x)= 6.9
differential_evolution step 642: f(x)= 6.9
differential_evolution step 643: f(x)= 6.9
differential_evolution step 644: f(x)= 6.9
differential_evolution step 645: f(x)= 6.9
differential_evolution step 646: f(x)= 6.9
differential_evolution step 647: f(x)= 6.9
differential_evolution step 648: f(x)= 6.9
differential_evolution step 649: f(x)= 6.9
differential_evolution step 650: f(x)= 6.9
differential_evolution step 651: f(x)= 6.9
differential_evolution step 652: f(x)= 6.9
differential_evolution step 653: f(x)= 6.9
differential_evolution step 654: f(x)= 6.9
differential_evolution step 655: f(x)= 6.9
differential_evolution step 656: f(x)= 6.9
differential_evolution step 657: f(x)= 6.9
differential_evolution step 658: f(x)= 6.9
differential_evolution step 659: f(x)= 6.9
differential_evolution step 660: f(x)= 6.9
differential_evolution step 661: f(x)= 6.9
differential_evolution step 662: f(x)= 6.9
differential_evolution step 663: f(x)= 6.9
differential_evolution step 664: f(x)= 6.9
differential_evolution step 665: f(x)= 6.9
differential_evolution step 666: f(x)= 6.9
differential_evolution step 667: f(x)= 6.9
differential_evolution step 668: f(x)= 6.9
differential_evolution step 669: f(x)= 6.9
differential_evolution step 670: f(x)= 6.9
differential_evolution step 671: f(x)= 6.9
differential_evolution step 672: f(x)= 6.9
differential_evolution step 673: f(x)= 6.9
differential_evolution step 674: f(x)= 6.9
differential_evolution step 675: f(x)= 6.9
differential_evolution step 676: f(x)= 6.9
differential_evolution step 677: f(x)= 6.9
differential_evolution step 678: f(x)= 6.1
differential_evolution step 679: f(x)= 6.1
differential_evolution step 680: f(x)= 6.1
differential_evolution step 681: f(x)= 6.1
differential_evolution step 682: f(x)= 6.1
differential_evolution step 683: f(x)= 6.1
differential_evolution step 684: f(x)= 6.1
differential_evolution step 685: f(x)= 6.1
differential_evolution step 686: f(x)= 6.1
differential_evolution step 687: f(x)= 6.1
differential_evolution step 688: f(x)= 6.1
differential_evolution step 689: f(x)= 6.1
differential_evolution step 690: f(x)= 6.1
differential_evolution step 691: f(x)= 6.1
differential_evolution step 692: f(x)= 6.1
differential_evolution step 693: f(x)= 6.1
differential_evolution step 694: f(x)= 6.1
differential_evolution step 695: f(x)= 6.1
differential_evolution step 696: f(x)= 6.1
differential_evolution step 697: f(x)= 6.1
differential_evolution step 698: f(x)= 6.1
differential_evolution step 699: f(x)= 6.1
differential_evolution step 700: f(x)= 6.1
differential_evolution step 701: f(x)= 6.1
differential_evolution step 702: f(x)= 6.1
differential_evolution step 703: f(x)= 6.1
differential_evolution step 704: f(x)= 6.1
differential_evolution step 705: f(x)= 6.1
differential_evolution step 706: f(x)= 6.1
differential_evolution step 707: f(x)= 6.1
differential_evolution step 708: f(x)= 6.1
differential_evolution step 709: f(x)= 6.1
differential_evolution step 710: f(x)= 6.1
differential_evolution step 711: f(x)= 6.1
differential_evolution step 712: f(x)= 6.1
differential_evolution step 713: f(x)= 6.1
differential_evolution step 714: f(x)= 6.1
differential_evolution step 715: f(x)= 6.1
differential_evolution step 716: f(x)= 6.1
differential_evolution step 717: f(x)= 6.1
differential_evolution step 718: f(x)= 6.1
differential_evolution step 719: f(x)= 6.1
differential_evolution step 720: f(x)= 6.1
differential_evolution step 721: f(x)= 6.1
differential_evolution step 722: f(x)= 6.1
differential_evolution step 723: f(x)= 6.1
differential_evolution step 724: f(x)= 6.1
differential_evolution step 725: f(x)= 6.1
differential_evolution step 726: f(x)= 6.1
differential_evolution step 727: f(x)= 6.1
differential_evolution step 728: f(x)= 6.1
differential_evolution step 729: f(x)= 6.1
differential_evolution step 730: f(x)= 6.1
differential_evolution step 731: f(x)= 6.1
differential_evolution step 732: f(x)= 6.1
differential_evolution step 733: f(x)= 6.1
differential_evolution step 734: f(x)= 6.1
differential_evolution step 735: f(x)= 6.1
differential_evolution step 736: f(x)= 6.1
differential_evolution step 737: f(x)= 6.1
differential_evolution step 738: f(x)= 6.1
differential_evolution step 739: f(x)= 6.1
differential_evolution step 740: f(x)= 6.1
[ 7.  0.  2.  0.  8.  8.  0.  0.  6.  9. 12. 10.  0.  0.  0.  0.]

Yopt = result.x
OptInterv = cumulateFreq(Yopt.reshape((nt, nb)).T)
X_bt = Frequency2Schedule(Yopt, nb, nt, Thoriz)

print(f'Min Cost = {result.fun:.2f}')
ninterv = np.sum(CountNinterv(Yopt, nt, nb))
print(f'Number of interventions = {ninterv}')
c = nonlcon(Yopt)
BudgetIssues = np.sum(c[-Thoriz - 1:-1] > 0)
print(f'Number of time intervals with exceeding budget = {BudgetIssues}')
for b in range(1, nb + 1):
    probfulfilment = 0.90 - c[-nb + b - 1]
    print(f'Prob(reliability>0) for bridge {b} = {probfulfilment:.4f}')

plt.figure()
h = []
leg = []
for i in range(nMCS):
    beta_bt = MaintenanceApplication(beta_bti0[:, :, i], X_bt, Thoriz)  # optimal maintenance
    for b in range(nb):
        h.append(plt.plot(range(1, Thoriz + 1), beta_bt[b, :], f'-{col[b]}', linewidth=1))
        if i == 0:
            leg.append(f'Bridge {b + 1}')
plt.legend(leg, loc='best')
plt.grid()
plt.xlim([1, Thoriz])
plt.ylim([0, 1])
plt.xlabel('Time, t (years)')
plt.ylabel('Reliability index, β')
plt.title('Reliability Index Over Time')
plt.show()

Min Cost = 6.10
Number of interventions = 8
Number of time intervals with exceeding budget = 0
Prob(reliability>0) for bridge 1 = 1.0000
Prob(reliability>0) for bridge 2 = 0.9000
Prob(reliability>0) for bridge 3 = 1.0000
Prob(reliability>0) for bridge 4 = 1.0000