Use of Mathematica for symbolic and numerical dependability analysis of systems modelled by state diagrams

Version 0.4
2008-01-24
© Bjarne E. Helvik

Introduction

This Mathematica notebook exemplifies the use of the package StateDiagrams.m for analysis of dependability models based on continuous time discrete state Markov chains, i.e. so called state models.  (A separate plotting package is also provided, see subsection below.) Assuming the package is installed on your computer, the package is loaded by

<<"StateDiagrams.m"

A list of the available functions in the package is obtained by the command:

Names["Dependability`StateDiagrams`*"]

As a general reference, a short description of a function is obtained by putting a question mark in front of its name. For instance:

?ProbTransient

The package may be used for symbolic as well as numerical analysis. In most of this notebook symbolic evaluation will be applied. The size of the system that can be evaluated, i.e. the number of states that can be handled, is dependent on whether a numerical or symbolic analysis is performed as well as the structure of the transition matrix. In general, numerical evaluation allows far lager systems.

Warning: This package is an early version and bugs may be present. There is no check on the correctness of the function arguments.

It is required that the user has some familiarity with dependability modelling and analysis. For an introduction see for instance [EHH00, Helv07].

Definition of a simple system model

The system modelled is described by a discrete state continuous time Markov model. It is implicitly assumed that the model has a finite number of states.

As an example, a model of a duplicated system with one repair facility is chosen. The system is simple enough to be easy to grasp, and at the same time non trivial. This system is dealt with in Section 7.2.2 of [EHH07] in the subsection on "Permanenthardware failures". The state diagram of the system is depicted below.

λ/:Positive[λ]:=True
μ/:Positive[μ]:=True

Head[λ]

For exemplification we introduce numerical parameter values for the failure, λ, and repair, μ, rates, where the scale is 1/hour.

Parameters = {λ -> 0.0038, μ -> 1.34};

Transition matrix

The transition matrix contains the transition rates between the states of the model. For the sake of convenience, the transposed matrix is first established, where element i,j gives the transition rate from state i to j.

Mtemp= {{,2 λ,0},{μ,xxx,λ},{0,μ,}};

Note that in this matrix, the diagonal elements are undefined or rubbish. We now have the incomplete transition matrix

Transpose[Mtemp] //MatrixForm

In general, the diagonal of the transition matrix may be derived from the other elements of the matrix. Hence, the complete transition matrix Mx may be obtained:

Chooses a simple two state model for exemplification

(Mx = SetDiagonal[Transpose[Mtemp]])//MatrixForm

Operational mode

The operational mode is given by a Boolean vector showing whether the system is in a working (UP) state or not.

WorkingQ = {True, True, False};

State labels

The state labels are given as a list of strings describing the states of the model. The sequence should correspond to the sequence above. State labels are optional. Example:

Labels = {"OK", "1 Failed", "2 Failed"};

Plot of the resulting model

Simple functionality to plot the state transition diagram is implemented in the StateDiagram package.

PlotDiagram[Mx,WorkingQ,Labels]

State probabilities and (un)availability

Stationary probabilities and stationary unavailability

The stationary probabilities (i.e. the  steady state probabilities) are obtained by solving the set of linear equations after the normalization constant is introduced into Mx.

Ps = ProbStationary[Mx] // Simplify

The unavailability is obtained by summing over the failed  states of the system.

UnAvailability[Ps,WorkingQ]

Numerical example

The numerical value for the unavailability may be obtained by introducing the parameters in the above closed for result, i.e.

NUnAvail= % /. Parameters

Alternatively, the entire computation may be carried out numerically, by first obtaining the numerical transition matrix. This is the advised approach when Mx is large.

Mx/. Parameters

and obtaining the unavailability as above:

UnAvailability[ProbStationary[%],WorkingQ]

Transient state probabilities and unavailability

The transient state probabilities for a system may be obtained by solving the differential equation ((dP(t))/(dt))==Mx.P(t)  for the system defined in the previous section with a given initial state probability, P0==P(0). Note that solving the transient state probabilities is more demanding than the stationary, and hence, symbolic evaluation is only feasible for small systems. Numerical evaluation extends the applicability to larger systems.

Assume that the system is in the OK state at t=0.

P0 ={1,0,0};

P[t_] =ProbTransient[Mx,P0,t]

%//FullSimplify

Checking the plausibility of the solution by summing up the state probabilities:

P[t] .{1,1,1} // Simplify

Comparison with the stationary state probabilities

If we assume that the system is stationary at t=0, i.e. P0 = Ps, we obviously obtain

ProbTransient[Mx,Ps,t]

which is equal to the stationary state probabilities as expected.

Transient unavailability

P[t] may be used to determine the transient unavailability.

UnAv[t_] = UnAvailability[P[t], WorkingQ] //FullSimplify

