'''

 

In a workshop there are two production lines in parallel. Each production line 
has a critical machine with constant failure rate λ = 0.01 failures per hour. 
There is one (common) spare machine that can replace a failed machine.
We assume that switching time can be ignored. The repair rate of the machines 
is assumed constant and equal to µ = 0.2 per hour. If a production line
is down the loss is assumed to be c_U = 10 000 NOKs per hour. 
Only one repair man is available.

States:

3: Two machines are functioning
2: One machine is failed, the spare machine takes over, i.e., two machines are running
1: Two machines are failed, only one is running
0: All machines are failed


A = | −µ   µ    0    0 |
    |  λ −λ−µ   µ    0 |
    |  0  2λ  −2λ−µ  µ |
    |  0   0   2λ  −2λ |


First we find the steady state solution, then the time dependent after maxTime = 2 

Time unit is hours

'''
import numpy as np

mu = 0.2
lmbda = 0.01
c_U = 10000
r = 3   # States are 0,1,2,...,r=3
N = r+1 # Number of states
A = np.zeros([N, N])  # cerates a matrix of dimension (0 to N-1) * (0 to N-1), i.e., (0 to r) * (0 to r)
#
# The failures are specified directly into the A-matrix
#
A[3,2] = 2*lmbda    # Two machines can fail
A[2,1] = 2*lmbda    # Two machines can fail
A[1,0] = lmbda      # One machine can fail

for i in range(0,N-1):  # mu values above diagonal   
    A[i,i+1] = mu
    
def fixA(A):  
    N = len(A)     
    for i in range(0,N):    # Fill in diagonal elements
        s = 0
        for j in range(0,N):
            if i!=j: 
                s += A[i,j]
        A[i,i] = -s

def P_asymptotic(A):  # Find asymptotic solution
    N = len(A) 
    A1=np.copy(A)
    for i in range(0,N):
        A1[i,0] = 1
    B=np.zeros([N])    
    B[0] = 1
    P=np.linalg.solve(np.transpose(A1),B)
    return P

def P_timeDependent(A,t,dt,initState): # Time dependent solution at time t, starts in InitState, integration step = dt 
    N = len(A)
    P=np.zeros([N])
    P[initState] = 1 
    IM = np.eye(N) + np.dot(A,dt) # Create integration matrix, i.e., IM = (I + A*dt)
    for t in np.arange(0,t,dt): # note, the arange() is a range() function where increments may be decimals
        P= np.dot(P, IM )
    return P
    
def P_average(A,t,dt,initState):
    cnt=0
    N = len(A)
    P=np.zeros([N])
    P[initState] = 1 
    pAvg=P*0.5
    IM = np.eye(N) + np.dot(A,dt) # Create integration matrix, i.e., IM = (I + A*dt)
    for t in np.arange(0,t,dt): # note, the arange() is a range() function where increments may be decimals
        P= np.dot(P, IM )
        pAvg += P
        cnt += 1
    pAvg -= P*0.5  # Trapezoidal formula
    return pAvg/cnt    
    
def printVec(lab,P):
    print("\n"+lab)
    for i , value in enumerate(P): 
        print(f'P_{i}:','{:.3E}'.format(value))
        
leadTxt = "\nOne repair man *****************"
uCost=[0,0]
for i in range(2) :   # Repeat twice
    fixA(A)
    P=P_asymptotic(A)
    Pt=P_timeDependent(A,8,0.01,3)
    Pa=P_average(A,8,0.01,3)
    print(leadTxt)
    printVec("Steady state:",P)
    print("Average U-cost",P[0]*2*c_U+P[1]*c_U)
    print("Average U-cost (per day)",(P[0]*2*c_U+P[1]*c_U)*8)
    uCost[i] = P[0]*2*c_U+P[1]*c_U
    printVec ("P(t=8 hours):",Pt)
    printVec ("Average P(t):",Pa)
    print("Average U-cost 7-15",Pa[0]*2*c_U+Pa[1]*c_U)
    print("Average U-cost (per day)",(Pa[0]*2*c_U+Pa[1]*c_U)*8)
    leadTxt = "\nTwo repair men *****************"
    '''
    With two repair men, we have the following matrix
    A = | −2µ  2µ     0     0 |
        |  λ −λ−2µ   2µ     0 |
        |  0  2λ    −2λ−2µ  µ |
        |  0   0     2λ   −2λ |
    '''
    A[0,1] = 2*mu
    A[1,2] = 2*mu
    
print("\nSummary results for decision:")    
print("Willing to pay per hour, 24/7:",uCost[0]-uCost[1])