# A Hidden Semi-Markov Model with Duration-Dependent State Transition Probabilities for Prognostics

Authors: Ning Wang; Shu-dong Sun; Zhi-qiang Cai; Shuai Zhang; Can Saygin
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Realistic prognostic tools are essential for effective condition-based maintenance systems. In this paper, a Duration-Dependent Hidden Semi-Markov Model (DD-HSMM) is proposed, which overcomes the shortcomings of traditional Hidden Markov Models (HMM), including the Hidden Semi-Markov Model (HSMM): (1) it allows explicit modeling of state transition probabilities between the states; (2) it relaxes observations’ independence assumption by accommodating a connection between consecutive observations; and (3) it does not follow the unrealistic Markov chain’s memoryless assumption and therefore it provides a more powerful modeling and analysis capability for real world problems. To facilitate the computation of the proposed DD-HSMM methodology, new forward-backward algorithm is developed. The demonstration and evaluation of the proposed methodology is carried out through a case study. The experimental results show that the DD-HSMM methodology is effective for equipment health monitoring and management.


---

## Body

## 1. Introduction

Fault is a change from the normal operating condition of a system to an abnormal condition, which occurs as a result of system performance degradation over time [1]. Diagnostics indicates the occurrence of a fault and its root cause. Prognostics is fault prediction method; it involves detection of a pending fault before it occurs, identifying its root cause and estimating the remaining useful life (RUL), which is also known as time-to-failure [2]. Condition-based Maintenance (CBM) is a maintenance program that recommends maintenance actions based on the information collected through condition monitoring.A CBM program can be used to do diagnostics or prognostics; however, regardless of the application, it follows three steps [3–5]. First, data relevant to events and system health are collected through data acquisition techniques. Data acquisition in CBM includes event-type data (i.e., information of what happened) and condition monitoring data, which are the measurements related to system health. Second, event and condition monitoring data are interpreted for better understanding in the data processing step. Finally, maintenance decisions are made based on the interpretation and analysis of data. In particular, to identify the weakest components and states and improve the efficiency of CBM, the integrated importance measure of multistate system was introduced by Si et al. [6, 7]. An extensive survey on machine diagnostics and prognostics implementing condition-based monitoring can be found in Jardine et al. [8] and Heng et al. [9].Data analysis for event data only is reliability analysis, which maps the event data over a time axis to determine the probability of events and uses the probability distribution to predict failures. On the other hand, data acquisition in CBM provides event and condition monitoring data. Therefore, it is more effective to combine events and conditions in a model in order to do diagnostics or prognostics. Hidden Markov model (HMM) is a technique for modeling and analyzing event and condition monitoring data together. It consists of two stochastic processes:(1) a Markov chain with finite number of states describing an underlying mechanism and (2) an observation process depending on the hidden state [10–12]. An HMM contains finite states connected by transitions. Each state is characterized by a transition probability and an observation probability [13].Researchers have proposed a number of techniques to address these limitations. Continuous variable duration HMM is adopted in the speech recognition. Compared to standard HMM, results show that the absence of a correct duration model increases the error rate by 50% [14–16]. Another example is in handwritten word recognition area; due to the inherent ambiguity related to the segmentation process in handwritten words, it is a practical idea to use the variable duration model for the states in a HMM-based handwritten word recognition system [17, 18]. Recently, some researchers apply HMM in the area of diagnostics and prognostics in machining process [19, 20]. However, these studies use only ordinary HMM technique. The inherent limitation of HMM as mentioned above still exists in these models.Prognostic methods used in CBM are often a combination of statistical inference and machine learning methods [4, 21–24]. Model-based methods assume that measured information is stochastically correlated with the actual machine condition. HMM identifies the actual machine conditions from observable monitored data through a statistical approach. HMM has been very effective in various applications ranging from speech recognition [10, 14, 25–27] to tool wear monitoring and machining [3, 20, 28, 29].The primary advantage of HMM is its robust mathematical foundation that can allow for many practical applications and different areas of use. An added benefit of employing HMMs is the ease of model interpretation in comparison with pure “black-box” modeling methods such as artificial neural networks that are often employed in advanced diagnostic models [28]. However, an inherent limitation of HMM approach is that its state duration follows an exponential distribution. In other words, HMM does not provide adequate representation of temporal structure.To overcome the limitations of HMM in prognosis, Dong and He [30, 31] propose a Hidden Semi-Markov Model-based (HSMM) methodology by adding an explicit temporary structure into HMM to predict RUL of equipment. In this model, the states of HSMM are used to represent the health status of equipment. The trained HSMM can be used to diagnose the health state of equipment. Through parameter estimation of the health-state duration probability distribution and the proposed backward recursive equations, the RUL of the equipment can be predicted [32]. Although the results from HSMM are promising, the deterioration in the same state of the system is not taken into consideration in this model. It assumed that the state transition probabilities stay the same in the same state, which assumes all observations are independent, which typically does not hold in real world applications.This paper presents a new approach that expands the HSMM methodology [30, 31, 33] with duration-dependent state transition probabilities. Different from HSMM, the proposed Duration-Dependent Hidden Semi-Markov Model (DD-HSMM) does not follow the unrealistic Markov chain’s memoryless assumption and therefore provides a more realistic and powerful modeling and analysis capability. The major contribution of the DD-HSMM methodology is that it allows explicit modeling of the transition probabilities, which (1) do not only depend on the state but also (2) vary with the duration of each state, and (3) it provides the capability to relax observations’ independence assumption by accommodating a link between consecutive observations, which makes it more realistic in real world applications.


## 2. Theoretical Background

### 2.1. Description of General HSMM