Plot[UnAv[t]/.Parameters,{t,0,(8/μ)/.Parameters}]

Interval availability during the first period of operation, with duration equal to one expected service restoration time

μ Integrate[UnAv[x], {x,0,1/μ}]/.Parameters //Simplify

Computation of system times

In this section,  the mean time to first failure (MTFF), mean time to failure (MTTF), mean up time (MUT), mean down time (MDT) and mean time between failures (MTBF) are obtained from the transition matrix, Mx,  of the example system. See [Buza70] for how all moments of the system times may be obtained from the transition matrix.

NB! The solutions are not optimized with respect to speed. Repetitive computations of the same quantities are carried out if more than one system time is determined.

Mean time to first failure, MTFF

Assuming that state 1 is the initial state, the mean time to failure is obtained as

MTFF[Mx, WorkingQ]

For comparison with other system times we also obtain its numerical value:

% /. Parameters

The system may be "started" in other states at t=0 by specifying the initial state. Let say we start the system with one server failed:

Labels[[2]]

MTFF[Mx, WorkingQ,2]

Note the negligible effect on the numerical value:

% /. Parameters

If the system is in a down state at t=0 the MTFF calculation is meaningless.

MTFF[Mx, WorkingQ,3]

Mean time to failure (MTTF)

This is the time from the system is found working  in a stationary state an until it fails. Hence, we should expect the MTFF to be between the two values for the MTFF determined above:

MTTF[Mx, WorkingQ]

% /. Parameters

Mean up time, MUT

Similarly the mean up time will be in the same magnitude as the quantities above when we have a high availability system

MUT[Mx, WorkingQ]

% /. Parameters

Note also that since the system is restored to state 2 after a system failure

Labels[[2]]

MUT is identical for the MTFF of the system initiated in state 2 found above.

Mean down time, MDT

Since the system has only one down state, the mean down time should be equal to the sojourn time of this state, i.e.:

MDT[Mx, WorkingQ]

% /. Parameters

Mean time between failures, MTBF)

MTBF[Mx, WorkingQ]

% /. Parameters

Check that the various system times are consistent:

MTBF[Mx, WorkingQ] ==MUT[Mx, WorkingQ]+MDT[Mx, WorkingQ]//Simplify

The reliability function

Let us regard the duplicated repaired example system modelled in the beginning of this notebook. 
We would like to obtain the reliability function (survival probability) of this system. First it is implicitly assumed that the everything OK state (or initial state) is state 1.

(Subscript[R, 1])[τ_]=RelFunc[Mx, WorkingQ,τ]

Comparison with the previously computed MTFF

The MTFF may be obtained by integrating the reliability function. It is seen that this yields the same result as above for well behaved system parameters.

Integrate[%,{τ,0, Infinity}]

% /. Parameters

Another initial state

The system may be started at an arbitrary operational state at t=0, for instance the state where only one processor is operational.

(Subscript[R, 2])[τ_]=RelFunc[Mx,WorkingQ,τ,2];

The two plots below compare these two reliability functions, where  (Subscript[R, 2])[τ] is represented by the dotted line. It is seen that initiating the system in state 2 gives a lower reliability.

Plot[Evaluate[{(Subscript[R, 1])[τ],(Subscript[R, 2])[τ]}/.Parameters],{τ,0,(1/10) MTFF[Mx,WorkingQ]/.Parameters},PlotStyle->{{Dashing[{1,1}]},Dashing[{0.015`,0.01`}]}]

Studying the behaviour close to  t=0 it is seen that lower reliability corresponds to a transition into the failed state before both processors become operational.

Plot[Evaluate[{(Subscript[R, 1])[τ],(Subscript[R, 2])[τ]}/.Parameters],{τ,0,(10/μ)/.Parameters},PlotStyle->{{Dashing[{1,1}]},Dashing[{0.015`,0.01`}]}]

Relation to Eigenvalues (Advanced topic)

It may be shown that for a system  that is repaired as long as the system is in an up state (working state) the reliability function may be closely approximated by (e^(Λ τ)) where Λ is the (nonzero, negative) eigenvalue of the system closest to zero.

This is illustrated below without dealing with the theory. Regard the example system, and make the system failed state absorbing by modifying the transition matrix of the system, i.e.:

Mxx =Mx; Mxx[[2,3]]=Mxx[[3,3]]=0;Mxx//MatrixForm

The eigenvalues for this system is

Eigenvalues[Mxx]

% /. Parameters

where it is seen that the absolute value of the third of these is much less than the others. (These eigenvalues is also recognized in the expression for (Subscript[R, 1])[τ] above.) The eigenvalues represent the inverse of the time constants of the system, and hence, as long as the absolute value of one eigenvalue is much less than the others, this will dominate long term behaviour of the system.

