An application of graph theory in Markov chains reliability analysis

STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE An Application of Graph Theory in Markov Chains Reliability Analysis Pavel SKALNY Department of Applied Mathematics, Faculty of Electrical Engineering and Computer Science, VSB–Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava-Poruba, Czech Republic pavel.skaln[email protected] Abstract. The paper presents reliability analysis which was realized for an industrial company. The aim of the paper is to present the usage of discrete time Markov chains and the flow in network approach. Discrete Markov chains a well-known method of stochastic modelling describes the issue. The method is suitable for many systems occurring in practice where we can easily distinguish various amount of states. Markov chains are used to describe transitions between the states of the process. The industrial process is described as a graph network. The maximal flow in the network corresponds to the production. The Ford-Fulkerson algorithm is used to quantify the production for each state. The combination of both methods are utilized to quantify the expected value of the amount of manufactured products for the given time period. Keywords Discrete time Markov chain, flow in network, reliability. 1. Introduction The reliability of production plays the fundamental role in an industrial sphere. Nowadays the reliability of industry process is on a high level. It is increased by improving the quality of each component or by redundancy of the production process. Even though it is the top reliability process, there is still a chance that system fails. In our case we analyse the process which has no redundancy. Thus the information about the probability of the systems failures is very valuable. In the previous work an [1] a reliability analysis of the part of an industry process was realized. In the former application there was analysed a part of the industry process distributed in parallel. The main aim of the cited paper was to estimate the probability, that the production will be equal or greater than the industrial partner demand. To deal with the task a Monte Carlo simulation of Discrete time Markov chains was used. The goal of this paper is to present the Discrete time Markov chain analysis realized on the system, which is not distributed in parallel. Since the analysed process is more complicated the graph theory approach seems to be an appropriate tool. Graph theory has previously been applied to reliability, but for different purposes than we intend. One of example of application of graph theory in reliability is a reliability polynomial [3] or network reliability. In this approach a graph describes a network where each component (edge) has the same probability of fail. The most related work to our paper is probably the research of Christopher Dabrowski [4], [5]. As well as in our paper Dabrowski uses graph theory as a tool for counting with discrete time Markov chains. In contrast to our work the Markov chains described in Dabrowski one are used for different purposes. Dabrowski analyses a large scale grid. His study describe the application for finding states of the grid, which could lead to the system degradation. Dabrowski uses the graph theory to describe transitions of the Markov chain model between the initial state and the absorbing state. We use the graph theory to describe the analysed process on which the discrete time Markov chain is be applied. 2. Markov Chains Markov chain is a random process with a discrete time set T⊂N∪ {0}, which satisfies the so called “Markov property”. The Markov property means that the future evolution of the system depends only on the current state of the system and not on its past history. P{Xn+1 =xn+1|X0=x0, . . . , Xn=xn}= =P{Xn+1 =xn+1|Xn=xn},(1) c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 154 STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE where X1, . . . , Xnis a sequence of random variables. The index denotes certain time t∈T x1,...xnis a sequence of states in time t∈T. As a transition probability pij we regard probability, that the system changes from the state ito the state j: pij =P{Xn+1 =xj|Xn=xi}.(2) Matrix P, where pij is placed in row iand column j, is for all admissible iand jcalled transition probability matrix: P=     p11 p12 · · · p1n p21 p22 · · · p2n . . .. . .. . .. . . pn1pn2· · · pnn      .(3) Clearly all elements of the matrix Psatisfy the following property: ∀i∈ {1,2,· · · , n}: n X j=1 pij = 1.(4) As a probability vector we will understand a vector ~v where, an amount of its elements viare equal to the amount of states in Markov chain and the following equations holds: ∀vi∈~v :vi≤1, and X alli = 1.(5) By ~v(t)we denote vector of probabilities of all the states in time t∈T. By ~v(0) we denote an initial probability vector. Usually, all its coordinates are equal to zero except the first, which is equal to 1. The vector ~v(0) denotes the probability in which state the system occurs in time t= 0. It is easy to proof that [6]: ~v(t) = ~v(0) ·Pt.(6) As a stationary distribution we will understand a vector ~z which satisfies the property: z=~z P.(7) Suppose that the limit π= limn→∞ ~π(n)exists, then the probability vector ~π is called limiting distribution (limiting vector): π= lim n→∞ ~π(0)Pn.(8) In some literature there is written that, the stationary distribution is equal to the limiting distribution. Actually it is true only if the discrete time Markov chain (further as DTMC) satisfies the condition of ergodicity (irreducibility, aperiodicity). To be able to calculate the limiting distribution πwe will describe further properties of DTMC. A Markov chain jis periodic with period p, if on leaving state ja return is possible only in a number of transitions that is a multiple of the integer p > 1. For example Markov chain with following transition probability matrix: 0 1 1 0  is periodic with period p= 2. The DTMC is said to be irreducible if every state can be rached from every other state. The probability of ever returning to state jis denoted by fjj and is given by: fjj = ∞ X n=1 fn jj.(9) If fjj <1, then state jis said to be transient. If fjj = 1, the state is recurrent and we can define the mean recurrent time. As mean recurrent time Mjj we will understand: Mjj = ∞ X n=1 nfn jj.(10) If Mjj <∞(or equal to ∞), we say that the state is positive recurrent (or null recurrent). In a finite Markov chain, each state is either positive-recurrent or transient, and furthermore, at least one state must be positive-recurrent [6]. Theorem Let a Markov chain Cbe irreducible. Then Cis positive-recurrent or null recurrent or transient, i.e.: •all the states are positive-recurrent, or •all the states are null-recurrent, or •all the states are transient. If all the states of DTMC are positive-recurrent and aperiodic, then the probability distribution ~v(n)converges to a limiting distribution π, which is independent of initial distribution ~v(0). The πcan be obtained be solving the following equations: π=~π P, π > 0and X for all i πi= 1.(11) 3. Preliminaries In this section we will establish fundamental terminology of graph theory which will be used further in this paper. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 155 STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE 3.1. Definition of an Oriented Graph An oriented graph is an ordered pair G= (V;E), where Vis the set of vertices and Eis the set of edges. E⊆ V×V. 3.2. Definition of Network Network is a four-tupple S= (G;s;t;x), where: •Gis an oriented graph, •vertices s∈V(G),t∈V(G)are the source and sink, •x:E(G)→R+is a a positive labelling on edges, called edge capacities. 3.3. Definition of Flow in Network Flow in a network S= (G;s;t;x)is a function f: E(G)→R+ 0, where: •no edge capacity is exceeded: ∀e∈E(G):0≤≤ f(e)≤x(e), •the conservation of flow equation holds: ∀v∈V(G), v 6=z, s :Pe→vf(e) = Pe←vf(e). The value of a flow fis kfk=Pe←zf(e)−Pe→zf(e). 3.4. Definition of Cut Cut in the network S= (G;s;t;x)is such a set of edges C,C⊆E(G), such that in the factor of graph G,G−C, no oriented path remains. As minimal cut we will understand the cut where each proper subset of Cis not a cut in the network. 4. Problem Formulation The research presented here, was motivated by the practical problem. Analysed company was asked to quantify the probability that production fails was. Knowledge of risk, that the order won’t be delivered in time, is important for the partner’s firm to establish sufficient gods supplies. In the previous application see [1] a reliability analysis of the part of an industry process was realized. For each machine we could distinguish two modes – in order ’1’ or in fail ’0’. Thus the whole system could occur in one of the 2nstates where nis an amount of machines of analysed industry process. Since machines were organized in parallel it was easy to calculate a whole production of each state. The production of the certain state was calculated as a sum of production of functional states. More complicated situation occurs when the system is not connected in parallel. The aim of this paper is to present a calculation of a maximal production of each of 2nstates by usage a graph flow in network theory. In the following application we will simplify the process that, even at the beginning the gods can be reached from the begin to end of the process. We assume that for each state the process will produce a maximum possible production w. 5. Application In this section we will present an example how to use a flow in network theory in a reliability analysis. First, we will define states of the industry process and calculate a maximum production of each state by using the Ford Fulkerson algorithm. At the end we will demonstrate on the certain data how to calculate an expected value of production for a given time t. The analyzed industry process consists of six machines. For each machine we distinguish 2 different states one - work, zero - in fail. Thus the industrial process consists of 26different states. State denotes an ordered six -tuple of ones and zeroes. Let us describe the industry process of a firm with an oriented network. Every machine is represented as a vertex. The begin and the end of the processes represented by a source and sink of the network. The begin of a process consists of acceptance of gods and division between machines of the process. The end of a process consists of product inspection and expediting to the customers. Oriented edges describes the direction of the production process. Labelling of vertices represents the maximum amount of gods processed by the certain machine. To be able to work with labelled edges each vertex V1, V2, . . . V6from the structure of the process (Fig. 1) is replaced by two vertices connected by the edge labelled by the same value as former vertex. For our purposes all edges except newly created edges (with original labelling of vertices ) incident with source and sink are labelled by ∞. For each of 26state we will find the maximum flow in a network. For each state the edges incident with the vertices representing the machines occurred in failure will be removed from the network. To calculate the maximal production of each state we will find a max flow in network. In our application we will assume that the maximal flow is achieved for any possible state. Without this simplification, there would be nearly impossible to estimate the production of certain state. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 156 STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE Begin w(V1) = 60 w(V2) = 150 w(V3) = 150 w(V4) = 80 w(V5) = 80 w(V6) = 170 End Fig. 1: Structure of process. To find a maximal flow we will use a well known Ford-Fulkerson algorithm. For searching the graph edges a "Breadth-first" search was used. In the Tab. 1 there is a demonstration of several states and their production wi. To be able, to compute with DTMC, we need to calculate elements of transition probability matrix. To calculate the transition probability matrix Pwe need to calculate probabilities of failure pf, and probability of repair prfor each machine. To estimate probability pfthat the system fails during one hour we calculated the expected value as an average length of period fail Q: pf=1 Q.(12) Using the maximum likelihood method we estimated the probability prwhich says that a machine will be repaired within one hour: pr=∆V V,(13) where Vis an amount of all repairs that were realized and ∆Vis an amount of all repairs that lasted less than one hour. The calculated probabilities for given machines V1, V2,. . . V6are presented in the Tab. 2. Since there are 26different states the probability matrix will consist of 26+6 elements. With presented probabilities pf,pr we can calculate the probability transition matrix. For example p11 =Q6 i=1 pfi and p89 probability that the system changes from the state 8 to the state 9 in the Eq. 14, see the Tab. 1: p89 = (1 −pr1)·pr2·pf3·(1 −pf4)· (1 −pf5)·(1 −pf6)=0.64 ·0.36 ·0.005 · ·0.997 ·0.996 ·0.997 = 0.0011.(14) Tab. 1: Production of states. state V1V2V3V4V5V6prod. w 1 1 1 1 1 1 1 210 2 0 1 1 1 1 1 150 3 101111 60 4 110111 80 . . .. . .. . .. . .. . .. . .. . . 8 001111 0 9 010111 80 . . .. . .. . .. . .. . .. . .. . . 23 000111 0 . . .. . .. . .. . .. . .. . .. . . 64 000000 0 Each element of transition probability matrix is a result of multiplication of non-zero probabilities. Since the probability matrix consists of non-zero elements (some of them re close to zero), every state can be reached from every other in time t= 1 the analyzed DTMC is irreducible and aperiodic. Because of the DTMC is finite and irreducible it is also positiverecurrent. Then we can use the Eq. 6 to calculate the limiting distribution π. The Eq. 11 were solved numerically by using the "backslash"(implemented solver of linear equations) from Matlab program. In Tab. 4 there are presented few Elements of the 1st row of transition probability matrix and limiting distribution π. After calculating the transition probability matrix Pand the production of each state we can quantify the reliability of the industry process. To describe the reliability of the process we will calculate the expected value of production. The expected value of production Wwithin time Tis a sum of all expected values for each time step t∈ {1,2,· · · , T }. In our case, the expected value E(W)(t)of production for the certain time tis equal to: E(W)(t) = t X i=1 ~v(0)Pt~ W , (15) where ~ Wis a vector of productions wof all 26states. The results for calculating an expected value of production for time step 1 (time t) are presented in the Tab. 3 (maximal production for one time step is 170). According to the Tab. 3 we can conclude, that the expected value of production for the certain time quickly converges to the limit expected value (π·w). c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 157 STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE Tab. 2: Input probabilities. prob. V1V2V3V4V5V6 pf0.100 0.006 0.005 0.003 0.004 0.003 pr0.360 0.360 0.400 0.375 0.240 0.190 Tab. 3: Calculation of expected value. time expected value relative expected value 1 201.8095 0.9610 2 197.2069 0.9390 3 194.5720 0,9265 4 193.0340 0.9192 5 192.1160 0.9148 6 191.5550 0.9121 7 191.2030 0.9105 8 190.9750 0.9094 ∞190.4157 0.9067 Tab. 4: First row of transition probability matrix, limiting distribution π. p1,10.8810 π10.7300 p1,20.0980 π20.2030 p1,30.0050 π30.0120 p1,40.0040 π40.0090 p1,50.0030 π50.0050 p1,60.0040 π60.0120 p1,70.0010 π70.0110 p1,85.1905e−5π80.0030 . . .. . .. . .. . . p1,64 1.0800e−13 π64 8.8980e−11 6. Complexity of the Algorithm In this section we aim to estimate the complexity of the Ford-Fulkerson (F-F) algorithm for calculating the production ~ Wof all states. As a T(n)we will denote the worst-case time complexity of an algorithm with the input of size n. The complexity of F-F algorithm is T(n) = O(n3). In our case the size of input nis equal to the number of vertices |V|of the network. The algorithm should run up to 2|V|. Thus T(n)< O P|V| i=1 |V| i(i)3. The given formula is an upper bound of complexity of the algorithm for searching the vector of production ~ Win a network. To be able to estimate a complexity of the whole application we have to add a complexity of matrix exponentiation which is less than T(n)<2|V|3tand multiplication of first row of transition probability matrix by the vector ~ Wwhich is equal tu 2|V|. The last form is not very significant In comparison with complexity of matrix exponentiation and calculating the vector ~ Wof production. The complexity of calculating the production ~ Wwill also depend on the amount of vertices of the network and on the character of the network. Another way how to decrease the amount of computation is to use the minimal cut. Let Sffail state be a state with zero production wequal to zero. Let Cdenotes the set of vertices of the state Sf. Clearly the edges incident with Cform the cut of the network. Any other state where the vertices Care in failure has the zero production. For example in the Tab. 1 we can see that the production of the 8th (0,0,1,1,1,1) state is zero. Thus any state (for example the 23rd state), where the first and the second vertex is in fail has also the zero production. For those states, there is no need to use the FordFulkerson algorithm to calculate the maximal flow. For the large time tthere would be more suitable to to use limit probability distribution than to exponentiate the matrix. A problem could occur with saving the transition probability matrix in the computer memory. An amount of elements of the matrix grow exponentially (22n) with an amount of elements of analysed systems. The problem could by partially solved by rounding small values to zero and save the matrix as a sparse c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 158 STATISTICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE matrix. In this case there should be more complicated to proof the convergence of the DTMC. Another way how to work with large matrices is to dived the DTMC into more separate systems and calculate with their transition probability matrix separately. Other possibility how reduce the transition probability matrix is to shorter the time step and define only one possible change within the time step. It has the effect of simplifying the transition matrix. 7. Conclusion In the paper we have presented the application of graph flow in the network in the reliability analysis. The main worth of this work is an innovative usage of discrete time Markov chains and the flow in network theory in reliability analysis. The DTMC we used to describe the changes of the analysed process. For calculating for larger time step t, there is possible to use the limiting distribution instead of exponentiate the probability matrix. The Ford-Fulkerson algorithm was used to compute the industry production of analysed system. The main disadvantage of the presented approach is a high computational complexity. In further research there is possible to modify the Ford-Fulkerson algorithm to be more suitable for our purposes. Acknowledgment This work was supported by the Grant of SGS No. SP2013/116, VSB–Technical University of Ostrava, Czech Republic. References [1] SKALNY, P. and B. KRAJC. Discrete-Time Markov Chains in Reliability Analysis-Case Study. In: International Joint Conference CISIS’12ICEUTE’12-SOCO’12 Special Sessions. Ostrava: Springer Berlin Heidelberg, 2012, pp. 421–427. ISBN 978-3-642-33017-9. DOI: 10.1007/978-3-64233018-6_43. [2] WEST, D. B. Introduction to graph theory. N.J.: Prentice Hall, 2001. ISBN 01-301-4400-2. [3] PAGE, L. B. and J. E. PERRY. Reliability polynomials and link importance in networks. IEEE Transactions on Reliability. 1994, vol. 43, iss. 1, pp. 51–58. ISSN 0018–9529. DOI: 10.1109/24.285108. [4] DABROWSKI, Ch., F. HUNT and Katherine MORRISON. Improving the Efficiency of Markov Chain Analysis of Complex Distributed Systems. National Institute of Standards and Technology. Gaithersburg MD, 2010. Available at: http://www.nist.gov/customcf/get_pdf. cfm?pub_id=906149. [5] DABROWSKI, Ch. and F. HUNT. Using Markov Chain and Graph Theory Concepts to Analyze Behavior in Complex Distributed Systems. In: The 23rd European Modeling and Simulation Symposium. Rome: Curran Associates, Inc., 2011. Available at: http://www.nist.gov/itl/antd/upload/ DabrowskiHunt-Paper129-updated.pdf. [6] Stewart, W. Introduction to the Numerical Solution of Markov Chains. Princeton University Press, 1994. ISBN 0-691-03699-3. [7] CORMEN, T. H. Introduction to algorithms. Cambridge: MIT Press, 2001. ISBN 02-620-32937. About Authors Pavel SKALNY was born in Ostrava. He received his M.Sc. from Mathematics and History in Silesian University in Opava in 2009. His research interests include Stochastic modelling in reliability and fracture mechanics modelling. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 159