A Hidden Semi-Markov Model (HSMM) is an extension of HMM by allowing the underlying process to be a semi-Markov chain with a variable duration or sojourn time for each state. The HSMM model is an ideal mathematical model for estimating the unobservable health states with observable sensor signal. For example, a small change in a bearing alignment could cause a small nick in the bearing, which could cause scratches in the bearing race and additional nicks, leading to complete bearing failure. This process can be well described by the HSMM. Letst be the hidden state at time t and O the observation sequence; a HSMM is characterized by its parameters. The parameters of a HSMM are as follows: the initial state distribution (denoted by π), the transition model (denoted by A), the observation matrix (denoted by C), and the state duration distribution (denoted by D). Thus, a HSMM can be written as λ=(π,A,C,D). For a given state S, C is the probability matrix of observation being ot at time t and ot+1 at time t+1.


### 2.2. State Transitions

In HSMM, there areN states, and the transitions between the states are according to the transition matrix A; that is, P(i→j)=aij. Similar to standard HMM, the state s0 at time t=0 is a special state “START.”Although the distinct health-state transitionsql-1→sql is Markov,
(1)P[sql=j∣sql-1=i]=aij,
the state transition st-1→st is usually not Markov. It is the reason why the model is called “semi-Markov” [32], which means in the HSMM case, the conditional independence between the past and the future is only ensured when the process moves from one health state to another distinct health state.


### 2.3. Inference Procedures

Similar to HMM, HSMM also has basic problems to deal with, that is, evaluation, recognition, and training problems:(1)
evaluation (also called classification): given the observation sequenceO=o1o2⋯oT and a HSMM λ, what is the probability of the observation sequence given the model, that is, P(O∣λ);(2)
decoding (also called recognition): given the observation sequenceO=o1o2⋯oT and a HSMM λ, what sequence of hidden states S=s1s2⋯sT most probably generates the given sequence of observations;(3)
learning (also called training): how do we adjust the model parametersλ=(π,A,D,B) to maximize P(O∣λ).Different algorithms have been developed for above three problems. The most straightforward way of solving the evaluation problem is enumerating every possible state sequence of lengthT (the number of observations). However, the computation burden for this exhaustive enumeration is prohibitively high. Fortunately, there is a more efficient algorithm that is based on dynamic programming, called forward-backward procedure. The goal for decoding problem is to find the optimal state sequence associated with the given observation sequence. The most widely used optimality criterion is to find the single best state sequence (path), that is, to maximize P(S∣O,λ) that is equivalent to maximizing P(S,O∣λ). Viterbi algorithm is used to find this single best state sequence, which is based on dynamic programming methods. For learning problem, there is no known way to obtain analytical solution. However, we can adjust the model parameter λ=(π,A,D,B) such that P(O∣λ) is locally maximized using an iterative procedure, such as the Baum-Welch method (or equivalently the Expectation-Maximization algorithm).


## 2.1. Description of General HSMM

A Hidden Semi-Markov Model (HSMM) is an extension of HMM by allowing the underlying process to be a semi-Markov chain with a variable duration or sojourn time for each state. The HSMM model is an ideal mathematical model for estimating the unobservable health states with observable sensor signal. For example, a small change in a bearing alignment could cause a small nick in the bearing, which could cause scratches in the bearing race and additional nicks, leading to complete bearing failure. This process can be well described by the HSMM. Letst be the hidden state at time t and O the observation sequence; a HSMM is characterized by its parameters. The parameters of a HSMM are as follows: the initial state distribution (denoted by π), the transition model (denoted by A), the observation matrix (denoted by C), and the state duration distribution (denoted by D). Thus, a HSMM can be written as λ=(π,A,C,D). For a given state S, C is the probability matrix of observation being ot at time t and ot+1 at time t+1.


## 2.2. State Transitions

In HSMM, there areN states, and the transitions between the states are according to the transition matrix A; that is, P(i→j)=aij. Similar to standard HMM, the state s0 at time t=0 is a special state “START.”Although the distinct health-state transitionsql-1→sql is Markov,
(1)P[sql=j∣sql-1=i]=aij,
the state transition st-1→st is usually not Markov. It is the reason why the model is called “semi-Markov” [32], which means in the HSMM case, the conditional independence between the past and the future is only ensured when the process moves from one health state to another distinct health state.


## 2.3. Inference Procedures

Similar to HMM, HSMM also has basic problems to deal with, that is, evaluation, recognition, and training problems:(1)
evaluation (also called classification): given the observation sequenceO=o1o2⋯oT and a HSMM λ, what is the probability of the observation sequence given the model, that is, P(O∣λ);(2)
decoding (also called recognition): given the observation sequenceO=o1o2⋯oT and a HSMM λ, what sequence of hidden states S=s1s2⋯sT most probably generates the given sequence of observations;(3)
learning (also called training): how do we adjust the model parametersλ=(π,A,D,B) to maximize P(O∣λ).Different algorithms have been developed for above three problems. The most straightforward way of solving the evaluation problem is enumerating every possible state sequence of lengthT (the number of observations). However, the computation burden for this exhaustive enumeration is prohibitively high. Fortunately, there is a more efficient algorithm that is based on dynamic programming, called forward-backward procedure. The goal for decoding problem is to find the optimal state sequence associated with the given observation sequence. The most widely used optimality criterion is to find the single best state sequence (path), that is, to maximize P(S∣O,λ) that is equivalent to maximizing P(S,O∣λ). Viterbi algorithm is used to find this single best state sequence, which is based on dynamic programming methods. For learning problem, there is no known way to obtain analytical solution. However, we can adjust the model parameter λ=(π,A,D,B) such that P(O∣λ) is locally maximized using an iterative procedure, such as the Baum-Welch method (or equivalently the Expectation-Maximization algorithm).