Λ = Last[Eigenvalues[Mxx]];
(Subscript[R, E])[τ_]= Exp[Λ τ ]

The approximation is notable only for times very close to zero,

Plot[Evaluate[{(Subscript[R, 1])[τ],(Subscript[R, E])[τ]}/.Parameters],{τ,0,(10/μ)/.Parameters},PlotStyle->{{Dashing[{1,1}]},Dashing[{0.02`,0.01`,0.005`,0.01`}]}]

and reflects the overall system reliability behaviour very well.

Plot[Evaluate[{(Subscript[R, 1])[τ],(Subscript[R, E])[τ]}/.Parameters],{τ,0,2 MTFF[Mx,WorkingQ]/.Parameters},PlotStyle->{{Dashing[{1,1}]},Dashing[{0.02`,0.01`,0.005`,0.01`}]}]

The absolute value of this eigenvalue |Λ| is often referred to as the failure rate of the system, and its inverse is close to the MTFF (and other times to failure).

1/Abs[Λ] /. Parameters

MTFF[Mx,WorkingQ] /. Parameters

Building system models from subsystem models (advanced topic)

Preliminaries

State models of systems with a realistic number of elements and failure causes tends to become large. Building them becomes tedious and is fault prone. Hence, it is advantageous to be able to build models for "parts" of the system and automatically compose them into system models. If it is feasible to split the system into independent subsystems, this may be done by combining state diagrams of subsystems with structural models (Reliability block diagrams or Fault trees) on the subsystem level. See for instance [EHH07] for a presentation. However, in most cases there are dependencies between subsystems caused by failure rates in one subsystem being dependent on the state of another subsystem, common repair resources, etc.

The objective of this chapter is to present Mathematica functions which may be used to build state models for systems from the state models of independent subsystems and exemplify how these may be modified to include dependencies between subsystems.

Use of Krönecker sum and supplementary functions

Regard a system made up by a number of subsystems. Each of the subsystems  are defined according to the section "Definition of a simple system model".  Let x denote the transition matrix of subsystem X, Xl denote the list of state labels for subsystem X and Xm denote the Boolean vector of operational modes in subsystem X, i.e., Xm[[i]] = True when i is a working state of subsystem X.

For the time being, consider these subsystems as independent with respect to failures and repairs.

To exemplify the functions, a system with three subsystems is regarded. Subsystem A with two states, subsystem B with three states, and subsystem C with two states. The number of subsystems and states are chosen to keep the example readable. In this case we get twelve system states. The Array function is used to generate the transition matrix and lists for the subsystems, e.g.

Array[a,{2,2}] // MatrixForm

Transition matrix

The transition matrix of the complete system may be derived from the subsystem transition matrixes by the use of Krönecker algebra. For a short introduction see for instance [Lesa88]. The Krönecker sum of two transition matrixes yields the transition matrix of the composed system, e.g.,

KroenSum[Array[a,{2,2}],Array[c,{2,2}]]//MatrixForm

The function KroenSum takes an arbitrary number of arguments (transition matrixes).

Mabc = KroenSum[Array[a,{2,2}],Array[b,{3,3}],Array[c,{2,2}]];

Mabc //MatrixForm

Labels

The function LabelList composes the labels of the subsystem states into a labels of the corresponding system states, e.g.,

TableForm[
LabelList[Array[Al,2],Array[Bl,3] ,Array[Cl,2]], 
TableDepth->1, TableHeadings->Automatic]

For instance, in system state 11, subsystem A is in state 2, subsystem B is in state 3 and subsystem C is in state 1.

System structure

The series structure requires that all subsystems must be in a working state if the corresponding system state shall be working, i.e.,

SeriesMode[Array[Am,2],Array[Bm,3],Array[Cm,2]]

The parallel structure requires that at least one of the subsystems must be in a working state if the corresponding system state shall be working, i.e.,

ParallelMode[Array[Am,2],Array[Bm,3],Array[Cm,2]]

The 2 out of N structure requires that at least two out of the subsystems must be in a working state if the corresponding system state shall be working, i.e.,

KoutofNMode[Array[Am,2],Array[Bm,3],Array[Cm,2],2]

It is easily demonstrated that a N out of N system is identical to a series system and a 1 out of N system is identical to a parallel system, e.g.

KoutofNMode[Array[Am,2],Array[Bm,3],Array[Cm,2],3] == SeriesMode[Array[Am,2],Array[Bm,3],Array[Cm,2]]

Composition of two simple systems with common repair

(Mdual = KroenSum[Mx, Mx])//MatrixForm

The labels may be prefixed with system name to improve readability, e.g.

