(* ::Package:: *) (*:Mathematica Version: 6.0 *) (*:Package Version: 0.4.1 *) (*:Name: Dependability`StateDiagrams` *) (*:Context: Dependability`StateDiagrams` *) (*:Title: *) (*:Author: Bjarne E. Helvik *) (*:History: Preliminary version issued as an exercise notebook with solutions in PhD (Dr. Ing.) course " P\[ARing]litelighet i telematikk og datamaskinsystem" at NTH, Trondheim Norway , 1996-10-01. First Package version (0.1) issued 1999-12-23. System composition functions included 2000-05-24 (version 0.2), Introduction of eigenvalues for computation of transient solutions 2003-07-07 (version 0.3) 2008-01-24 Inclusion of a new version of the PlotDiagram function based on version 6.0 functionality. This ha previous been separate, based on Combinatorica. 2008-11-13 Corrected cosmetic bug in Plotdiagram *) (*:Copyright: All rights reserved Bjarne E. Helvik. *) (*:Reference: Usage messages. The notebook Analyse-StateDiagrams.nb *) (*:Summary: This package evaluates dependability models based on continuous time discrete state Markov chains, i.e. so called state models. It also includes functionality to compose system models from state models of subsystems. Makes a plot for the state diagram. *) (*:Keywords: Dependability, Availability, Reliability, MTBF, MUT, MDT, MTFF *) (*:Requirements: No special system requirements. *) (*:Sources: Buzacott, "Markov Approach to finding Failure Times of repairable systems", IEEE Trans. on Reliability, Vol. R-19, No. 6, pp. 128-133, 1970 *) BeginPackage["Dependability`StateDiagrams`"] SetDiagonal::usage = "SetDiagonal[m] produces a complete transition matrix by introducing diagonal elements. m is an incomplete transition matrix where m[[j,i]] i=/=j is the transition rates from state i to j. NB! i is the matrix column no., j is the row no." ProbStationary::usage = "ProbStationary[m] determines the stationary state probability vector for the system defined by the complete transition matrix m." UnAvailability::usage = "UnAvailability[ps,working] computes the unavailability for a system with state probability vector ps and a corresponding Boolean vector working, where working[[i]]=True if i is a working state." ProbTransient::usage = "ProbTransient[m,P0,t] determines the transient state probabilities for a system at time t. The system is defined by the complete transition matrix m. P0 is the initial state probability vector at time 0." RelFunc::usage = "RelFunc[m, WorkingQ, t, Ini] determines the reliability function (the probability that the system has not failed in [0,t]) at time t. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state. The system is in state Ini at time t=0. The parameter Ini is optional with default value 1." MTFF::usage= "MTFF[m, WorkingQ, Ini] determines the Mean Time First Failure. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state. The system is in state Ini at time t=0. The parameter Ini is optional wit default value 1." MTTF::usage= "MTTF[m, WorkingQ] determines the Mean Time To Failure when the system is found working in a stationary state. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state." MUT::usage= "MUT[m, WorkingQ] determines the Mean Up Time. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state." MDT::usage= "MDT[m, WorkingQ] determines the Mean Down Time. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state." MTBF::usage= "MTBF[m, WorkingQ] determines the Mean Time Between Failures. The system is defined by the complete transition matrix m and the working states is given by the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state." KroenSum::usage= "KroenSum[M] where M = m1, m2, ..., mn yields the transition matrix of a system composed if independent subsystems 1 to n. mi is the complete transition matrix of subsystem i. Note the order of subsystems is significant and order used in functions KroenSum, LabelList, SeriesMode, ParallelMode and KoutofNMode must be consistent. (The mathematical operation performed by KroenSum is denoted the Kroenecker sum.)" LabelList::usage= "LabelList[L] where L = l1, l2, ..., ln yields a list of state labels for a system composed if independent subsystems 1 to n. li is the list of state labels for subsystem i. (li may be prefix with subsystem name to avoid ambiguity.) Note the order of subsystems is significant and order used in functions KroenSum, LabelList, SeriesMode, ParallelMode and KoutofNMode must be consistent." SeriesMode::usage= "SeriesMode[WQ] where WQ = WorkingQ1, WorkingQ2, ..., WorkingQn yields a Boolean vector where an element is True if the corresponding system state is working a state, otherwise False. The system is composed if independent subsystems 1 to n, forming a SERIES RELIABILITY STRUCTURE, i.e., all subsystems must work to have a working system. WorkingQi is a Boolean vector (list) where WorkingQi[[j]]=True if j is a working state of subsystem i. Note the order of subsystems is significant and order used in functions KroenSum, LabelList, SeriesMode, ParallelMode and KoutofNMode must be consistent." ParallelMode::usage= "ParallelMode[WQ] where WQ = WorkingQ1, WorkingQ2, ..., WorkingQn yields a Boolean vector where an element is True if the corresponding system state is working a state, otherwise False. The system is composed if independent subsystems 1 to n, forming a PARALLEL RELIABILITY STRUCTURE, i.e., at least one subsystem must work to have a working system. WorkingQi is a Boolean vector (list) where WorkingQi[[j]]=True if j is a working state of subsystem i. Note the order of subsystems is significant and order used in functions KroenSum, LabelList, SeriesMode, ParallelMode and KoutofNMode must be consistent." KoutofN::usage= "KoutofN[koutofn, k] generates a Boolean expression which is True if at least k out of the elements in the Boolean vector (list) koutofn is true." KoutofNMode::usage= "KoutofNMode[WQ, k] where WQ = WorkingQ1, WorkingQ2, ..., WorkingQn yields a Boolean vector where an element is True if the corresponding system state is working a state, otherwise False. The system is composed if independent subsystems 1 to n, forming a K-OUT-OF-N RELIABILITY STRUCTURE, i.e., at least k out of the subsystems must work to have a working system. WorkingQi is a Boolean vector (list) where WorkingQi[[j]]=True if j is a working state of subsystem i. Note the order of subsystems is significant and order used in functions KroenSum, LabelList, SeriesMode, ParallelMode and KoutofNMode must be consistent." MergeStates::usage= "MergeStates[M,i,j1,..,jn] produces the complete transition matrix of a model where states i, j1, .., jn are merged into one state. M is the complete transition matrix of original system. The return transitions into the 'unmerged' states are those of the original state i. NOTE: For a merge which does not alter the model, it is required that the transition rates out of each of the merged states into the rest of the system are identical." (***************** Other functions *****************) StepPlot::usage= "StepPlot[x_List] Creates a step plot at the specified points. The x_List must contains (x, y)-values." PlotDiagram::usage = "PlotDiagram[m, WorkingQ, Labels] Plots a dependability model represented by a state diagram. The parameters are the transition matrix m, the corresponding Boolean vector WorkingQ, where WorkingQ[[i]]=True if i is a working state, and the list of state labels denoted Labels. Labels must be stings or list of strings." (* *) Begin["`Private`"] (*Ensures that the "old" diagonal element is removed and ensures a row sum of zero *) SetDiagonal[m_] := Module[{mxx}, mxx = m (Table[1,{Length[m]},{Length[m]}] - IdentityMatrix[Length[m]]); mxx = mxx - (Transpose[mxx] . Table[1,{Length[m]}]) IdentityMatrix[Length[m]]] ProbStationary[m_] := Module[{mxx, bxx}, bxx = Append[Table[0,{Length[m]-1}],1]; mxx = Append[Drop[m,-1],Table[1,{Length[m]}]]; LinearSolve[mxx, bxx]] UnAvailability[ps_,working_] := ps . (If[#,0,1]& /@ working) (* The transient state probabilities are obtained by the built in mma routine DSolve. DSolve operates on set of explicitly defined equations. Hence these must first be defined from the transition matrix. Furthermore, the initial condition is defined as a set of explicit equations. The solution is given as a rule and is transformed to a linear vector. *) ProbTransient[m_,P0_,t_] := Module[{p, pt, eqn}, p = Table[pt[i][t],{i,Length[m]}]; eqn = Flatten[Append[ Inner[Equal ,m . p , D[p,t], List], Inner[Equal ,(p/.t->0) , P0, List]]]; p /. DSolve[eqn, p, t] //Simplify //Flatten] (* Alternative to ProbTransient by using eigenvalues. Not thourghly tested*) ProbTransientEigen[m_,P0_,t_]:=Module[{Vals,Vects, Ck}, (* a set of distinkt eienvalues are required *) {Vals,Vects}=Eigensystem[m]; (* Test if eigenvalues are distinct *) If[Length[Union[Vals]] < Length[Vals], Print["Repeated eigenvalues, Not implemented"]; Abort[]]; (* Finds constans by using the initial condition*) Ck=LinearSolve[Transpose[Vects],P0]; (* Generates the solution *) Ck. (Exp /@ (Vals * t) Vects )] RelFunc[m_, WorkingQ_,t_, Ini_:1] := Module[{M,Mx,My}, If[!WorkingQ[[Ini]], Print["Initial state is a down state"]; Return[]]; M = ArrangeMatrix[RotateLeft[ Transpose[RotateLeft[Transpose[m],Ini-1]],Ini-1], RotateLeft[WorkingQ,Ini-1]]; Mx = Append[M[[1,1]], Transpose[M[[2,1]]] . Table[1,{Length[M[[2,1]]]}]]; NullVec = Table[0, {Length[Mx]}]; My = Transpose[Append[Transpose[Mx],NullVec]]; NullVec[[1]] = 1; (1- Last[ProbTransient[My,NullVec,t]]) // Simplify ] (* Based on the paper "Buzacott [Markov Approach to finding Failure Times of repairable systems, IEEE Trans. on Reliability, Vol. R-19, No. 6, pp. 128-133]" which shows how the moments of the system times may be obtained from the transition matrix of a repairable system. *) (* Splits a matrix columnwise into two submatrixes, dependent on whether the states are up or down states *) SplitWrkng[mm_,WorkingQ_] := Module[{ColsW ={}, ColsF ={}}, index = Table[i,{i,Length[mm]}]; If[WorkingQ[[#]],ColsW=Append[ColsW,mm[[#]]], ColsF=Append[ColsF,mm[[#]]]]& /@ index; {ColsW ,ColsF}] (* Sorts the matrix so columns and rows for up and down states are collected so the up-state transitions (EX) are in the upper left corner and the down-state transitions are in the lower left (HX). Furthermore the matrix is diveded into the four submatrixes. Note that the paper operates with transposed matrixes.*) ArrangeMatrix[m_, WorkingQ_] := Module[{XY,XZ,EX,HX,GX,FX}, (* sorts columns *) XY = SplitWrkng[m,WorkingQ]; (* from the working states*) XZ = SplitWrkng[Transpose[XY[[1]]],WorkingQ]; EX = Transpose[XZ[[1]]]; GX = Transpose[XZ[[2]]]; (* from the failed states *) XZ = SplitWrkng[Transpose[XY[[2]]],WorkingQ]; FX = Transpose[XZ[[1]]]; HX = Transpose[XZ[[2]]]; {{EX,GX},{FX,HX}}] (* NB! Deal with an arbitrary initial state; default value is 1*) MTFF[m_, WorkingQ_, Ini_:1] := Module[{EX}, If[!WorkingQ[[Ini]], Print["Initial state is a down state"]; Return[]]; EX = ArrangeMatrix[RotateLeft[ Transpose[RotateLeft[Transpose[m],Ini-1]],Ini-1], RotateLeft[WorkingQ,Ini-1]][[1,1]]; - LinearSolve[EX, Prepend[Table[0,{Length[EX]-1}],1]] . Table[1,{Length[EX]}] ] MTTF[m_, WorkingQ_] := Module[{bz,pp, h1}, pp = SplitWrkng[ProbStationary[m],WorkingQ][[1]]; h1 = Table[1,{Length[pp]}]; bz = ArrangeMatrix[m,WorkingQ]; (h1 . Inverse[- bz[[1,1]]] . pp ) / (pp . h1) //Simplify] MUT[m_, WorkingQ_] := Module[{bz,pp, hw,hd}, pp = SplitWrkng[ProbStationary[m],WorkingQ][[1]]; bz = ArrangeMatrix[m,WorkingQ]; hw = Table[1,{Length[bz[[1,1]]]}]; hd = Table[1,{Length[bz[[2,2]]]}]; (hw . pp ) / (hd . bz[[2,1]] . pp) //Simplify] MDT[m_, WorkingQ_] := Module[{bx,pp, hw, hd}, pp = SplitWrkng[ProbStationary[m],WorkingQ][[1]]; bz = ArrangeMatrix[m,WorkingQ]; hw = Table[1,{Length[bz[[1,1]]]}]; hd = Table[1,{Length[bz[[2,2]]]}]; (1 - hw . pp ) / (hd . bz[[2,1]] . pp) //Simplify] MTBF[m_, WorkingQ_] := Module[{bx,pp, hw, hd}, pp = SplitWrkng[ProbStationary[m],WorkingQ][[1]]; bz = ArrangeMatrix[m,WorkingQ]; hd = Table[1,{Length[bz[[2,2]]]}]; 1 / (hd . bz[[2,1]] . pp) //Simplify] (*The Kroenecker product is most easily obtained by the outer product. The outer product yields a tensor which must be appropriately flattened to a two-dimentional structure. *) KroenSum[a_,b_] := Flatten[Transpose[Flatten[Transpose[ Outer[Times, a,IdentityMatrix[Length[b]]] + Outer[Times, IdentityMatrix[Length[a]], b], {4,1,3,2}], 1],{3,2,1}],1] ; KroenSum[a_,b_,c__] := KroenSum[KroenSum[a,b],c] LabelList[C__] := Flatten[Outer[List,C],Length[{C}]-1] SeriesMode[C__] := Apply[And, Flatten[Outer[List,C],Length[{C}]-1], {1,1}]; ParallelMode[C__] := Apply[Or, Flatten[Outer[List,C],Length[{C}]-1], {1,1}]; (* The function KoutofN uses the function KSubsets from < DiscreteMath`Combinatorica`>. The function is included in this package in order to reduce memory usage. Alternatively the Combinatorica package may be loaded. *) KSubsets[l_List,0] := { {} } KSubsets[l_List,1] := Partition[l,1] KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l]) KSubsets[l_List,k_Integer?Positive] := {} /; (k > Length[l]) KSubsets[l_List,k_Integer?Positive] := Join[ Map[(Prepend[#,First[l]])&, KSubsets[Rest[l],k-1]], KSubsets[Rest[l],k] ] (* < Or[x]}/.{List -> And} KoutofNMode[C__, k_] :=KoutofN[#,k]& /@ Flatten[Outer[List,C],Length[{C}]-1]; MergeStates[M_,i_,j__] := Module[{II,JJ}, II = IdentityMatrix[Length[M]]; II = Fold[ Drop,II,Map[List,{j}]]; (* In an "ideal merge" the return transitions from the merged states should, be identical. If not return transitions from state i dominates *) JJ= M .Transpose[ II]; (* Add entrance rate to the merged states *) II = IdentityMatrix[Length[M]]; II[[i]] = II[[i]] + Table[If[ MemberQ[{j},x],1,0],{x, Length[II[[i]]]}]; II = Fold[ Drop,II,Map[List,{j}]]; Transpose[ Transpose[JJ] .Transpose[ II]] ] (***************** Other functions *****************) StepPlot[x_List] := Module[{temp1, temp2, TempC}, {temp1, temp2} = Transpose[x]; TempC = Transpose[{Drop[temp1, 1], Drop[temp2, -1]}]; ListPlot[Sort[Join[x, TempC], (#1[[1]] < #2[[1]]) &], Joined -> True]] PlotDiagram[Mx_, WorkingQ_, Labels_] := Module[{My, transitions, full, final, sl}, (* ensures that composite labels are strings*) sl = ToString /@ Labels; (* sets diagonal elements to zero*) My = Mx Table[If[i==j,0,1], {i, Length[Mx]},{j,Length[Mx]}] (* This gives kosmetic failures with reals. My = Mx - DiagonalMatrix[Diagonal[Mx]]*); (*Build label to label transition*) transitions = Outer[(#2 -> #1) &, sl, sl]; (*Merges the input for Graphplot*) full = Flatten[ MapThread[{#1, #2} &, {transitions, My}, 2], {1, 2}]; (*Remove the null transitions*) final = Reverse[Complement[full, Select[full, (#[[2]] == 0) &]]]; GraphPlot[final,(* VertexLabeling->True,*) DirectedEdges -> True, VertexRenderingFunction -> ({White, EdgeForm[Black], If[WorkingQ[[Position[sl, #2][[1, 1]]]], Disk[#, .15], Rectangle[# - {0.15, 0.1}, # + {0.15, 0.1}]], Black, Text[#2, #1]} &)]] End[ ] EndPackage[ ]