## 3. Inference and Learning Mechanisms of DD-HSMM

### 3.1. Model Structure

Although HSMM has explicit state duration probability distributionPi(d), the sate transition probabilities aij are duration invariant. In this paper, we replace duration-invariant state transition probabilities with duration-dependent state transition probabilities. The parameters for a DD-HSMM are as follows: the initial state distribution denoted by π={πi}1≤i≤N(πi=P[s1=i]1≤i≤N); the transition model denoted by A={aij(d)}(1≤i,j≤N,1≤d≤D); the observation matrix denoted by B={bi(k)}  (bi(k)=P[vk∣st=i], where 1≤i≤N,1≤k≤M,M is the observation number in i, and vi={v1,v2,…,vM} is the observation symbols in state i); and the state duration distribution denoted by D(Pi(d)(i=1,2,…,N)). Thus, a DD-HSMM can be written as λ′=(π,A,D,B). Here, for the duration in given state S is dt(i)=d, the state transition probability aij(d)=P(qt+1=j∣qt=i,dt(i)=d),1≤i,j≤N,1≤d≤Di (N is the state number, Di is the max staying time in state i). And the state transition probabilities satisfy the constraint: ∑j=1Naij(d)=1,1≤d≤Di,1≤i≤N.


### 3.2. Duration-Dependent State Transition Probability

In DD-HSMM, the state transition probability distributionA={aij(d)},1≤i,j≤N,1≤d≤D. We define duration-dependents state transition probabilities as follows:
(2)aij(d)=P(qt+1=jqt=i,dt(i)=d),
where N and D are the number of states and the maximum duration in any states, respectively. Equation (2) represents the transition from state i to state j, given that the duration in state i at time t is dt(i)=d. It indicated that, in the DD-HSMM case, the state transition probability is not only state dependent, but also duration variant.


### 3.3. Inference Procedures

Similar to HSMM, DD-HSMM also has basic problems to deal with, that is, evaluation, recognition, and training problems. To facilitate the computation in the proposed DD-HSMM-based health prediction model, in the following, new forward-backward variables are defined and modified forward-backward algorithm is developed.A dynamic programming scheme is employed for the efficient computation of the inference procedures. To implement the inference procedures, a forward variableαt(i,d) is defined as the probability of generating o1,o2,…,ot and ending in state i and the duration dt(i)=d:
(3)αt(i,d)=P(o1,o2,…,ot,qt=i,dt(i)=d∣λ).
The initial conditions are established at time t=1 as follows:
(4)α1(1,d1)=1,α1(j,1)=α1(1,d1)a1j(d1)bj(o1),α1(N,dN)=∑i=2Nα1(i,1)Pi(1)aiN(1).
All unspecified α values are zero. For time t=2,…,T,
(5)αt(j,1)=∑i=2i≠jN-1∑d=1Diαt-1(i,d)Pi(d)aij(d)bj(ot),αt(j,d)=αt-1(j,d-1)bj(ot),αt(N,dN)=∑i=2N-1∑d=1Diαt(i,d)Pi(d)aiN(d),
where aij(d) is the state transition probability from state i to state j, given that the duration in state i at time t is dt(i)=d. bj(ot) is the output probability of observation vector ot from state j and pi(d) is the state duration probability of state i. N is the number of states in DD-HSMM and Di is the maximum duration in state i.Similar to the forward variable, the backward variable can be written as(6)βt(i,d)=P(ot,…,oTqt=i,dt(i)=d,λ).For the backward probability, the initial conditions are set at timet=T as follows:
(7)βT(N,dN)=1,βT(i,d)=Pi(d)aiN(d),βT(1,d1)=∑j=2N-1α1j(d1)bj(oT)βT(j,1).For time0<t<T,
(8)βt(i,d)=Pi(d)∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+bj(ot+1)βt+1(i,d+1)1<i<N,(9)βt(1,d1)=∑j=2N-1a1j(d1)bj(ot)βt(j,1).Then the total probability can be computed by(10)Pr=P(O∣λ)=∑i=1N∑d=1Diαt(i,d)βt(i,d).


### 3.4. Modified Forward-Backward Algorithm for DD-HSMM