StringJoin["SysA:", #]& /@ Labels

DualLabels = LabelList[StringJoin["SysA:", #]& /@ Labels, StringJoin["SysB:", #]& /@ Labels];

WorkingDualQ=SeriesMode[WorkingQ,WorkingQ]

The states and the corresponding operational mode of this generated system model are shown in the table below. The state indexes are the states indexes of the new system.

TableForm[
{DualLabels,WorkingDualQ} , 
TableDepth ->2, TableDirections->{Row, Column},
TableHeadings->{{"State description", "Working ?"}, Automatic}]

It is seen from the transition matrix, as well as intuitively, that in the states where both system are failed, i.e. 5, 6, 8, and 9, both systems are repaired. This violates the assumption of just one repair facility. Hence, the repair rates in these states must be modified. We may choose to 
•	maintain symmetry and reduce the rates to μ/2, or
•	give priority to system A, and set the repair rate of system B to zero when A has failed.
In this example the latter is chosen. Note that it is necessary to recompute the diagonal. (Keep in mind that the syntax in the Mathematica indexing is Mdual[[<to-state>, <From-state>]].)

Mdual[[4,5]] = Mdual[[5,6]] = Mdual[[7,8]]=Mdual[[8,9]]=0;

(Mdual = SetDiagonal[Mdual])//MatrixForm

We may now for instance compute the numerical value of the unavailability of the system.

UnAvailability[ProbStationary[Mdual /. Parameters],WorkingDualQ]

This is slightly higher than the unavailability of a series structure of the two systems

1-(1-NUnAvail)^2

due to the repair limitations. The unavailability of a system, NUnAvail, is determined above.

Composition of a HW and SW a subsystem

The "simple system" used as an example above is an ideal duplicated system with one repair facility. Let say that this processor run synchronously.  This may be regarded as the hardware subsystem of a system. We like to add a model of the software subsystem. For the sake of simplicity we use the subsystem model presented below.

(Msw = SetDiagonal[Transpose[{{0, γ},{ν,0}}]])//MatrixForm

SWlabels ={"Working SW", "Crashed SW"};

WorkingSWQ={True,False};

PlotDiagram[Msw,WorkingSWQ,SWlabels]

The software behaves identically whether it runs on one or two synchronous processors. (For the sake of keeping the example simple, we omit the important modelling issues concerning the change of configuration.) Hence, the HW and SW subsystems are almost independent, and we may generate a composite model.

(Msys = KroenSum[Mx, Msw])//MatrixForm

SYSlabels = LabelList[Labels, SWlabels];

SYSlabels//FullForm

??LabelList

WorkingSYSQ = SeriesMode[WorkingQ, WorkingSWQ];

The states and the corresponding operational mode of this generated system model are shown in the table and figure below. The state indexes are the states indexes of the new system.

TableForm[
{SYSlabels,WorkingSYSQ} , 
TableDepth ->2, TableDirections->{Row, Column},
TableHeadings->{{"State description", "Working ?"}, Automatic}]

PlotDiagram[Msys,WorkingSYSQ,SYSlabels]

Msys//MatrixForm

SYSlabels

It is obvious that the when both processors have failed, it is uninteresting whether the software works or not. Hence, we may "merge" the last two states. However, this will not be a mathematically correct merge since the repairs out of the two states differ. However, it is physically correct to let the return into the "1 Failed, Working Sw" state be the return transition from the merged state.

(MsysX = MergeStates[Msys, 5,6])// MatrixForm

(The function MergeStates allows us to merge an arbitrary number of states, cf.)

?MergeStates

WorkingSYSQ= Drop[WorkingSYSQ,{6}]

SYSlabels = Drop[SYSlabels,{6}]; SYSlabels[[5]] = First[SYSlabels[[5]]];

The table below displays the final transition matrix as a table with state descriptions as labels in order to improve readability.

TableForm[ MsysX, TableHeadings->{SYSlabels, SYSlabels} ]

And the corresponding diagram

PlotDiagram[MsysX,WorkingSYSQ,SYSlabels]

The unavailability, as exemplified below, and other dependability measures may now be obtained.

UnAvailability[ProbStationary[MsysX],WorkingSYSQ] //Simplify

References

[Buza70] Buzacott, "Markov approach to finding failure times of repairable systems", IEEE Trans. on Reliability, Vol. R-19, No. 6, pp. 128-133, 1970

[EHH07] Peder J. Emstad, Poul E. Heegaard, Bjarne E. Helvik, Laurent Paquereau, “Dependability and Performance in Information and Communication Systems; Fundamentals”, Lecture notes, Department of Telematics, NTNU, July 2007.

[Helv07] Bjarne E. Helvik: "Dependable Computing Systems and Communication Networks; Design and Evaluation", Draft  lecture notes - version per December 2007, Dept. of Telematics, NTNU. Tapir.

[Lesa88] A. Lesanovsky,  "Multistate Markov Models for systems with dependent units",  IEEE Trans. On Reliability, Vol. 37, No.5, pp. 505 - 511, Dec.1988