Linear Algebraic Queueing Theory Approach to Reliability

Mike Ulrey

Boeing

Overview

Main themes

-- Main ideas taken from and inspired by “Queueing Theory: A Linear Algebraic Approach” by Lester Lipsky, Springer, 2nd Ed, 2009.

-- Method to compute system survival (reliability) functions based on
	-- Stochastic behavior(s) of the individual components (of the system)
	-- Topology of the system
-- Method produces a representation that is
	-- Symbolic
	-- Closed-form
	-- Fast and easy to compute
	-- Applicable to a wide range of systems
-- “Duality” between reliability and queueing theory
-- “Magic” formula involving the matrix exponential
-- Symmetry appears in many different guises

Assumptions

-- Components fail independently of each other
-- Component waiting times to failure follow a (negative) exponential distribution (constant failure rate)

Disclaimer

-- Nothing is novel here except perhaps highlighting a method for unifying several known results using a single technique

Component reliability

Single component -- two states

Exponential distribution

Unreliability (time to failure) and reliability (survival) functions

Grid[{Map[Text[#]&,{"Exponential PDF","Exponential CDF (Unreliability)","Survival Function (Reliability)"}],{PDF[ExponentialDistribution[λ],t],unReliability=CDF[ExponentialDistribution[λ],t],reliability=SurvivalFunction[ExponentialDistribution[λ],t]}},Frame->All]
Show[Plot[reliability/.{λ->0.1},{t,0,25},PlotStyle->{Green}],Plot[unReliability/.{λ->0.1},{t,0,25},PlotStyle->{Red}],AxesLabel->{"Time (hr)","Reliability"},PlotRange->{Automatic,All},Epilog->{Inset[Text["Unreliability (λ=0.1)"],{19,0.8}],Inset[Text["Reliability (λ=0.1)"],{19,0.2}]}]

Exponential distribution characterized by a constant failure rate (hazard function)

HazardFunction[ExponentialDistribution[λ],t]

“Mean time before failure (MTBF)” = Inverse of failure rate. 
E.g. if λ=0.1 failures per hour, MTBF=10 hrs.

Mean[ExponentialDistribution[λ]]

Memoryless property: If component survives up to time s, probability that it survives until time t+s is the same as the probability that it survives up to time t starting at 0. 
Component does not “age”.

Assuming[s>0&&t>0,Probability[X>s+t\[Conditioned]X>s,X\[Distributed]ExponentialDistribution[λ]]]

Matrix Exponential Formulation

Queueing paradigm: A customer enters the system and sequentially visits a subset of the servers before leaving the system

s0="Enter";
s1="State 1\nRate (Subscript[μ, 1])";
si="State i\n Rate (Subscript[μ, i])";
sj="State j\n Rate (Subscript[μ, j])";
sm="State m\n Rate (Subscript[μ, m])";
sl="Leave";
t01="(Subscript[p, 1])";
t0i="(Subscript[p, i])";
t0j="(Subscript[p, j])";
t0m="(Subscript[p, m])";
tjl="(Subscript[q, j])";
til="(Subscript[q, i])";
tml="(Subscript[q, m])";
tij="(Subscript[p, ij])";
t1j="(Subscript[p, 1  j])";

{s0,s1,si,sj,sm,sl}=Map[Style[#,Bold,Red]&,{s0,s1,si,sj,sm,sl}];
{t01,t0i,t0j,t0m,tjl,til,tml,tij,t1j}=Map[Style[#,Bold,Blue]&,{t01,t0i,t0j,t0m,tjl,til,tml,tij,t1j}];

lgp={{s0->s1,t01},{s0->si,t0i},{s0->sj,t0j},{s0->sm,t0m},{s1->sj,t0j},{si->sj,tij},{sj->sl,tjl},{sm->sl,tml}};

LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05]

The “magic” Matrix Exponential” (ME) method

Define the rate (M) matrix and the transition probability (P) matrices:
M=Diagonal[{(Subscript[μ, 1,])(Subscript[μ, 2]),...,(Subscript[μ, m])}],
 P=SparseArray[{...,{i,j}->(Subscript[p, ij]),...}]
 I=IdentityMatrix[m],
 p={(Subscript[p, 1]), (Subscript[p, 2]),...,(Subscript[p, m])}, where (Subscript[p, i]) = probability of starting in state i  (t=0),
 q={(Subscript[q, 1]), (Subscript[q, 2]),...(Subscript[q, m])}, where (Subscript[q, j]) = probability of leaving system from state j,
 ϵ={1, 1,...,1}.

Compute B=M.(I-P) .

Then the probability of still being in the system at time t, having started in state i with probability (Subscript[p, i]) at t=0, is given by

p.(e^(-t B)).ϵ

In particular, the probability of still being in the system at time t, having started in state 1 at time t=0, is given by the sum of all of the elements of the first row of (e^(-t B)). This also means that 1 minus the sum of the elements in the first row is the probability of leaving the system by time t.

Cyclic Markov chain  example

No “absorbing” states

Repairable component

s0="Good\nRate λ";
s1="Under repair\nRate μ";
t01="1";
t10="1";
{s0,s1}=Map[Style[#,Bold,Red]&,{s0,s1}];
{t01,t10}=Map[Style[#,Bold,Blue]&,{t01,t10}];
lgp={{s0->s1,t01},{s1->s0,t10}};

LayeredGraphPlot[lgp,
Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->"Repairable component"]

The ME method applied to the repairable component example

Define the M matrix (component failure and repair rates)

m2Repair=DiagonalMatrix[{λ,μ}];
MatrixForm[m2Repair]

Define the P matrix (state transition probabilities)

p2Repair={{0,1},{1,0}};
MatrixForm[p2Repair]

Compute the B matrix (system process rates)

b2Repair=m2Repair.(IdentityMatrix[2]-p2Repair);
MatrixForm[b2Repair]

Compute the ME distribution matrix

sol2Repair=MatrixExp[-t b2Repair];
MatrixForm[sol2Repair]

Interpretation of the entries of the ME distribution matrix

ss11=Assuming[λ>0&&μ>0,Limit[sol2Repair[[1,1]],t->∞]];
ss12=Assuming[λ>0&&μ>0,Limit[sol2Repair[[1,2]],t->∞]];
ss21=Assuming[λ>0&&μ>0,Limit[sol2Repair[[2,1]],t->∞]];
ss22=Assuming[λ>0&&μ>0,Limit[sol2Repair[[2,2]],t->∞]];
mEvents={{"Being in Good state \n having started in Good state","Being in Under repair state \n having started in Good state"},{"Being in Good state \n having started in Under repair state","Being in Under repair state \n having started in Under repair state"}};Grid[Prepend[Transpose[{Flatten[mEvents,1],Flatten[sol2Repair,1],Flatten[{{ss11,ss12},{ss21,ss22}},1]}],Map[Text,{"Event","Probability at time t","Steady-state probability"}]],Frame->All,Background->LightYellow]

State occupancy plots

Plot[{sol2Repair[[1,1]]/.{λ->0.1,μ->0.2},sol2Repair[[1,2]]/.{λ->0.1,μ->0.2}},{t,0,20},AxesLabel->{"Time","Probability of being in state:"},PlotLabel->Text["Repairable Component state occupancy probabilities:\nStarting in Good state and λ=0.1, μ=0.2"],Epilog->{Inset[Text["Good"],{3,0.88}],Inset[Text[sol2Repair[[1,1]]],{3,0.68}],Inset[Text["Under repair"],{3,0.15}],Inset[Text[sol2Repair[[1,2]]],{3,0.3}],Inset[Text["Good steady state"],{14,0.75}],Inset[Text[ss11],{18,0.75}],Inset[Text["Under repair steady state"],{13,0.25}],Inset[Text[ss12],{18,0.25}]},ImageSize->600]

Markov chain examples with absorbing states

Simplex system (1 component)

One component: System fails when component fails.

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])";
s2="All failed";
t01="1";
t12="1";
{s0,s1,s2}=Map[Style[#,Bold,Red]&,{s0,s1,s2}];
{t01,t12}=Map[Style[#,Bold,Blue]&,{t01,t12}];
lgp={{s0->s1,t01},{s1->s2,t12}};
LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->Text["Simplex System (1 component)"]]

Parallel system (3 components)

System fails when all 3 components fails (3-out-of-3F system)

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])";
s2="State 2 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])";
s3="State 3 \n Rate (Subscript[λ, 1])+(Subscript[λ, 3])";
s4="State 4 \n Rate (Subscript[λ, 2])+(Subscript[λ, 3])";
s5="State 5 \n Rate (Subscript[λ, 1])";
s6="State 6 \n Rate (Subscript[λ, 2])";
s7="State 7 \n Rate (Subscript[λ, 3])";
s8="All failed";

t01="1";
t12="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t13="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t14="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t25="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2]))";
t26="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2]))";;
t35="\n\n(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 3]))";;
t37="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 3]))";;
t46="\n\n(Subscript[λ, 3])/((Subscript[λ, 2])+(Subscript[λ, 3]))";;
t47="(Subscript[λ, 2])/((Subscript[λ, 2])+(Subscript[λ, 3]))";;
t58="(Subscript[λ, 1])/(Subscript[λ, 1])";
t68="(Subscript[λ, 2])/(Subscript[λ, 2])";
t78="(Subscript[λ, 3])/(Subscript[λ, 3])";

{s0,s1,s2,s3,s4,s5,s6,s7,s8}=Map[Style[#,Bold,Red]&,{s0,s1,s2,s3,s4,s5,s6,s7,s8}];
{t01,t12,t13,t14,t25,t26,t35,t37,t46,t47,t58,t68,t78}=Map[Style[#,Bold,Blue]&,{t01,t12,t13,t14,t25,t26,t35,t37,t46,t47,t58,t68,t78}];

lgp={{s0->s1,t01},{s1->s2,t12},{s1->s3,t13},{s1->s4,t14},{s2->s5,t25},{s2->s6,t26},{s3->s5,t35},{s3->s7,t37},{s4->s6,t46},{s4->s7,t47},{s5->s8,t58},{s6->s8,t68},{s7->s8,t78}};
LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->"Parallel System (3 components)"]

Series system (3 components)

System fails when at least one component fails (1-out-of-3F system)

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])";
s2="One component failed";
t01="1";
t121="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t122="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t123="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";

{s0,s1,s2}=Map[Style[#,Bold,Red]&,{s0,s1,s2}];
{t01,t121,t122,t123}=Map[Style[#,Bold,Blue]&,{t01,t121,t122,t123}];

lgp={{s0->s1,t01},{s1->s2,t121},{s1->s2,t122},{s1->s2,t123}};

LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->"Series System (3 components)"]

Triple Modular Redundant (TMR) system

System fails when at least 2 components fail (2-out-of-3F system)

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])";
s2="State 2 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])";
s3="State 3 \n Rate (Subscript[λ, 1])+(Subscript[λ, 3])";
s4="State 4 \n Rate (Subscript[λ, 2])+(Subscript[λ, 3])";
s5="Two failed";

t01="1";
t12="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t13="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t14="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";

t121="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t122="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t123="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t251="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2]))";
t252="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2]))";
t351="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 3]))";
t352="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 3]))";
t451="(Subscript[λ, 2])/((Subscript[λ, 2])+(Subscript[λ, 3]))";
t452="(Subscript[λ, 3])/((Subscript[λ, 2])+(Subscript[λ, 3]))";

{s0,s1,s2,s3,s4,s5}=Map[Style[#,Bold,Red]&,{s0,s1,s2,s3,s4,s5}];
{t01,t12,t13,t14,t251,t252,t351,t352,t451,t452}=Map[Style[#,Bold,Blue]&,{t01,t12,t13,t14,t251,t251,t351,t352,t451,t452}];

lgp={{s0->s1,t01},{s1->s2,t12},{s1->s3,t13},{s1->s4,t14},{s2->s5,t251},{s2->s5,t252},{s3->s5,t351},{s3->s5,t352},{s4->s5,t451},{s4->s5,t452}};
LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->"Triple Modular Redundant (TMR) System"]

Triple modular redundant-simplex system

Similar to 2-out-of-3 fail system, except when first component fails, both the failed component and one of the unfailed components is discarded. System fails when second failure occurs.

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])";
s2="State 2 \n Rate (Subscript[λ, 1])";
s3="State 3 \n Rate (Subscript[λ, 2])";
s4="State 4 \n Rate (Subscript[λ, 3])";
s5="Two failed";

t01="1";
t12="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t13="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t14="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";

t25="(Subscript[λ, 1])/(Subscript[λ, 1])";
t35="(Subscript[λ, 2])/(Subscript[λ, 2])";
t45="(Subscript[λ, 3])/(Subscript[λ, 3])";

{s0,s1,s2,s3,s4,s5}=Map[Style[#,Bold,Red]&,{s0,s1,s2,s3,s4,s5}];
{t01,t12,t13,t14,t25,t35,t45}=Map[Style[#,Bold,Blue]&,{t01,t12,t13,t14,t25,t35,t45}];

lgp={{s0->s1,t01},{s1->s2,t12},{s1->s3,t13},{s1->s4,t14},{s2->s5,t25},{s3->s5,t35},{s4->s5,t45}};
LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->"TMR-Simplex System"]

The “magic” Matrix Exponential (ME) method -- parallel system example

Define component rates matrix M, transition probablities matrix P and system process rate matrix B:
(Note: The “extra” state 8 is the absorbing state.)

mMatrix=DiagonalMatrix[{μ1,μ2,μ3,μ4,μ5,μ6,μ7,0}];
pMatrix=SparseArray[{{1,2}->p12,{1,3}->p13,{1,4}->p14,{2,5}->p25,{2,6}->p26,{3,5}->p35,{3,7}->p37,{4,6}->p46,{4,7}->p47,{5,8}->p58,{6,8}->p68,{7,8}->p78,{8,8}->1}];
bMatrix=mMatrix.(IdentityMatrix[8]-pMatrix);
meMatrix=MatrixExp[-t bMatrix];
Grid[Transpose[{Map[Text,{"State holding times matrix M","Transition Prob matrix P","System process rate matrix B=M.(I-P)"}],Map[MatrixForm,{mMatrix,pMatrix,bMatrix}]}],Frame->All]

Extract the distribution corresponding to being in the absorbing state 8 (“All failed”) at time t, having started in the initial state 1 (“All good”) at time t=0. It is in position (1,8).
Substitute values for the transition probabilities in terms of λ’s and μ’s.
Subtract from 1 to get the system survival function:

Simplify[1-meMatrix[[1,8]]/.{p12->λ3/μ1,p13->λ2/μ1,p14->λ1/μ1,p25->λ2/μ2,p26->λ1/μ2,p35->λ3/μ3,p37->λ1/μ3,p46->λ3/μ4,p47->λ2/μ4,p58->1,p68->1,p78->1}/.{μ1->λ1+λ2+λ3,μ2->λ1+λ2,μ3->λ1+λ3,μ4->λ2+λ3,μ5->λ1,μ6->λ2,μ7->λ3}]

Sanity check  -- should be equal to 1 minus the probability that all three components fail before time t. Check!

1-Expand[(1-(E^(-t λ1)))(1-(E^(-t λ2)))(1-(E^(-t λ3)))]

System survival function calculations

mSimplex=DiagonalMatrix[{μ1,0}];
pSimplex={{0,p12},{0,1}};
bSimplex=mSimplex.(IdentityMatrix[2]-pSimplex);
solSimplex=Simplify[1-MatrixExp[-t bSimplex][[1,2]]/.{p12->1}/.{μ1->λ1}]

mParallel=DiagonalMatrix[{μ1,μ2,μ3,μ4,μ5,μ6,μ7,0}];
pParallel={{0,p12,p13,p14,0,0,0,0},{0,0,0,0,p25,p26,0,0},{0,0,0,0,p35,0,p37,0},{0,0,0,0,0,p46,p47,0},{0,0,0,0,0,0,0,p58},{0,0,0,0,0,0,0,p68},{0,0,0,0,0,0,0,p78},{0,0,0,0,0,0,0,1}};
bParallel=mParallel.(IdentityMatrix[8]-pParallel);
solParallel=Simplify[1-MatrixExp[-t bParallel][[1,8]]/.{p12->λ3/μ1,p13->λ2/μ1,p14->λ1/μ1,p25->λ2/μ2,p26->λ1/μ2,p35->λ3/μ3,p37->λ1/μ3,p46->λ3/μ4,p47->λ2/μ4,p58->1,p68->1,p78->1}/.{μ1->λ1+λ2+λ3,μ2->λ1+λ2,μ3->λ1+λ3,μ4->λ2+λ3,μ5->λ1,μ6->λ2,μ7->λ3}]

mTMRSimplex=DiagonalMatrix[{μ1,μ2,μ3,μ4,0}];
pTMRSimplex={{0,p12,p13,p14,0},{0,0,0,0,p25},{0,0,0,0,p35},{0,0,0,0,p45},{0,0,0,0,1}};
bTMRSimplex=mTMRSimplex.(IdentityMatrix[5]-pTMRSimplex);
solTMRSimplex=Simplify[1-MatrixExp[-t bTMRSimplex][[1,5]]/.{p12->λ3/μ1,p13->λ2/μ1,p14->λ1/μ1,p25->1,p35->1,p45->1}/.{μ1->λ1+λ2+λ3,μ2->λ1,μ3->λ3,μ4->λ2}]

mSeries=DiagonalMatrix[{μ1,0}];
pSeries={{0,p12},{0,1}};
bSeries=mSeries.(IdentityMatrix[2]-pSeries);
solSeries=Simplify[1-MatrixExp[-t bSeries][[1,2]]/.{p12->1}/.{μ1->λ1+λ2+λ3}]

mTMR=DiagonalMatrix[{μ1,μ2,μ3,μ4,0}];
pTMR={{0,p12,p13,p14,0},{0,0,0,0,p25},{0,0,0,0,p35},{0,0,0,0,p45},{0,0,0,0,1}};
bTMR=mTMR.(IdentityMatrix[5]-pTMR);
solTMR=Simplify[1-MatrixExp[-t bTMR][[1,5]]/.{p12->λ3/μ1,p13->λ2/μ1,p14->λ1/μ1,p25->1,p35->1,p45->1}/.{μ1->λ1+λ2+λ3,μ2->λ1+λ2,μ3->λ1+λ3,μ4->λ2+λ3}]

System survival functions

Note that all of the survival (reliability) functions, except for TMR-simplex, are symmetric in their (Subscript[λ, i]) parameters. Symmetric means formula is unchanged by any permutation of the indices of the (Subscript[λ, i])’s.

subs1={λ1->λ,λ2->λ,λ3->λ};
subs2={λ t->x};
Needs["PlotLegends`"];
Grid[Transpose[{Map[Text,{"Parallel","TMR-Simplex","Simplex","TMR (2-out-of-3)","Series"}],Map[Text,{solParallel,solTMRSimplex,solSimplex,solTMR,solSeries}]}],Frame->All,Background->LightBlue]

System survival function plots

Plot[{solParallel/.subs1/.subs2,solTMRSimplex/.subs1/.subs2,solSimplex/.subs1/.subs2,solTMR/.subs1/.subs2,solSeries/.subs1/.subs2},{x,0,5},AxesLabel->{"λ t","Reliability"},PlotStyle->{Red,Blue,Green,Purple,Brown},PlotLabel->Text["Survival Function Comparison"],PlotLegend->{"Parallel","TMR-Simplex","Simplex","TMR (2/3)","Series"},LegendPosition->{0.1,-0.3},ImageSize->600]

Alternative approach: Sum over all paths

Use the ME method to compute the survival function associated with each individual path through the system state space from initial state to final (absorbing) state. Then sum these individual path functions over all paths. Can do this because they are mutually exclusive events.

This approach is more general in the sense that it is possible to distinguish paths which lead to system failure from those which don’t, in case there are any “sequence-dependent” failure sequences.

This approach may be more suitable for complex systems in which the “absorbing” states are not so obvious, and must be determined by an external simulation connected to the reliability execution engine.

s0="All good";
s1="State 1 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])";
s2="State 2 \n Rate (Subscript[λ, 1])+(Subscript[λ, 2])";
s3="State 3 \n Rate (Subscript[λ, 1])+(Subscript[λ, 3])";
s4="State 4 \n Rate (Subscript[λ, 2])+(Subscript[λ, 3])";
s5="State 5 \n Rate (Subscript[λ, 1])";
s6="State 6 \n Rate (Subscript[λ, 2])";
s7="State 7 \n Rate (Subscript[λ, 3])";
s8="All failed";

t01="1";
t12="(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t13="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t14="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3]))";
t25="(Subscript[λ, 2])/((Subscript[λ, 1])+(Subscript[λ, 2]))";
t26="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 2]))";;
t35="\n\n(Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 3]))";;
t37="(Subscript[λ, 1])/((Subscript[λ, 1])+(Subscript[λ, 3]))";;
t46="\n\n(Subscript[λ, 3])/((Subscript[λ, 2])+(Subscript[λ, 3]))";;
t47="(Subscript[λ, 2])/((Subscript[λ, 2])+(Subscript[λ, 3]))";;
t58="(Subscript[λ, 1])/(Subscript[λ, 1])";
t68="(Subscript[λ, 2])/(Subscript[λ, 2])";
t78="(Subscript[λ, 3])/(Subscript[λ, 3])";

{s0,s1,s2,s3,s4,s5,s6,s7,s8}=Map[Style[#,Bold,Red]&,{s0,s1,s2,s3,s4,s5,s6,s7,s8}];
{t01,t12,t13,t14,t25,t26,t35,t37,t46,t47,t58,t68,t78}=Map[Style[#,Bold,Blue]&,{t01,t12,t13,t14,t25,t26,t35,t37,t46,t47,t58,t68,t78}];

lgp={{s0->s1,t01},{s1->s4,t14},{s4->s7,t47},{s7->s8,t78}};
LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->Text["Parallel System (3 components)\nOne path only:\nComponent 1 fails->Component 2 fails->Component 3 fails"]]

Parallel system single path calculations

Pick a “typical” path from “All good” to “All failed” -- say,
Component 1 fails -> Component 2 fails -> Component 3 fails.

Define the M, P, and B matrices

mSinglePath=DiagonalMatrix[{μ1,μ2,μ3,0}];
pSinglePath=SparseArray[{{1,2}->p12,{2,3}->p23,{3,4}->p34,{4,4}->1}];
bSinglePath=mSinglePath.(IdentityMatrix[4]-pSinglePath);
meDistSinglePath=MatrixExp[-t bSinglePath];
Grid[{Map[Text,{"Rate matrix M","Transition Prob matrix P","System process rate matrix B=M.(I-P)"}],Map[MatrixForm,{mSinglePath,pSinglePath,bSinglePath}]},Frame->All]

Extract the unreliability function for this path as the (1,4) entry in the ME distribution matrix.
This is the probability of entering the absorbing state (“All failed”) before time t, given that initial state was “All good” at time t=0, along the specifc path.

meDistSinglePath=MatrixExp[-t bSinglePath];sf=meDistSinglePath[[1,4]]/.{p12->λ1/μ1,p23->λ2/μ2,p34->λ3/μ3}/.{μ1->λ1+λ2+λ3,μ2->λ2+λ3,μ3->λ3}

Now sum over all paths. Do this simply by summing this formula over all 6 permutations of the rates {λ1, λ2, λ3}.

permList=Permute[{λ1,λ2,λ3},SymmetricGroup[3]];
formulaList=sf/.{λ1->#[[1]],λ2->#[[2]],λ3->#[[3]]}& /@ permList;
formulaTotal=Total[formulaList];
formulas=AppendTo[formulaList,FullSimplify[formulaTotal]];
Grid[Transpose[{Map[Text,{"Path 1","Path 2","Path 3","Path 4","Path 5","Path 6","Total"}],formulas}],Frame->All]

Compare answer to the original ME formulation. Check!

Simplify[1-formulaTotal==solParallel]

Symmetric polynomials

Traditional generating function for (complete) symmetric functions (Subscript[h, r]):

H(t)=(∑_(r=1)^n)(Subscript[h, r])(t^r)=(∏_(i=1)^n)(((1-(Subscript[μ, i])t))^(-1))

Example: generate complete symmetric functions in three indeterminates {μ1, μ2, μ3}.

h[t_,n_]:=(∏_(i=1)^n)(((1-(Subscript[μ, i])t))^(-1));
Table[Expand[SeriesCoefficient[h[t,3],{t,0,k}]],{k,0,5}]//Column

These same complete symmetric functions can be generated by

G(t)=μ1 μ2 ...μn(∑_(j=0)^∞)(((-1))^j)(Subscript[h, j])((t^(j+n))/((j+n)!))

where G is the unreliability function (1 - survival function) of an n-stage path with successive rates μ1, μ2, ...μn and transition probabilites all equal to 1.

Example for n=3. Note how the symmetric functions appear in the upper-right hand corners (position (1,4)) of the powers of B:

mMat=DiagonalMatrix[{μ1,μ2,μ3,0}];
pMat=SparseArray[{{1,2}->1,{2,3}->1,{3,4}->1,{4,4}->1}];
bMat=mMat.(IdentityMatrix[4]-pMat);
Expand[Table[-1/(μ1 μ2 μ3)MatrixPower[bMat,k][[1,4]],{k,3,8}]]//Column

Recall that

exp(-t B)=I-t B+(1/(2!))(((t B))^2)-(1/(3!))(((t B))^3)+...

Therefore, it is not surprising that the same complete symmetric functions appear as (part of) the coefficients of the series representing the path unreliability function (1 - survival function):

g[t_]:=MatrixExp[-t bMat][[1,4]];
Table[Together[(((-1))^j)(j+3)!/(μ1 μ2 μ3)SeriesCoefficient[g[t],{t,0,j+3}]],{j,0,5}]//Column

Another Markov chain example

From “Performance and Reliability Analysis of of Computer Systems”, 
Sahner, Trivedi, Puliafito, 
Kluwer Academic Publishers, 1996, pp. 74-76

System description

We wish to compute the probabilities of being in (absorbing) states 5 or 6 at time t, given that the process started in either state 1, 2, or 3, weighted by initial probabilities 0.4, 0.4, 0.2, respectively.

s0="Start";
s1="State 1 \n Rate (Subscript[λ, 1])";
s2="State 2 \n Rate (Subscript[λ, 2])";
s3="State 3 \n Rate (Subscript[λ, 3])+(Subscript[λ, 4])+(Subscript[λ, 5])";
s4="State 4 \n Rate (Subscript[λ, 6])";
s5="State 5\nAbsorbing";
s6="State 6\nAbsorbing";
t01="(Subscript[p, 1])";
t02="(Subscript[p, 2])";
t03="(Subscript[p, 3])";
t13="1";
t23="1";
t34="(Subscript[λ, 3])/((Subscript[λ, 3])+(Subscript[λ, 4])+(Subscript[λ, 5]))";
t35="(Subscript[λ, 4])/((Subscript[λ, 3])+(Subscript[λ, 4])+(Subscript[λ, 5]))";
t36="(Subscript[λ, 5])/((Subscript[λ, 3])+(Subscript[λ, 4])+(Subscript[λ, 5]))";
t45="1";

{s0,s1,s2,s3,s4,s5,s6}=Map[Style[#,Bold,Red]&,{s0,s1,s2,s3,s4,s5,s6}];
{t01,t02,t03,t13,t23,t34,t35,t36,t45}=Map[Style[#,Bold,Blue]&,{t01,t02,t03,t13,t23,t34,t35,t36,t45}];


lgp={{s0->s1,t01},{s0->s2,t02},{s0->s3,t03},{s1->s3,t13},{s2->s3,t23},{s3->s4,t34},{s3->s5,t35},{s3->s6,t36},{s4->s5,t45}};

LayeredGraphPlot[lgp,Left,VertexLabeling->True,AspectRatio-> 1/2,PlotRangePadding->0.05,PlotLabel->Text["Markov Chain Example"]]

Initial state probabilities :
 (Subscript[p, 1]) = 4/10, 
 (Subscript[p, 2]) = 4/10, 
 (Subscript[p, 3]) = 2/10;

Transition rates:
 (Subscript[μ, 1]) = (Subscript[λ, 1]) = 4, 
 (Subscript[μ, 2]) = (Subscript[λ, 2]) = 1,  
 (Subscript[μ, 3]) = (Subscript[λ, 3]) + (Subscript[λ, 4]) + (Subscript[λ, 5]) = 10, where (Subscript[λ, 3]) = 2, (Subscript[λ, 4]) = 5, (Subscript[λ, 5])= 3,
 (Subscript[μ, 4]) = (Subscript[λ, 6]) = 2;

Transition probabilities: 
(Subscript[p, 13]) = 1, 
(Subscript[p, 23]) = 1, 
(Subscript[p, 34]) = 2/10, 
(Subscript[p, 35]) = 5/10, 
(Subscript[p, 36]) = 3/10, 
(Subscript[p, 45]) = 1;

Survival function (reliability) calculation using the “magic” Matrix Exponential (ME) method:

Note that the solution in the text (cited above) takes one and one-half pages!

STPsubs={μ1->4,μ2->1,μ3->10,μ4->2,p13->1,p23->1,p34->2/10,p35->5/10,p36->3/10,p45->1};
pSTP={4/10,4/10,2/10,0,0,0};
ϵSTP={0,0,0,0,1,1};

pMatSTP={{0,0,p13,0,0,0},{0,0,p23,0,0,0},{0,0,0,p34,p35,p36},{0,0,0,0,p45,0},{0,0,0,0,1,0},{0,0,0,0,0,1}};
mMatSTP=DiagonalMatrix[{μ1,μ2,μ3,μ4,0,0}];

expMatSTP=MatrixExp[-t mMatSTP.(IdentityMatrix[6]-pMatSTP)];
formulaSTP=Collect[solSTP=(1-pSTP. expMatSTP.ϵSTP)/.STPsubs,{(E^(-10 t)),(E^(-4 t)),(E^(-2 t)),(E^(-t))},Simplify];

Plot[solSTP,{t,0,5},PlotRange->All,AxesLabel->{"Time","Reliability"},PlotLabel->Text["Markov Chain Example Survival Function"],Epilog->Inset[Text[formulaSTP],{1.7,0.6}]]