In order to give reestimation formulas for all variable of the DD-HSMM, one DD-HSMM-featured forward-backward variable is defined:(11)ξt(i,j,d)=p(t=qn,st=i,t′=qn+1,st′=j,dt(i)=pki=d∣O,λ′).
In this equation, ξt(i,j,d) is the probability of state transition from state i to state j at time t+1 after being in state i for a duration of dt(i)=d, given the model λ and observation O. From the definition of the forward-backward variables, we can derive ξt(i,j,d) as follows:
(12)ξt,t′(i,j,d)=P(t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d∣O0T,λ)=P(t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d,O0T∣λ)×(P(O0T∣λ))-1=P(Ot+1t′t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d,O0t,Ot+1t′,Ot′+1T∣λ)×(P(O0T∣λ))-1=P(t=qn,st=i,dt(i)=d,O0t∣λ)×P(t′=qn+1,st′=j,dt(i)=d,Ot+1t′,=PKKOt′+1T∣t=qn,st=i,dt(i)=d,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j,Ot+1t′,Ot′+1T∣t=αt(i,d)Pkii=qn,st=i,dt(i)=d,O0t,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j,Ot+1t′,Ot′+1T∣t=αt(i,d)Pkii=qn,st=i,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j∣t=qn,=αt(i,d)Pkiist=i,dt(i)=d,λ)×P(Ot+1t′,Ot′+1T∣t=qn,st=i,=αt(d)KKt′=qn+1,st′=j,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)aij(d)P(Ot+1t′,Ot′+1T∣t=qn,st=i,=αt(i,d)aij(d)PiKKt′=qn+1,st′=j,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1={αt(1,d1)a1j(d1)bj(ot)βt(j,1)P(O0T∣λ),i=1αt(i,d)pi(d)aij(d)×bj(ot+1)βt+1(j,1)×(P(O0T∣λ))-1,2≤i,j<Nαt(i,d)pi(d)×aiN(d)βt(i,d)×(P(O0T∣λ))-1,2≤i<N,j=N.
Then, we have
(13)ξt(i,j,d)=p(t=qn,st=i,t′=qn+1,st′=pki=j,dt(i)=d∣O,λ′)={1Prαt(1,d1)a1j(d1)bj(ot)βt(j,1),i=11Prαt(i,d)pi(d)aij(d)×bj(ot+1)βt+1(j,1),2≤i,j<N1Prαt(i,d)pi(d)2≤i<N,j=N.×aiN(d)βt(i,d),
and the probability in state i at time t with duration of d is defined as γt(i,d), and, from the definition of the forward-backward variables, we can easily derive γt(i,d) as follows:
(14)γt(i,d)=p(qt=i,di=dO,λ)=1Prαt(i,d)βt(i,d).The forward-Backward algorithm computes the following probabilities.Forward Pass. The forward pass of the algorithm computes αt(i,d).Step 1 (initialization (t=1)).
The forward variable is shown as follows:(15)α1(i,d)={1j=1,d=d1α1(1,d1)a1j(d1)bj(o1)1<j<N,d=1∑i=2Nα1(i,1)Pi(1)aiN(1)j=N,d=dN.Step 2 (forward recursion (t>1)).
Fort=2,…,T,
(16)αt(j,d)={∑i=2i≠jN-1∑d=1Diαt-1(i,d)Pi(d)aij(d)bj(ot)d=1αt-1(j,d-1)bj(ot)∑i=2N-1∑d=1Diαt(i,d)Pi(d)aiN(d)j=N.Backward Pass. The backward pass computes βt(i,d)Step 1 (initialization (t=T)).
The backward variable is shown as follows:(17)βT(i,d)={1i=N,d=dNPi(d)aiN(d)1<i<N,d1<d<dN∑j=2N-1α1j(d1)bj(oT)βT(j,1)i=1,d=d1.Step 2 (backward recursion (t<T)).
Fort=2,…,T,
(18)βt(i,d)={Pi(d)∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+bj(ot+1)βt+1(i,d+1)1<i<N∑j=2N-1a1j(d1)bj(ot)βt(j,1)i=1,d=di.


### 3.5. Parameter Reestimation for DD-HSMM

The reestimation formula for initial state distribution is the probability that stateiwas the first state, given O:
(19)π¯i=πi[∑d=1Dβ(i,d)P(d∣i)(co0o1ico0o1i⋯cod-1odi)]P(O∣λ).
The reestimation formula of state transition probabilities is the ratio of expected number of transition from state i to state j, to the expected number of transitions from state i:
(20)aij¯(d)=∑t=1Tξt(i,j,d)∑t=1Tγt(i,d).


### 3.6. Training of State Duration Models Using Parametric Probability Distributions

In this paper, state duration densities are modeled by single Gaussian distribution estimated from training data. The existing state duration estimation method is through the simultaneous training DD-HSMM and their duration densities. However, these techniques are inefficient because the training process requires huge storage and computational load. Therefore, a new approach is adopted for training state duration models. In this approach, state duration probabilities are estimated on the lattice (or trellis) of observations and states, which is obtained in the DD-HSMM training stage.The meanμi and variance σ(i) of duration probability of state i are determined by
(21)μ(i)=∑t=1T∑d=1Diχt(i,d)d∑t=1T∑d=1Diχt(i,d).σ-(i)=∑t=1T∑d=1Diχt(i,d)d2∑t=1T∑d=1Diχt(i,d)-[μ-(i)]2,
In these equations, χt(i,d) is the probability of state i at time t with the duration of dt(i)=d and χt(i,d) can present as
(22)χt(i,d)=1prαt(i,d)pi(d)[∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+aiN(d)βt(N,dN)∑j=2j≠iN-1aij(d)].


## 3.1. Model Structure

Although HSMM has explicit state duration probability distributionPi(d), the sate transition probabilities aij are duration invariant. In this paper, we replace duration-invariant state transition probabilities with duration-dependent state transition probabilities. The parameters for a DD-HSMM are as follows: the initial state distribution denoted by π={πi}1≤i≤N(πi=P[s1=i]1≤i≤N); the transition model denoted by A={aij(d)}(1≤i,j≤N,1≤d≤D); the observation matrix denoted by B={bi(k)}  (bi(k)=P[vk∣st=i], where 1≤i≤N,1≤k≤M,M is the observation number in i, and vi={v1,v2,…,vM} is the observation symbols in state i); and the state duration distribution denoted by D(Pi(d)(i=1,2,…,N)). Thus, a DD-HSMM can be written as λ′=(π,A,D,B). Here, for the duration in given state S is dt(i)=d, the state transition probability aij(d)=P(qt+1=j∣qt=i,dt(i)=d),1≤i,j≤N,1≤d≤Di (N is the state number, Di is the max staying time in state i). And the state transition probabilities satisfy the constraint: ∑j=1Naij(d)=1,1≤d≤Di,1≤i≤N.


## 3.2. Duration-Dependent State Transition Probability

In DD-HSMM, the state transition probability distributionA={aij(d)},1≤i,j≤N,1≤d≤D. We define duration-dependents state transition probabilities as follows:
(2)aij(d)=P(qt+1=jqt=i,dt(i)=d),
where N and D are the number of states and the maximum duration in any states, respectively. Equation (2) represents the transition from state i to state j, given that the duration in state i at time t is dt(i)=d. It indicated that, in the DD-HSMM case, the state transition probability is not only state dependent, but also duration variant.


## 3.3. Inference Procedures

Similar to HSMM, DD-HSMM also has basic problems to deal with, that is, evaluation, recognition, and training problems. To facilitate the computation in the proposed DD-HSMM-based health prediction model, in the following, new forward-backward variables are defined and modified forward-backward algorithm is developed.A dynamic programming scheme is employed for the efficient computation of the inference procedures. To implement the inference procedures, a forward variableαt(i,d) is defined as the probability of generating o1,o2,…,ot and ending in state i and the duration dt(i)=d:
(3)αt(i,d)=P(o1,o2,…,ot,qt=i,dt(i)=d∣λ).
The initial conditions are established at time t=1 as follows:
(4)α1(1,d1)=1,α1(j,1)=α1(1,d1)a1j(d1)bj(o1),α1(N,dN)=∑i=2Nα1(i,1)Pi(1)aiN(1).
All unspecified α values are zero. For time t=2,…,T,
(5)αt(j,1)=∑i=2i≠jN-1∑d=1Diαt-1(i,d)Pi(d)aij(d)bj(ot),αt(j,d)=αt-1(j,d-1)bj(ot),αt(N,dN)=∑i=2N-1∑d=1Diαt(i,d)Pi(d)aiN(d),
where aij(d) is the state transition probability from state i to state j, given that the duration in state i at time t is dt(i)=d. bj(ot) is the output probability of observation vector ot from state j and pi(d) is the state duration probability of state i. N is the number of states in DD-HSMM and Di is the maximum duration in state i.Similar to the forward variable, the backward variable can be written as(6)βt(i,d)=P(ot,…,oTqt=i,dt(i)=d,λ).For the backward probability, the initial conditions are set at timet=T as follows:
(7)βT(N,dN)=1,βT(i,d)=Pi(d)aiN(d),βT(1,d1)=∑j=2N-1α1j(d1)bj(oT)βT(j,1).For time0<t<T,
(8)βt(i,d)=Pi(d)∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+bj(ot+1)βt+1(i,d+1)1<i<N,(9)βt(1,d1)=∑j=2N-1a1j(d1)bj(ot)βt(j,1).Then the total probability can be computed by(10)Pr=P(O∣λ)=∑i=1N∑d=1Diαt(i,d)βt(i,d).


## 3.4. Modified Forward-Backward Algorithm for DD-HSMM

In order to give reestimation formulas for all variable of the DD-HSMM, one DD-HSMM-featured forward-backward variable is defined:(11)ξt(i,j,d)=p(t=qn,st=i,t′=qn+1,st′=j,dt(i)=pki=d∣O,λ′).
In this equation, ξt(i,j,d) is the probability of state transition from state i to state j at time t+1 after being in state i for a duration of dt(i)=d, given the model λ and observation O. From the definition of the forward-backward variables, we can derive ξt(i,j,d) as follows:
(12)ξt,t′(i,j,d)=P(t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d∣O0T,λ)=P(t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d,O0T∣λ)×(P(O0T∣λ))-1=P(Ot+1t′t=qn,st=i,t′=qn+1,st′=j,=Pkidt(i)=d,O0t,Ot+1t′,Ot′+1T∣λ)×(P(O0T∣λ))-1=P(t=qn,st=i,dt(i)=d,O0t∣λ)×P(t′=qn+1,st′=j,dt(i)=d,Ot+1t′,=PKKOt′+1T∣t=qn,st=i,dt(i)=d,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j,Ot+1t′,Ot′+1T∣t=αt(i,d)Pkii=qn,st=i,dt(i)=d,O0t,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j,Ot+1t′,Ot′+1T∣t=αt(i,d)Pkii=qn,st=i,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)P(t′=qn+1,st′=j∣t=qn,=αt(i,d)Pkiist=i,dt(i)=d,λ)×P(Ot+1t′,Ot′+1T∣t=qn,st=i,=αt(d)KKt′=qn+1,st′=j,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1=αt(i,d)aij(d)P(Ot+1t′,Ot′+1T∣t=qn,st=i,=αt(i,d)aij(d)PiKKt′=qn+1,st′=j,dt(i)=d,ot,λOt+1t′)×(P(O0T∣λ))-1={αt(1,d1)a1j(d1)bj(ot)βt(j,1)P(O0T∣λ),i=1αt(i,d)pi(d)aij(d)×bj(ot+1)βt+1(j,1)×(P(O0T∣λ))-1,2≤i,j<Nαt(i,d)pi(d)×aiN(d)βt(i,d)×(P(O0T∣λ))-1,2≤i<N,j=N.
Then, we have
(13)ξt(i,j,d)=p(t=qn,st=i,t′=qn+1,st′=pki=j,dt(i)=d∣O,λ′)={1Prαt(1,d1)a1j(d1)bj(ot)βt(j,1),i=11Prαt(i,d)pi(d)aij(d)×bj(ot+1)βt+1(j,1),2≤i,j<N1Prαt(i,d)pi(d)2≤i<N,j=N.×aiN(d)βt(i,d),
and the probability in state i at time t with duration of d is defined as γt(i,d), and, from the definition of the forward-backward variables, we can easily derive γt(i,d) as follows:
(14)γt(i,d)=p(qt=i,di=dO,λ)=1Prαt(i,d)βt(i,d).The forward-Backward algorithm computes the following probabilities.Forward Pass. The forward pass of the algorithm computes αt(i,d).Step 1 (initialization (t=1)).
The forward variable is shown as follows:(15)α1(i,d)={1j=1,d=d1α1(1,d1)a1j(d1)bj(o1)1<j<N,d=1∑i=2Nα1(i,1)Pi(1)aiN(1)j=N,d=dN.Step 2 (forward recursion (t>1)).
Fort=2,…,T,
(16)αt(j,d)={∑i=2i≠jN-1∑d=1Diαt-1(i,d)Pi(d)aij(d)bj(ot)d=1αt-1(j,d-1)bj(ot)∑i=2N-1∑d=1Diαt(i,d)Pi(d)aiN(d)j=N.Backward Pass. The backward pass computes βt(i,d)Step 1 (initialization (t=T)).
The backward variable is shown as follows:(17)βT(i,d)={1i=N,d=dNPi(d)aiN(d)1<i<N,d1<d<dN∑j=2N-1α1j(d1)bj(oT)βT(j,1)i=1,d=d1.Step 2 (backward recursion (t<T)).
Fort=2,…,T,
(18)βt(i,d)={Pi(d)∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+bj(ot+1)βt+1(i,d+1)1<i<N∑j=2N-1a1j(d1)bj(ot)βt(j,1)i=1,d=di.


## 3.5. Parameter Reestimation for DD-HSMM

The reestimation formula for initial state distribution is the probability that stateiwas the first state, given O:
(19)π¯i=πi[∑d=1Dβ(i,d)P(d∣i)(co0o1ico0o1i⋯cod-1odi)]P(O∣λ).
The reestimation formula of state transition probabilities is the ratio of expected number of transition from state i to state j, to the expected number of transitions from state i:
(20)aij¯(d)=∑t=1Tξt(i,j,d)∑t=1Tγt(i,d).


## 3.6. Training of State Duration Models Using Parametric Probability Distributions

In this paper, state duration densities are modeled by single Gaussian distribution estimated from training data. The existing state duration estimation method is through the simultaneous training DD-HSMM and their duration densities. However, these techniques are inefficient because the training process requires huge storage and computational load. Therefore, a new approach is adopted for training state duration models. In this approach, state duration probabilities are estimated on the lattice (or trellis) of observations and states, which is obtained in the DD-HSMM training stage.The meanμi and variance σ(i) of duration probability of state i are determined by
(21)μ(i)=∑t=1T∑d=1Diχt(i,d)d∑t=1T∑d=1Diχt(i,d).σ-(i)=∑t=1T∑d=1Diχt(i,d)d2∑t=1T∑d=1Diχt(i,d)-[μ-(i)]2,
In these equations, χt(i,d) is the probability of state i at time t with the duration of dt(i)=d and χt(i,d) can present as
(22)χt(i,d)=1prαt(i,d)pi(d)[∑j=2j≠iN-1aij(d)bj(ot+1)βt+1(j,1)+aiN(d)βt(N,dN)∑j=2j≠iN-1aij(d)].


## 4. DD-HSMM Based Health Prognostic

Many applications in the actuarial, econometric, engineering, and medical literature involve the use of the hazard rate function [33]. The mathematical properties of HR function can reveal a variety of features in the data.LetT denote the time to failure of an item under consideration, with lifetime distribution function F(t) and reliability function R(t), where F(t)+R(t)=1 and F(0)=0. Assume that F(0)=0 and density function f(t)=F′(t) exist; then the HR function can be defined as:
(23)λ(t)=lim⁡N→∞Δt→0Δm(t)[N-m(t)]Δt=dm(t)[N-m(t)]dt=dm(t)/M1-F(t)=f(t)R(t).
In which, M is the total number of sample items, m(t) is the number of items that fail before time t, and Δm(t) is the number of items that fail during the time interval (t,t+Δt). The ERL function μ(t) is the expected time remaining to failure, given that the system has survived to time t; then μ(t)=E(T-t∣T>t)=(1/R(t))∫t∞R(x)dx, for t such that R(t)>0. Therefore, λ(t) can be approximated as the conditional probability of failure during the time interval (t,t+Δt) given survival to time t.Suppose that a machine will go through health statesi,(i=1,2,…,N-1) before entering failure state N. Let D(i) denote the expected duration of the machine staying at health state i; based on the parameters estimated above, we can get D(i) as follows:
(24)D(i)=μ(i)+ρδ2(i).
And ρ can be denoted by
(25)ρ=(T-∑l=0L-1μ(i))∑l=0L-1δ2(i).
Then, once the machine has entered the health state i, its expected residual life equals the summation of the expected residual duration of the machine staying at health state i and the total remaining staying in the future health states before failure. Denote D^(id) as the expected residual duration of the machine staying in the health state i for d. When the equipment entered state i at time ti, the conditional probability of failure during (ti+d,ti+(d+Δt)) can be defined as the probability that the machine will transit to any other state during the coming Δt and the probability that the machine still stay at state i. It can be seen from (9) and (10) that λ^(t+d)Δt can be denoted as follows:
(26)λ^(t+d)Δt=ξt(i,j,d)γt(i,d).
Then
(27)D^(id)=D(i)(1-λ^(t+d)Δt)=D(i)(1-ξt(i,j,d)γt(i,d)).
The DD-HSMM equipment health prediction procedure is given as follows.Step  1. From the DD-HSMM training procedure (i.e., parameter estimation), the state transition probability for the DD-HSMM can be obtained.Step   2. Through the DD-HSMM parameter estimation, the duration probability density function for each health-state can be obtained. Therefore, the duration mean and variance can be calculated.Step   3. By classification, identify the current health status of the equipment.Step   4. The remaining useful life (RUL) of equipment can be predicted by the following formula (suppose that the equipment currently stays at health state i with duration of d):
(28)RULi(d)=D^(id)+∑j=i+1N-1D(j).


## 5. Case Study

In this case study, long-term wear experiments on rolling element bearings were conducted [1]. In order to collect adequate amount of data sets for the validation of the proposed scheme, three experiments with normal operating conditions, three experiments with cage defect fault, and three experiments each of inner and outer race defect faults were performed until the bearing reached a complete failure state and stopped operating. Bearing characteristic frequencies in the frequency domain are extracted from the vibration signals corresponding to different degrees of the health states of the bearing acquired during experiments.During the test running, under each condition, vibration signals were collected. These signals were extracted using a Mahalanobis-Taguchi System (MTS) based model in the original paper [1] and used for the proposed DD-HSMM methodology in this paper. The expert judgment is made of four integer numbers ranging from 0 to 3, representing 4 system states, as follows:0→
the bearing is operating normally;1→
the bearing is operating and shows signs of deterioration; it is advisable to take some preventive action at the next planned maintenance;2→
the bearing is operating but requires immediate attention;3→
the bearing has failed.


### 5.1. Operation State Identification

In order to identify the accuracy of the operation state identification method proposed in this paper, experimental data with normal operating condition were obtained. The experimental data set included 50 samples for each state (denoted by 0, 1, 2, and 3). Of these data points, 20 of them were used to train the model, and the remaining 30 samples were used to validate the model.In the DD-HSMM, mixture Gaussian distribution and the single Gaussian distribution were used to model the output probability distribution and the state duration densities separately, in which the number of states is 4. The maximum number of iterations in training process is set to 100 and the convergence error to 0.000001.The DD-HSMM-based training model is shown as Figure1. Thex-axis shows the training steps and they-axis represents the likelihood probability of different states. As can be seen from Figure 1, the progression of the four states reaches the set error in less than 40 steps. This demonstrates the potential of the model to have a strong real-time signal processing capability.Figure 1
Training curve of the DD-HSMM model.The classification results obtained on the remaining 30 data samples are shown in Table1. As indicated in the results, the accuracy of the DD-HSMM method is 94.2%.Table 1
Prognostics results.
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Recognition accuracy
Normal 0
29
1
0
0
96.7%
Degradation 1
2
27
1
0
90%
Degradation 2
0
1
28
1
93.3%
Failure  3
0
0
1
29
96.7%
Total accuracy
94.2%


### 5.2. Health Prediction for RUL

As described before, a four-state DD-HSMM prediction model is constructed. In the training process, even if the device is in the same running condition, the dwell time is different, transition probabilities between states and the mean or variance of duration in each state are not the same. Tables2 and 3 show the state transition probability, the mean, and variance of duration in each state when dt(1)=1, representing the bearing in state 1 with duration of 1. Tables 4 and 5 show the state transition probability, the mean, and variance of duration in each state when dt(1)=4, representing the bearing in state 1 with duration of 4.Table 2
State transition probability(dt(1)=1).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Normal 0
0.8813
0.0454
0.0732
0.0001
Degradation 1
0.0000
0.6231
0.3473
0.0296
Degradation 2
0.0000
0.0000
0.9263
0.0737
Failure 3
0.0000
0.0000
0.0000
1.0000Table 3
Mean and variance of duration in each state(dt(1)=1).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Mean
7.9577
7.2105
7.1435
3.3435
Variance
1.3953
0.7429
0.7924
0.5452Table 4
State transition probability(dt(1)=4).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Normal 0
0.8813
0.0454
0.0732
0.0001
Degradation 1
0.0000
0.4728
0.4154
0.1118
Degradation 2
0.0000
0.0000
0.9263
0.0737
Failure 3
0.0000
0.0000
0.0000
1.0000Table 5
Mean and variance of duration in each state(dt(1)=4).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Mean
7.9577
5.3105
7.1435
3.3435
Variance
1.3953
1.1328
0.7924
0.5452First, the statei of the current operating state based on the recognition results is determined; then the residence time ∑j=i+1N-1D(j) is calculated according to the duration parameters of the operating state in training process. Then, the remaining effective life in the current operational state is calculated using (25). Finally, the RUL of the bearing can be calculated using (26). Suppose that the bearing is now at state 1 with a duration of 1; then the following can be obtained: D(2)+D(3)=10.9426, D^(11)=6.0875 by (25), and RUL1(1)=17.0211 by (26).


### 5.3. Prediction Comparison

In order to compare the prognostic method based on the DD-HSMM with the prognostic method based on the HSMM, (29) is used to evaluate the life error. In (29), RULactual represents the actual life of the component, and RULforecasted represents the expected life predicted by DD-HSMM or HSMM:
(29)Error=100%×|RULactual-RULforecasted|RULactual.Table6 shows the prediction comparison of DD-HSMM versus HSMM. Failure prediction of the HSMM method is only state dependent, while the DD-HSMM method uses both state dependency and duration dependency. The DD-HSMM method has a self-updating capability, in which the historical data on states are used in the calculation of state transition probability matrix. As indicated in the results, the DD-HSMM method is more accurate than the HSMM method.Table 6
Comparison of DD-HSMM versus HSMM.
Actual RUL
DD-HSMM model
HSMM model
Predicted RUL
Error (%)
Predicted RUL
Error (%)
27.0000
26.4027
2.212
26.9852
0.0548
24.0000
23.6596
1.418
12.438
22.0000
22.0856
0.389
22.66
17.0000
17.6142
3.613
17.0734
0.437
15.0000
15.7246
4.831
13.829
12.0000
11.6945
2.546
10.5199
12.334
11.0000
10.0602
8.544
4.365
9.0000
9.1944
2.16
16.888
5.0000
5.1253
2.506
3.6104
27.792
3.0000
3.0122
0.407
20.347


## 5.1. Operation State Identification

In order to identify the accuracy of the operation state identification method proposed in this paper, experimental data with normal operating condition were obtained. The experimental data set included 50 samples for each state (denoted by 0, 1, 2, and 3). Of these data points, 20 of them were used to train the model, and the remaining 30 samples were used to validate the model.In the DD-HSMM, mixture Gaussian distribution and the single Gaussian distribution were used to model the output probability distribution and the state duration densities separately, in which the number of states is 4. The maximum number of iterations in training process is set to 100 and the convergence error to 0.000001.The DD-HSMM-based training model is shown as Figure1. Thex-axis shows the training steps and they-axis represents the likelihood probability of different states. As can be seen from Figure 1, the progression of the four states reaches the set error in less than 40 steps. This demonstrates the potential of the model to have a strong real-time signal processing capability.Figure 1
Training curve of the DD-HSMM model.The classification results obtained on the remaining 30 data samples are shown in Table1. As indicated in the results, the accuracy of the DD-HSMM method is 94.2%.Table 1
Prognostics results.
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Recognition accuracy
Normal 0
29
1
0
0
96.7%
Degradation 1
2
27
1
0
90%
Degradation 2
0
1
28
1
93.3%
Failure  3
0
0
1
29
96.7%
Total accuracy
94.2%


## 5.2. Health Prediction for RUL

As described before, a four-state DD-HSMM prediction model is constructed. In the training process, even if the device is in the same running condition, the dwell time is different, transition probabilities between states and the mean or variance of duration in each state are not the same. Tables2 and 3 show the state transition probability, the mean, and variance of duration in each state when dt(1)=1, representing the bearing in state 1 with duration of 1. Tables 4 and 5 show the state transition probability, the mean, and variance of duration in each state when dt(1)=4, representing the bearing in state 1 with duration of 4.Table 2
State transition probability(dt(1)=1).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Normal 0
0.8813
0.0454
0.0732
0.0001
Degradation 1
0.0000
0.6231
0.3473
0.0296
Degradation 2
0.0000
0.0000
0.9263
0.0737
Failure 3
0.0000
0.0000
0.0000
1.0000Table 3
Mean and variance of duration in each state(dt(1)=1).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Mean
7.9577
7.2105
7.1435
3.3435
Variance
1.3953
0.7429
0.7924
0.5452Table 4
State transition probability(dt(1)=4).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Normal 0
0.8813
0.0454
0.0732
0.0001
Degradation 1
0.0000
0.4728
0.4154
0.1118
Degradation 2
0.0000
0.0000
0.9263
0.0737
Failure 3
0.0000
0.0000
0.0000
1.0000Table 5
Mean and variance of duration in each state(dt(1)=4).
System state
Normal 0
Degradation level 1
Degradation level 2
Failure 3
Mean
7.9577
5.3105
7.1435
3.3435
Variance
1.3953
1.1328
0.7924
0.5452First, the statei of the current operating state based on the recognition results is determined; then the residence time ∑j=i+1N-1D(j) is calculated according to the duration parameters of the operating state in training process. Then, the remaining effective life in the current operational state is calculated using (25). Finally, the RUL of the bearing can be calculated using (26). Suppose that the bearing is now at state 1 with a duration of 1; then the following can be obtained: D(2)+D(3)=10.9426, D^(11)=6.0875 by (25), and RUL1(1)=17.0211 by (26).


## 5.3. Prediction Comparison

In order to compare the prognostic method based on the DD-HSMM with the prognostic method based on the HSMM, (29) is used to evaluate the life error. In (29), RULactual represents the actual life of the component, and RULforecasted represents the expected life predicted by DD-HSMM or HSMM:
(29)Error=100%×|RULactual-RULforecasted|RULactual.Table6 shows the prediction comparison of DD-HSMM versus HSMM. Failure prediction of the HSMM method is only state dependent, while the DD-HSMM method uses both state dependency and duration dependency. The DD-HSMM method has a self-updating capability, in which the historical data on states are used in the calculation of state transition probability matrix. As indicated in the results, the DD-HSMM method is more accurate than the HSMM method.Table 6
Comparison of DD-HSMM versus HSMM.
Actual RUL
DD-HSMM model
HSMM model
Predicted RUL
Error (%)
Predicted RUL
Error (%)
27.0000
26.4027
2.212
26.9852
0.0548
24.0000
23.6596
1.418
12.438
22.0000
22.0856
0.389
22.66
17.0000
17.6142
3.613
17.0734
0.437
15.0000
15.7246
4.831
13.829
12.0000
11.6945
2.546
10.5199
12.334
11.0000
10.0602
8.544
4.365
9.0000
9.1944
2.16
16.888
5.0000
5.1253
2.506
3.6104
27.792
3.0000
3.0122
0.407
20.347


## 6. Conclusion

This paper presents a Duration-Dependent Hidden Semi-Markov Model (DD-HSMM) for prognostics. As opposed to the Hidden Semi-Markov Model (HSMM), failure prediction capability of the DD-HSMM method uses state dependency and duration dependency. The two important aspects of equipment health monitoring, which are the stages and the rate of aging, are taken into consideration in an integrated manner in the proposed DD-HSMM model. The duration-dependent state transition probability in the Hidden Semi-Markov model makes the decision-making more relevant to real world applications.In order to facilitate the computational procedure, a new forward-backward algorithm and reestimation approach are developed. By using autoregression, the interdependency between observations is established in the model. By incorporating an explicitly defined temporal structure into the model, the DD-HSMM is capable of predicting the remaining useful life of equipment more accurately.The demonstration of the proposed model is carried out using experimental data on rolling element bearings. The proposed model provides a powerful state recognition capability and very accurate results in terms of remaining useful life prediction. In order to draw general conclusion on the capabilities of the proposed DD-HSMM, more experimental data in various prognostics areas are needed.



---

*Source: 632702-2014-04-14.xml*