# Hidden Semi-Markov Models for Predictive Maintenance

Authors: Francesco Cartella; Jan Lemeire; Luca Dimiccoli; Hichem Sahli
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Realistic predictive maintenance approaches are essential for condition monitoring and predictive maintenance of industrial machines. In this work, we propose Hidden Semi-Markov Models (HSMMs) with (i) no constraints on the state duration density function and (ii) being applied to continuous or discrete observation. To deal with such a type of HSMM, we also propose modifications to the learning, inference, and prediction algorithms. Finally, automatic model selection has been made possible using the Akaike Information Criterion. This paper describes the theoretical formalization of the model as well as several experiments performed on simulated and real data with the aim of methodology validation. In all performed experiments, the model is able to correctly estimate the current state and to effectively predict the time to a predefined event with a low overall average absolute error. As a consequence, its applicability to real world settings can be beneficial, especially where in real time the Remaining Useful Lifetime (RUL) of the machine is calculated.


---

## Body

## 1. Introduction

Predictive models that are able to estimate the current condition and the Remaining Useful Lifetime of an industrial equipment are of high interest, especially for manufacturing companies, which can optimize their maintenance strategies. If we consider that the costs derived from maintenance are one of the largest parts of the operational costs [1] and that often the maintenance and operations departments comprise about 30% of the manpower [2, 3], it is not difficult to estimate the economic advantages that such innovative techniques can bring to industry. Moreover, predictive maintenance, where in real time the Remaining Useful Lifetime (RUL) of the machine is calculated, has been proven to significantly outperforms other maintenance strategies, such as corrective maintenance [4]. In this work, RUL is defined as the time, from the current moment, that the systems will fail [5].Failure, in this context, is defined as a deviation of the delivered output of a machine from the specified service requirements [6] that necessitate maintenance.Models like Support Vector Machines [7], Dynamic Bayesian Networks [8], clustering techniques [9], and data mining approaches [10] have been successfully applied to condition monitoring, RUL estimation, and predictive maintenance problems [11, 12]. State space models, like Hidden Markov Models (HMMs) [13], are particularly suitable to be used in industrial applications, due to their ability to model the latent state which represents the health condition of the machine.Classical HMMs have been applied to condition assessment [14, 15]; however, their usage in predictive maintenance has not been effective due to their intrinsic modeling of the state duration as a geometric distribution.To overcome this drawback, a modified version of HMM, which takes into account an estimate of the duration in each state, has been proposed in the works of Tobon-Mejia et al. [16–19]. Thanks to the explicit state sojourn time modeling, it has been shown that it is possible to effectively estimate the RUL for industrial equipment. However, the drawback of their proposed HMM model is that the state duration is always assumed as Gaussian distributed and the duration parameters are estimated empirically from the Viterbi path of the HMM.A complete specification of a duration model together with a set of learning and inference algorithms has been given firstly by Ferguson [20]. In his work, Ferguson allowed the underlying stochastic process of the state to be a semi-Markov chain, instead of a simple Markov chain of a HMM. Such model is referred to as Hidden Semi-Markov Model (HSMM) [21]. HSMMs and explicit duration modeles have been proven beneficial for many applications [22–25]. A complete overview of different duration model classes has been made by Yu [26]. Most state duration models, used in the literature, are nonparametric discrete distributions [27–29]. As a consequence, the number of parameters that describe the model and that have to be estimated is high, and consequently the learning procedure can be computationally expensive for real complex applications. Moreover, it is necessary to specify a priori the maximum duration allowed in each state.To alleviate the high dimensionality of the parameter space, parametric duration models have been proposed. For example, Salfner [6] proposed a generic parametric continuous distribution to model the state sojourn time. However, in their model, the observation has been assumed to be discrete and applied to recognize failure-prone observation sequence. Using continuous observation, Azimi et al. [30–32] specified an HSMM with parametric duration distribution belonging to the Gamma family and modeled the observation process by a Gaussian.Inspired by the latter two approaches, in this work we propose a generic specification of a parametric HSMM, in which no constraints are made on the model of the state duration and on the observation processes. In our approach, the state duration is modeled as a generic parametric density function. On the other hand, the observations can be modeled either as a discrete stochastic process or as continuous mixture of Gaussians. The latter has been shown to approximate, arbitrarily closely, any finite, continuous density function [33]. The proposed model can be generally used in a wide range of applications and types of data. Moreover, in this paper we introduce a new and more effective estimator of the time spent by the system in a determinate state prior to the current time. To the best of our knowledge, a part from the above referred works, the literature on HSMMs applied to prognosis and predictive maintenance for industrial machines is limited [34]. Hence, the present work aims to show the effectiveness of the proposed duration model in solving condition monitoring and RUL estimation problems.Dealing with state space models, and in particular of HSMMs, one should define the number of states and correct family of duration density, and in case of continuous observations, the adequate number of Gaussian mixtures. Such parameters play a prominent role, since the right model configuration is essential to enable an accurate modeling of the dynamic pattern and the covariance structure of the observed time series. The estimation of a satisfactory model configuration is referred to asmodel selection in literature.While several state-of-the-art approaches use expert knowledge to get insight on the model structure [15, 35, 36], an automated methodology for model selection is often required. In the literature, model selection has been deeply studied for a wide range of models. Among the existing methodologies, information based techniques have been extensively analyzed in literature with satisfactory results. Although Bayesian Information Criterion (BIC) is particularly appropriate to be used in finite mixture models [37, 38], Akaike Information Criterion (AIC) has been demonstrated to outperform BIC when applied to more complex models and when the sample size is limited [39, 40], which is the case of the target application of this paper.In this work AIC is used to estimate the correct model configuration, with the final goal of an automated HSMMs model selection, which exploits only the information available in the input data. While model selection techniques have been extensively used in the framework of Hidden Markov Models [41–43], to the best of our knowledge, the present work is the first that proposes their appliance to duration models and in particular to HSMMs.In summary, the present work contributes to condition monitoring, predictive maintenance, and RUL estimation problems by(i)
proposing a general Hidden Semi-Markov Model applicable for continuous or discrete observations and with no constraints on the density function used to model the state duration;(ii)
proposing a more effective estimator of the state duration variabledt(i), that is, the time spent by the system in the ith state, prior to current time t;(iii)
adapting the learning, inference and prediction algorithms considering the defined HSMM parameters and the proposeddt(i) estimator;(iv)
using the Akaike Information Criterion for automatic model selection.The rest of the paper is organized as follows: in Section2 we introduce the theory of the proposed HSMM together with its learning, inference, and prediction algorithms. Section 3 gives a short theoretical overview of the Akaike Information Criterion. Section 4 presents the methodology used to estimate the Remaining Useful Lifetime using the proposed HSMM. In Section 5 experimental results are discussed. The conclusion and future research directions are given in Section 6.


## 2. Hidden Semi-Markov Models

Hidden Semi-Markov Models (HSMMs) introduce the concept of variable duration, which results in a more accurate modeling power if the system being modeled shows a dependence on time.In this section we give the specification of the proposed HSMM, for which we model the state duration with a parametric state-dependent distribution. Compared to nonparametric modeling, this approach has two main advantages:(i)
the model is specified by a limited number of parameters; as a consequence, the learning procedure is computationally less expensive;(ii)
the model does not require the a priori knowledge of the maximum sojourn time allowed in each state, being inherently learnt through the duration distribution parameters.


### 2.1. Model Specification

A Hidden Semi-Markov Model is a doubly embedded stochastic model with an underlying stochastic process that is not observable (hidden) but can only be observed through another set of stochastic processes that produce the sequence of observations. HSMM allows the underlying process to be a semi-Markov chain with a variable duration or sojourn time for each state. The key concept of HSMMs is that the semi-Markov property holds for this model: while in HMMs the Markov property implies that the value of the hidden state at timet depends exclusively on its value of time t-1, in HSMMs the probability of transition from state Sj to state Si at time t depends on the duration spent in state Sj prior to time t.In the following we denote the number of states in the model asN, the individual states as S={S1,…,SN}, and the state at time t as st. The semi-Markov property can be written as
(1)Pst+1=Si∣st=Sj,…,s1=Sk=Pst+1=i∣st=j,dtj1≤i,j,k≤N,
where the duration variable dt(j) is defined as the time spent in state Sj prior to time t.Although the state duration is inherently discrete, in many studies [44, 45] it has been modeled with a continuous parametric density function. Similar to the work of Azimi et al. [30–32], in this paper, we use the discrete counterpart of the chosen parametric probability density function (pdf). With this approximation, if we denote the pdf of the sojourn time in state Si as f(x,θi), where θi represents the set of parameters of the pdf relative to the ith state, the probability that the system stays in state Si for exactly d time steps can be calculated as ∫d-1d‍f(x,θi)dx. Considering the HSMM formulation, we can generally denote the state dependent duration distributions by the set of their parameters relative to each state as Θ={θ1,…,θN}.Many related works on HSMMs [31, 32, 44, 45] consider f(x,θi) within the exponential family. In particular, Gamma distributions are often used in speech processing applications. In this work, we do not impose a type of distribution function to model the duration. The only requirement is that the duration should be modeled as a positive function, being negative durations physically meaningless.HSMMs require also the definition of a “dynamic” transition matrix, as a consequence of the semi-Markov property. Differently from the HMMs in which a constant transition probability leads to a geometric distributed state sojourn time, HSMMs explicitly define a transition matrix which, depending on the duration variable, has increasing probabilities of changing state as the time goes on. For convenience, we specify the state duration variable in a form of a vectordt with dimensions N×1 as
(2)dt=dtjifst=Sj1ifst≠Sj.
The quantity dt(j) can be easily calculated by induction from dt-1(j) as
(3)dtj=stj·st-1j·dt-1j+1,
where st(j) is 1 if st=Sj, 0 otherwise.If we assume that at timet the system is in state Si, we can formally define the duration-dependent transition matrix as Adt=[aij(dt)] with
(4)aijdt=Pst+1=Sj∣st=Si,dti1≤i,j≤N.The specification of the model can be further simplified by observing that, at each timet, the matrix Adt can be decomposed in two terms: the recurrent and the nonrecurrent state transition probabilities.The recurrent transition probabilitiesP(dt)=[pij(dt)], which depend only on the duration vector dt and the parameters Θ, take into account the dynamics of the self-transition probabilities. It is defined as the probability of remaining in the current state at the next time step, given the duration spent in the current state prior to time t:
(5)piidt=Pst+1=Si∣st=Si,dti=Pst+1=Si∣st-dti≠Sist=Si,st-1=Si,…,st-dti+1=Si,st-dti≠Si=Pst+1=Si,st=Si,…,st-dti+2=Si∣st-dti+1=Si,st-dti≠Si·Pst=Si,st-1=Si,…,st-dti+2=Si∣st-dti+1=Si,st-dti≠Si-1.
The denominator in (5) can be expressed as ∑k=1∞‍P(st+k≠Si,st+k-1=Si,…,st-dti+2=Si∣st-dt(i)+1=Si,st-dt(i)≠Si), which is the probability that the system, at time t, has been staying in state Si for at least dt(i)-1 time units. The above expression is equivalent to 1-F(dt(i)-1,θi), where F(·,θi) is the duration cumulative distribution function relative to the the state Si, that is, F(d,θ)=∫-∞d‍f(x,θ)dx. As a consequence, from (5) we can define the recurrent transition probabilities as a diagonal matrix with dimensions N×N, as
(6)P(dt)=[pij(dt)]=1-F(dt(i),θi)1-F(dt(i)-1,θi)ifi=j0ifi≠j.
The usage of the cumulative functions in (6), which tend to 1 as the duration tends to infinity, suggests that the probability of self-transition tends to decrease as the sojourn time increases, leading the model to always leave the current state if time approaches infinity.The nonrecurrent state transition probabilities,A0=[aij0], rule the transitions between two different states. It is represented by a N×N matrix with the diagonal elements equal to zero, defined as
(7)A0=[aij0]=0ifi=jPst+1=Sj∣st=Siifi≠j.A0 must be specified as a stochastic matrix; that is, its elements have to satisfy the constraint ∑j=1N‍aij0=1 for all i.As a consequence of the above decomposition, the dynamic of the underlying semi-Markov chain can be defined by specifying only the state-dependent duration parametersΘ and the nonrecurrent matrix A0, since the model transition matrix can be calculated, at each time t, using (6) and (7):
(8)Adt=Pdt+I-PdtA0,
where I is the identity matrix. If we denote the elements of the dynamic transition matrix Adt as aij(dt), the stochastic constraint ∑j=1N‍aij(dt)=1 for all i and t is guaranteed from the fact that P(dt) is a diagonal matrix and A0 is a stochastic matrix.For several applications it is necessary to model theabsorbing state which, in the case of industrial equipment, corresponds to the “broken” or “failure” state. If we denote the absorbing state as Sk with k∈[1,N], we must fix the kth row of the nonrecurrent matrix A0 to be akk0=1 and aki0=0 for all 1≤i≤N with i≠k. By substituting such A0 matrix in (8), it is easy to show that the element akk(dt)=1 and remains constant for all t, while the duration probability parameters θk are not influent for the absorbing state Sk. An example of absorbing state specification will be given in Section 5.With respect to the input observation signals, in this work we consider both continuous and discrete data, by adapting the suitable observation model depending on the observation nature. In particular, for the continuous case, we model the observations with a multivariate mixture of Gaussians distributions. This choice presents two main advantages: (i) a multivariate model allows to deal with multiple observations at the same time; this is often the case of industrial equipments modeling since, at each time, multiple sensors’ measurements are available, and (ii) mixture of Gaussians has been proved to closely approximate any finite and continuous density function [33]. Formally, if we denote by xt the observation vector at time t and the generic observation vector being modeled as x, the observation density for the jth state is represented by a finite mixture of M gaussians
(9)bjx=∑m=1M‍cjmNx,μjm,Ujm,1≤j≤N,
where cjm is the mixture coefficient for the mth mixture in state Sj, which satisfies the stochastic constraint ∑m=1M‍cjm=1 for 1≤j≤N and cjm≥0 for 1≤j≤N and 1≤m≤M, while N is the Gaussian density, with mean vector μjm and covariance matrix Ujm for the mth mixture component in state j.In case of discrete data, we model the observations within each state with a nonparametric discrete probability distribution. In particular, ifL is the number of distinct observation symbols per state and if we denote the symbols as X={X1,…,XL} and the observation at time t as xt, the observation symbol probability distribution can be defined as a matrix B=[bj(l)] of dimensions N×L where
(10)bj(l)=Pxt=Xl∣st=Sj1≤j≤N;1≤l≤L.
Since the system in each state at each time step can emit one of the possible L symbols, the matrix B is stochastic; that is, it is constrained to ∑l=1L‍bj(l)=1 for all 1≤j≤N.Finally, as in the case of HMMs, we specify the initial state distributionπ={πi} which defines the probability of the starting state as
(11)πi=P[s1=Si],1≤i≤N.From the above considerations, two different HSMM models can be considered. In the case of continuous observation,λ=(A0,Θ,C,μ,U,π), and in the case of discrete observation the HSMM is characterized by λ=(A0,Θ,B,π). An example of continuous HSMM with 3 states is shown in Figure 1.Figure 1
Graphical representation of an HSMM.


### 2.2. Learning and Inference Algorithms

Let us denote the generic sequence of observations, being indiscriminately continuous vectors or discrete symbols, asx=x1x2⋯xT; in order to use the defined HSMM model in practice, similarly to the HMM, we need to solve three basic problems.(1)
Given the observationx and a model λ, calculate the probability that the sequence x has been generated by the model λ, that is, P(x∣λ).(2)
Given the observationx and a model λ, calculate the state sequence S=s1s2⋯sT which have most probably generated the sequence x.(3)
Given the observationx find the parameters of the model λ which maximize P(x∣λ).As in case of HMM, solving the above problems requires using the forward-backward [13], decoding (Viterbi [46] and Forney [47]) and Expectation Maximization [48] algorithms, which will be adapted to the HSMM introduced in Section 2.1. In the following, we also propose a more effective estimator of the state duration variable dt(i) defined in (2).


#### 2.2.1. The Forward-Backward Algorithm

Given a generic sequence of observationsx=x1x2⋯xT, the goal is to calculate the model likelihood, that is, P(x∣λ). This quantity is useful for the training procedure, where the parameters that locally maximize the model likelihood are chosen, as well as for classification problems. The latter is the case in which the observation sequence x has to be mapped to one of a finite set of C classes, represented by a set of HSMM parameters L={λ1,…,λC}. The class of x is chosen such that λ(x)=argmax⁡λ∈LP(X∣λ).To calculate the model likelihood, we first define the forward variable at each timet as
(12)αt(i)=Px1x2⋯xt,st=Si∣λ1≤i≤N.Contrarily to HMMs, for HSMMs the state duration must be taken into account in the the forward variable calculation. Consequently, Yu [26] proposed the following inductive formula:
(13)αtj,d=∑d′=1D‍∑i=1N‍αt-d′i,d′aij0pjjd′∏k=t-d+1tbjxkiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤j≤N,1≤t≤T
that is, the sum of the probabilities of being in the current state Sj (at time t) for the past d′ time units, (with 1≤d′≤D and D the maximum allowed duration for each state), coming from all the possible previous states Si, 1≤i≤N, and i≠j.The disadvantage of the above formulation is that, as discussed in Introduction, the specification of the maximum durationD represents a limitation to the modeling generalization. Moreover, from (13), it is clear that the computation and memory complexities drastically increase with D, which can be very large in many applications, in particular for online failure prediction.To alleviate this problem, Azimi et al. [30–32] introduced a new forward algorithm for HSMMs that, by keeping track of the estimated average state duration at each iteration, has a computational complexity comparable to the forward algorithm for HMMs [13]. However, the average state duration represents an approximation. Consequently, the forward algorithm of Azimi, compared with (13), pays the price of a lower precision in favor of a (indispensable) better computational efficiency.To calculate the forward variableαt(j) using Azimi’s approach, the duration-dependent transition matrix, defined in (8), is taken in consideration in the induction formula of (13), which becomes [30]
(14)αtj=∑i=1N‍αt-1iaijd^t-1bjxt.
To calculate the above formula, the average state duration of (2) must be estimated, for each time t, by means of the variable d^t=[d^t(i)], defined as
(15)d^ti=Edti∣x1x2⋯xt,st=Si,λ1≤i≤N,
where E denotes the expected value. To calculate the above quantity, Azimi et al. [30–32] use the following formula:
(16)d^t=γt-1⊙d^t-1+1,
where ⊙ represents the element by element product between two matrices/vectors and the vector γt=γt(i) (the probability of being in state Si at time t given the observation sequence and the model parameters) with dimensions N×1 is calculated in terms of αt(i) as
(17)γti=Pst=Si∣x1x2⋯xt,λ=αti∑j=1N‍αtjiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.Equation (16) is based on the following induction formula [30–32] that rules the dynamics of the duration vector when the system’s state is known:
(18)dt(i)=st-1(i)·dt-1(i)+1,
where, for each t, st(i) is 1 if st=Si, 0 otherwise.A simple example shows that (18) is incorrect: assuming an HSMM with three states and considering the state sequence (S1,S1,S2,…) the correct sequence of the duration vector is d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,1,1]T, where the superscript T denotes vector transpose. If we apply (18), we obtain d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,2,1]T, which is in contradiction with the definition of the state duration vector given in (2).To calculate the average state duration variabled^t(i) we propose a new induction formula that estimates, for each time t, the time spent in the ith state prior to t as
(19)d^t(i)=Pst-1=Si∣st=Si,x1,…,xt·d^t-1i+1(20)=aiid^t-1·αt-1(i)·bixtαti·d^t-1i+1,iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.
The derivation of (20) is given in Appendix. The intuition behind (19) is that the current average duration is the previous average duration plus one, weighted with the “amount” of the current state that was already in state Si in the previous step.Using the proposed (20), the forward algorithm can be specified as follows:(1)
initialization, with1≤i≤N:
(21)α1i=πibix1,d^1(i)=1,Ad^1=Pd^1+I-Pd^1A0,whereP(d^i) is estimated using (6);(2)
induction, with1≤j≤N and 1≤t≤T-1:
(22)αt+1j=∑i=1N‍αtiaijd^tbjxt+1,(23)d^t+1i=aiid^t·αti·bixt+1αt+1i·d^ti+1,(24)Ad^t+1=Pd^t+1+I-Pd^t+1A0,whereaij(d^t) are the coefficients of the matrix Ad^t;(3)
termination:(25)Px∣λ=∑i=1N‍αTi.Similar considerations as the forward procedure can be made for the backward algorithm, which is implemented by defining the variableβt(i) as
(26)βt(i)=Pxt+1xt+2⋯xT∣st=Si,λ1≤i≤N.Having estimated the dynamic transition matrixAd^t for each 1≤t≤T using (24), the backward variable can be calculated inductively as follows.(1)
Initialization:(27)βT(i)=1,1≤i≤N.(2)
Induction:(28)βti=∑j=1N‍aijd^tbjxt+1βt+1j,it=T-1,T-2,…,1,1≤i≤N.Although the variableβt(i) is not necessary for the calculation of the model likelihood, it will be useful in the parameter reestimation procedure, as it will be explained in Section 2.2.3.


#### 2.2.2. The Viterbi Algorithm

The Viterbi algorithm [46, 47] (also known asdecoding) allows determining the best state sequence corresponding to a given observation sequence.Formally, given a sequence of observationx=x1x2⋯xT, the best state sequence S*=s1*s2*⋯sT* corresponding to x is calculated by defining the variable δt(i) as
(29)δt(i)=max⁡s1,s2,…,st-1Ps1s2,…,st=Si,x1x2⋯xt∣λ.The procedure to recursively calculate the variableδt(i) and to retrieve the target state sequence (i.e., the arguments which maximize the δt(i)’s) for the proposed HSMM is a straightforward extension of the Viterbi algorithm for HMMs [13]. The only change is the usage, in the recursive calculation of δt(i), of the dynamic transition matrix Ad^t=[aij(d^t)], calculated through (24). The Viterbi algorithm for the introduced parametric HSMMs can be summarized as follows:(1)
initialization, with1≤i≤N:
(30)δ1i=πibix1,ψ1i=0,(2)
recursion, with1≤j≤N and 2≤t≤T:
(31)δtj=max⁡1≤i≤Nδt-1iaijd^tbjxt,(32)ψtj=argmax1≤i≤N⁡δt-1(i)aij(d^t),(3)
termination:(33)P*=max⁡1≤i≤NδT(i),(34)sT*=argmax1≤i≤N⁡δT(i),where we keep track of the argument maximizing (31) using the vector ψt, which, tracked back, gives the desired best state sequence:
(35)st*=ψt+1st+1*,t=T-1,T-2,…,1.


#### 2.2.3. The Training Algorithm

The training algorithm consists of estimating the model parameters from the observation data. As discussed in Section2.1, a parametric HSMM is defined by λ=(A0,Θ,C,μ,U,π) if the observations are continuous, or λ=(A0,Θ,B,π) if the observations are discrete. Given a generic observation sequence x=x1x2⋯xT, referred to astraining set in the following, the training procedure consists of finding the model parameter set λ* which locally maximizes the model likelihood Px∣λ.We use the modified Baum-Welch algorithm of Azimi et al. [30–32]. However, in our implementation we do not make assumption on the density function used to model the state duration, and we consider both continuous and discrete observations.Being a variant of the more general Expectation-Maximization (EM) algorithm, Baum-Welch is an iterative procedure which consists of two steps: (i) the expectation step, in which the forward and backward variables are calculated and the model likelihood is estimated and (ii) the maximization step, in which the model parameters are updated and used in the next iteration. This process usually starts from a random guess of the model parametersλ0 and it is iterated until the likelihood function does not improve between two consecutive iterations.Similarly to HMMs, the reestimation formulas are derived by firstly introducing the variableξt(i,j), which represents the probability of being in state Si at time t, and in state Sj at time t+1, given the model and the observation sequence, as
(36)ξti,j=Pst=Si,st+1=Sj∣x,λ.
However, in the HSMM case, the variable ξt(i,j) considers the duration estimation performed in the forward algorithm (see Equation (24)). Formulated in terms of the forward and backward variables, it is given by
(37)ξti,j=Pst=Si,st+1=Sj∣x,λ=Pst=Si,st+1=Sj,x∣λPx∣λ=αtiaijd^tbjxt+1βt+1jPx∣λ=αtiaijd^tbjxt+1βt+1j∑i=1N‍∑j=1N‍αtiaijd^tbjxt+1βt+1j.Fromξt(i,j) we can derive the quantity γt(i) (already defined in (17)) representing the probability of being in state Si at time t given the observation sequence and the model parameters:
(38)γt(i)=∑j=1N‍ξt(i,j).Finally, the the reestimation formulas for the parametersπ and A0 are given by
(39)π¯i=γ1i,(40)a¯ij0=∑t=1T-1‍ξt(i,j)⊙G∑j=1N‍∑t=1T-1‍ξt(i,j)⊙G,
where G=[gij] is a square matrix of dimensions N×N where gij=0 for i=j and gij=1 for i≠j, ⊙ represents the element by element product between two matrices, ∑t=1T-1‍γt(i) is the expected number of transitions from state Si, and ∑t=1T-1‍ξt(i,j) is the expected number of transitions from state Si to state Sj.Equation (39) represents the expected number of times that the model starts in state Si, while (40) represents the expected number of transitions from state Si to state Sj with i≠j over the total expected number of transitions from state Si to any other state different from Si.For the matrixA¯0, being normalized, the stochastic constraints are satisfied at each iteration, that is, ∑j=1N‍a¯ij0=1 for each 1≤i≤N, while the estimation of the prior probability π¯i inherently sums up to 1 at each iteration, since it represents the expected frequency in state Si at time t=1 for each 1≤i≤N.With respect to the reestimation of the state duration parametersΘ, firstly, we estimate the mean μi,d and the variance σi,d2 of the ith state duration for each 1≤i≤N, from the forward and backward variables and the estimation of the state duration variable
(41)μi,d=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1jd^ti∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j,(42)σi,d2=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j∑j=1,j≠iN‍·d^t(i)-μi,d2·∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j-1,
where (41) can be interpreted as the probability of transition from state Si to Sj with i≠j at time t weighted by the duration of state Si at t, giving the desired expected value, while in (42) the same quantity is weighted by the squared distance of the duration at time t from its mean, giving the estimation of the variance.Then, the parameters of the desired duration distribution can be estimated fromμi,d and σi,d2. For example, if a Gamma distribution with shape parameter ν and scale parameter η is chosen to model the state duration, the parameters νi and ηi for each 1≤i≤N can be calculated as νi=μi,d2/σi,d2 and ηi=σi,d2/μi,d.Concerning the observation parameters, once the modified forward and backward variables, accounting for the state duration, are defined as in (22) and (28), the reestimation formulas are the same as for Hidden Markov Models [13].In particular, for continuous observations, the parameters of the Gaussians’ mixture defined in (9) are reestimated by firstly defining the probability of being in state Sj at time t with the probability of the observation vector xt evaluated by the kth mixture component, as
(43)γtj,k=αt(j)βt(j)∑j=1N‍αt(j)βt(j)·cjkNxt,μjk,Ujk∑m=1M‍cjmNxt,μjm,Ujm.
By using the former quantity, the parameters cjk, μjk, and Ujk are reestimated through the following formulas:
(44)c¯jk=∑t=1T‍γt(j,k)∑t=1T‍∑m=1M‍γt(j,m),μ¯jk=∑t=1T‍γtj,k·xt∑t=1T‍γt(j,k),U¯jk=∑t=1T‍γt(j,k)·(xt-μjk)xt-μjkT∑t=1T‍γt(j,k),
where superscript T denotes vector transpose.For discrete observations, the reestimation formula for the observation matrixbj(l) is
(45)b¯jl=∑t=1,withxt=XlT‍γtj∑t=1T‍γtj,
where the quantity γt(j), which takes into account the duration dependent forward variable αt(j), is calculated through (17).The reader is referred to Rabiner’s work [13] for the interpretation on the observation parameters reestimation formulas.


## 2.1. Model Specification

A Hidden Semi-Markov Model is a doubly embedded stochastic model with an underlying stochastic process that is not observable (hidden) but can only be observed through another set of stochastic processes that produce the sequence of observations. HSMM allows the underlying process to be a semi-Markov chain with a variable duration or sojourn time for each state. The key concept of HSMMs is that the semi-Markov property holds for this model: while in HMMs the Markov property implies that the value of the hidden state at timet depends exclusively on its value of time t-1, in HSMMs the probability of transition from state Sj to state Si at time t depends on the duration spent in state Sj prior to time t.In the following we denote the number of states in the model asN, the individual states as S={S1,…,SN}, and the state at time t as st. The semi-Markov property can be written as
(1)Pst+1=Si∣st=Sj,…,s1=Sk=Pst+1=i∣st=j,dtj1≤i,j,k≤N,
where the duration variable dt(j) is defined as the time spent in state Sj prior to time t.Although the state duration is inherently discrete, in many studies [44, 45] it has been modeled with a continuous parametric density function. Similar to the work of Azimi et al. [30–32], in this paper, we use the discrete counterpart of the chosen parametric probability density function (pdf). With this approximation, if we denote the pdf of the sojourn time in state Si as f(x,θi), where θi represents the set of parameters of the pdf relative to the ith state, the probability that the system stays in state Si for exactly d time steps can be calculated as ∫d-1d‍f(x,θi)dx. Considering the HSMM formulation, we can generally denote the state dependent duration distributions by the set of their parameters relative to each state as Θ={θ1,…,θN}.Many related works on HSMMs [31, 32, 44, 45] consider f(x,θi) within the exponential family. In particular, Gamma distributions are often used in speech processing applications. In this work, we do not impose a type of distribution function to model the duration. The only requirement is that the duration should be modeled as a positive function, being negative durations physically meaningless.HSMMs require also the definition of a “dynamic” transition matrix, as a consequence of the semi-Markov property. Differently from the HMMs in which a constant transition probability leads to a geometric distributed state sojourn time, HSMMs explicitly define a transition matrix which, depending on the duration variable, has increasing probabilities of changing state as the time goes on. For convenience, we specify the state duration variable in a form of a vectordt with dimensions N×1 as
(2)dt=dtjifst=Sj1ifst≠Sj.
The quantity dt(j) can be easily calculated by induction from dt-1(j) as
(3)dtj=stj·st-1j·dt-1j+1,
where st(j) is 1 if st=Sj, 0 otherwise.If we assume that at timet the system is in state Si, we can formally define the duration-dependent transition matrix as Adt=[aij(dt)] with
(4)aijdt=Pst+1=Sj∣st=Si,dti1≤i,j≤N.The specification of the model can be further simplified by observing that, at each timet, the matrix Adt can be decomposed in two terms: the recurrent and the nonrecurrent state transition probabilities.The recurrent transition probabilitiesP(dt)=[pij(dt)], which depend only on the duration vector dt and the parameters Θ, take into account the dynamics of the self-transition probabilities. It is defined as the probability of remaining in the current state at the next time step, given the duration spent in the current state prior to time t:
(5)piidt=Pst+1=Si∣st=Si,dti=Pst+1=Si∣st-dti≠Sist=Si,st-1=Si,…,st-dti+1=Si,st-dti≠Si=Pst+1=Si,st=Si,…,st-dti+2=Si∣st-dti+1=Si,st-dti≠Si·Pst=Si,st-1=Si,…,st-dti+2=Si∣st-dti+1=Si,st-dti≠Si-1.
The denominator in (5) can be expressed as ∑k=1∞‍P(st+k≠Si,st+k-1=Si,…,st-dti+2=Si∣st-dt(i)+1=Si,st-dt(i)≠Si), which is the probability that the system, at time t, has been staying in state Si for at least dt(i)-1 time units. The above expression is equivalent to 1-F(dt(i)-1,θi), where F(·,θi) is the duration cumulative distribution function relative to the the state Si, that is, F(d,θ)=∫-∞d‍f(x,θ)dx. As a consequence, from (5) we can define the recurrent transition probabilities as a diagonal matrix with dimensions N×N, as
(6)P(dt)=[pij(dt)]=1-F(dt(i),θi)1-F(dt(i)-1,θi)ifi=j0ifi≠j.
The usage of the cumulative functions in (6), which tend to 1 as the duration tends to infinity, suggests that the probability of self-transition tends to decrease as the sojourn time increases, leading the model to always leave the current state if time approaches infinity.The nonrecurrent state transition probabilities,A0=[aij0], rule the transitions between two different states. It is represented by a N×N matrix with the diagonal elements equal to zero, defined as
(7)A0=[aij0]=0ifi=jPst+1=Sj∣st=Siifi≠j.A0 must be specified as a stochastic matrix; that is, its elements have to satisfy the constraint ∑j=1N‍aij0=1 for all i.As a consequence of the above decomposition, the dynamic of the underlying semi-Markov chain can be defined by specifying only the state-dependent duration parametersΘ and the nonrecurrent matrix A0, since the model transition matrix can be calculated, at each time t, using (6) and (7):
(8)Adt=Pdt+I-PdtA0,
where I is the identity matrix. If we denote the elements of the dynamic transition matrix Adt as aij(dt), the stochastic constraint ∑j=1N‍aij(dt)=1 for all i and t is guaranteed from the fact that P(dt) is a diagonal matrix and A0 is a stochastic matrix.For several applications it is necessary to model theabsorbing state which, in the case of industrial equipment, corresponds to the “broken” or “failure” state. If we denote the absorbing state as Sk with k∈[1,N], we must fix the kth row of the nonrecurrent matrix A0 to be akk0=1 and aki0=0 for all 1≤i≤N with i≠k. By substituting such A0 matrix in (8), it is easy to show that the element akk(dt)=1 and remains constant for all t, while the duration probability parameters θk are not influent for the absorbing state Sk. An example of absorbing state specification will be given in Section 5.With respect to the input observation signals, in this work we consider both continuous and discrete data, by adapting the suitable observation model depending on the observation nature. In particular, for the continuous case, we model the observations with a multivariate mixture of Gaussians distributions. This choice presents two main advantages: (i) a multivariate model allows to deal with multiple observations at the same time; this is often the case of industrial equipments modeling since, at each time, multiple sensors’ measurements are available, and (ii) mixture of Gaussians has been proved to closely approximate any finite and continuous density function [33]. Formally, if we denote by xt the observation vector at time t and the generic observation vector being modeled as x, the observation density for the jth state is represented by a finite mixture of M gaussians
(9)bjx=∑m=1M‍cjmNx,μjm,Ujm,1≤j≤N,
where cjm is the mixture coefficient for the mth mixture in state Sj, which satisfies the stochastic constraint ∑m=1M‍cjm=1 for 1≤j≤N and cjm≥0 for 1≤j≤N and 1≤m≤M, while N is the Gaussian density, with mean vector μjm and covariance matrix Ujm for the mth mixture component in state j.In case of discrete data, we model the observations within each state with a nonparametric discrete probability distribution. In particular, ifL is the number of distinct observation symbols per state and if we denote the symbols as X={X1,…,XL} and the observation at time t as xt, the observation symbol probability distribution can be defined as a matrix B=[bj(l)] of dimensions N×L where
(10)bj(l)=Pxt=Xl∣st=Sj1≤j≤N;1≤l≤L.
Since the system in each state at each time step can emit one of the possible L symbols, the matrix B is stochastic; that is, it is constrained to ∑l=1L‍bj(l)=1 for all 1≤j≤N.Finally, as in the case of HMMs, we specify the initial state distributionπ={πi} which defines the probability of the starting state as
(11)πi=P[s1=Si],1≤i≤N.From the above considerations, two different HSMM models can be considered. In the case of continuous observation,λ=(A0,Θ,C,μ,U,π), and in the case of discrete observation the HSMM is characterized by λ=(A0,Θ,B,π). An example of continuous HSMM with 3 states is shown in Figure 1.Figure 1
Graphical representation of an HSMM.


## 2.2. Learning and Inference Algorithms

Let us denote the generic sequence of observations, being indiscriminately continuous vectors or discrete symbols, asx=x1x2⋯xT; in order to use the defined HSMM model in practice, similarly to the HMM, we need to solve three basic problems.(1)
Given the observationx and a model λ, calculate the probability that the sequence x has been generated by the model λ, that is, P(x∣λ).(2)
Given the observationx and a model λ, calculate the state sequence S=s1s2⋯sT which have most probably generated the sequence x.(3)
Given the observationx find the parameters of the model λ which maximize P(x∣λ).As in case of HMM, solving the above problems requires using the forward-backward [13], decoding (Viterbi [46] and Forney [47]) and Expectation Maximization [48] algorithms, which will be adapted to the HSMM introduced in Section 2.1. In the following, we also propose a more effective estimator of the state duration variable dt(i) defined in (2).


### 2.2.1. The Forward-Backward Algorithm

Given a generic sequence of observationsx=x1x2⋯xT, the goal is to calculate the model likelihood, that is, P(x∣λ). This quantity is useful for the training procedure, where the parameters that locally maximize the model likelihood are chosen, as well as for classification problems. The latter is the case in which the observation sequence x has to be mapped to one of a finite set of C classes, represented by a set of HSMM parameters L={λ1,…,λC}. The class of x is chosen such that λ(x)=argmax⁡λ∈LP(X∣λ).To calculate the model likelihood, we first define the forward variable at each timet as
(12)αt(i)=Px1x2⋯xt,st=Si∣λ1≤i≤N.Contrarily to HMMs, for HSMMs the state duration must be taken into account in the the forward variable calculation. Consequently, Yu [26] proposed the following inductive formula:
(13)αtj,d=∑d′=1D‍∑i=1N‍αt-d′i,d′aij0pjjd′∏k=t-d+1tbjxkiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤j≤N,1≤t≤T
that is, the sum of the probabilities of being in the current state Sj (at time t) for the past d′ time units, (with 1≤d′≤D and D the maximum allowed duration for each state), coming from all the possible previous states Si, 1≤i≤N, and i≠j.The disadvantage of the above formulation is that, as discussed in Introduction, the specification of the maximum durationD represents a limitation to the modeling generalization. Moreover, from (13), it is clear that the computation and memory complexities drastically increase with D, which can be very large in many applications, in particular for online failure prediction.To alleviate this problem, Azimi et al. [30–32] introduced a new forward algorithm for HSMMs that, by keeping track of the estimated average state duration at each iteration, has a computational complexity comparable to the forward algorithm for HMMs [13]. However, the average state duration represents an approximation. Consequently, the forward algorithm of Azimi, compared with (13), pays the price of a lower precision in favor of a (indispensable) better computational efficiency.To calculate the forward variableαt(j) using Azimi’s approach, the duration-dependent transition matrix, defined in (8), is taken in consideration in the induction formula of (13), which becomes [30]
(14)αtj=∑i=1N‍αt-1iaijd^t-1bjxt.
To calculate the above formula, the average state duration of (2) must be estimated, for each time t, by means of the variable d^t=[d^t(i)], defined as
(15)d^ti=Edti∣x1x2⋯xt,st=Si,λ1≤i≤N,
where E denotes the expected value. To calculate the above quantity, Azimi et al. [30–32] use the following formula:
(16)d^t=γt-1⊙d^t-1+1,
where ⊙ represents the element by element product between two matrices/vectors and the vector γt=γt(i) (the probability of being in state Si at time t given the observation sequence and the model parameters) with dimensions N×1 is calculated in terms of αt(i) as
(17)γti=Pst=Si∣x1x2⋯xt,λ=αti∑j=1N‍αtjiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.Equation (16) is based on the following induction formula [30–32] that rules the dynamics of the duration vector when the system’s state is known:
(18)dt(i)=st-1(i)·dt-1(i)+1,
where, for each t, st(i) is 1 if st=Si, 0 otherwise.A simple example shows that (18) is incorrect: assuming an HSMM with three states and considering the state sequence (S1,S1,S2,…) the correct sequence of the duration vector is d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,1,1]T, where the superscript T denotes vector transpose. If we apply (18), we obtain d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,2,1]T, which is in contradiction with the definition of the state duration vector given in (2).To calculate the average state duration variabled^t(i) we propose a new induction formula that estimates, for each time t, the time spent in the ith state prior to t as
(19)d^t(i)=Pst-1=Si∣st=Si,x1,…,xt·d^t-1i+1(20)=aiid^t-1·αt-1(i)·bixtαti·d^t-1i+1,iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.
The derivation of (20) is given in Appendix. The intuition behind (19) is that the current average duration is the previous average duration plus one, weighted with the “amount” of the current state that was already in state Si in the previous step.Using the proposed (20), the forward algorithm can be specified as follows:(1)
initialization, with1≤i≤N:
(21)α1i=πibix1,d^1(i)=1,Ad^1=Pd^1+I-Pd^1A0,whereP(d^i) is estimated using (6);(2)
induction, with1≤j≤N and 1≤t≤T-1:
(22)αt+1j=∑i=1N‍αtiaijd^tbjxt+1,(23)d^t+1i=aiid^t·αti·bixt+1αt+1i·d^ti+1,(24)Ad^t+1=Pd^t+1+I-Pd^t+1A0,whereaij(d^t) are the coefficients of the matrix Ad^t;(3)
termination:(25)Px∣λ=∑i=1N‍αTi.Similar considerations as the forward procedure can be made for the backward algorithm, which is implemented by defining the variableβt(i) as
(26)βt(i)=Pxt+1xt+2⋯xT∣st=Si,λ1≤i≤N.Having estimated the dynamic transition matrixAd^t for each 1≤t≤T using (24), the backward variable can be calculated inductively as follows.(1)
Initialization:(27)βT(i)=1,1≤i≤N.(2)
Induction:(28)βti=∑j=1N‍aijd^tbjxt+1βt+1j,it=T-1,T-2,…,1,1≤i≤N.Although the variableβt(i) is not necessary for the calculation of the model likelihood, it will be useful in the parameter reestimation procedure, as it will be explained in Section 2.2.3.


### 2.2.2. The Viterbi Algorithm

The Viterbi algorithm [46, 47] (also known asdecoding) allows determining the best state sequence corresponding to a given observation sequence.Formally, given a sequence of observationx=x1x2⋯xT, the best state sequence S*=s1*s2*⋯sT* corresponding to x is calculated by defining the variable δt(i) as
(29)δt(i)=max⁡s1,s2,…,st-1Ps1s2,…,st=Si,x1x2⋯xt∣λ.The procedure to recursively calculate the variableδt(i) and to retrieve the target state sequence (i.e., the arguments which maximize the δt(i)’s) for the proposed HSMM is a straightforward extension of the Viterbi algorithm for HMMs [13]. The only change is the usage, in the recursive calculation of δt(i), of the dynamic transition matrix Ad^t=[aij(d^t)], calculated through (24). The Viterbi algorithm for the introduced parametric HSMMs can be summarized as follows:(1)
initialization, with1≤i≤N:
(30)δ1i=πibix1,ψ1i=0,(2)
recursion, with1≤j≤N and 2≤t≤T:
(31)δtj=max⁡1≤i≤Nδt-1iaijd^tbjxt,(32)ψtj=argmax1≤i≤N⁡δt-1(i)aij(d^t),(3)
termination:(33)P*=max⁡1≤i≤NδT(i),(34)sT*=argmax1≤i≤N⁡δT(i),where we keep track of the argument maximizing (31) using the vector ψt, which, tracked back, gives the desired best state sequence:
(35)st*=ψt+1st+1*,t=T-1,T-2,…,1.


### 2.2.3. The Training Algorithm

The training algorithm consists of estimating the model parameters from the observation data. As discussed in Section2.1, a parametric HSMM is defined by λ=(A0,Θ,C,μ,U,π) if the observations are continuous, or λ=(A0,Θ,B,π) if the observations are discrete. Given a generic observation sequence x=x1x2⋯xT, referred to astraining set in the following, the training procedure consists of finding the model parameter set λ* which locally maximizes the model likelihood Px∣λ.We use the modified Baum-Welch algorithm of Azimi et al. [30–32]. However, in our implementation we do not make assumption on the density function used to model the state duration, and we consider both continuous and discrete observations.Being a variant of the more general Expectation-Maximization (EM) algorithm, Baum-Welch is an iterative procedure which consists of two steps: (i) the expectation step, in which the forward and backward variables are calculated and the model likelihood is estimated and (ii) the maximization step, in which the model parameters are updated and used in the next iteration. This process usually starts from a random guess of the model parametersλ0 and it is iterated until the likelihood function does not improve between two consecutive iterations.Similarly to HMMs, the reestimation formulas are derived by firstly introducing the variableξt(i,j), which represents the probability of being in state Si at time t, and in state Sj at time t+1, given the model and the observation sequence, as
(36)ξti,j=Pst=Si,st+1=Sj∣x,λ.
However, in the HSMM case, the variable ξt(i,j) considers the duration estimation performed in the forward algorithm (see Equation (24)). Formulated in terms of the forward and backward variables, it is given by
(37)ξti,j=Pst=Si,st+1=Sj∣x,λ=Pst=Si,st+1=Sj,x∣λPx∣λ=αtiaijd^tbjxt+1βt+1jPx∣λ=αtiaijd^tbjxt+1βt+1j∑i=1N‍∑j=1N‍αtiaijd^tbjxt+1βt+1j.Fromξt(i,j) we can derive the quantity γt(i) (already defined in (17)) representing the probability of being in state Si at time t given the observation sequence and the model parameters:
(38)γt(i)=∑j=1N‍ξt(i,j).Finally, the the reestimation formulas for the parametersπ and A0 are given by
(39)π¯i=γ1i,(40)a¯ij0=∑t=1T-1‍ξt(i,j)⊙G∑j=1N‍∑t=1T-1‍ξt(i,j)⊙G,
where G=[gij] is a square matrix of dimensions N×N where gij=0 for i=j and gij=1 for i≠j, ⊙ represents the element by element product between two matrices, ∑t=1T-1‍γt(i) is the expected number of transitions from state Si, and ∑t=1T-1‍ξt(i,j) is the expected number of transitions from state Si to state Sj.Equation (39) represents the expected number of times that the model starts in state Si, while (40) represents the expected number of transitions from state Si to state Sj with i≠j over the total expected number of transitions from state Si to any other state different from Si.For the matrixA¯0, being normalized, the stochastic constraints are satisfied at each iteration, that is, ∑j=1N‍a¯ij0=1 for each 1≤i≤N, while the estimation of the prior probability π¯i inherently sums up to 1 at each iteration, since it represents the expected frequency in state Si at time t=1 for each 1≤i≤N.With respect to the reestimation of the state duration parametersΘ, firstly, we estimate the mean μi,d and the variance σi,d2 of the ith state duration for each 1≤i≤N, from the forward and backward variables and the estimation of the state duration variable
(41)μi,d=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1jd^ti∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j,(42)σi,d2=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j∑j=1,j≠iN‍·d^t(i)-μi,d2·∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j-1,
where (41) can be interpreted as the probability of transition from state Si to Sj with i≠j at time t weighted by the duration of state Si at t, giving the desired expected value, while in (42) the same quantity is weighted by the squared distance of the duration at time t from its mean, giving the estimation of the variance.Then, the parameters of the desired duration distribution can be estimated fromμi,d and σi,d2. For example, if a Gamma distribution with shape parameter ν and scale parameter η is chosen to model the state duration, the parameters νi and ηi for each 1≤i≤N can be calculated as νi=μi,d2/σi,d2 and ηi=σi,d2/μi,d.Concerning the observation parameters, once the modified forward and backward variables, accounting for the state duration, are defined as in (22) and (28), the reestimation formulas are the same as for Hidden Markov Models [13].In particular, for continuous observations, the parameters of the Gaussians’ mixture defined in (9) are reestimated by firstly defining the probability of being in state Sj at time t with the probability of the observation vector xt evaluated by the kth mixture component, as
(43)γtj,k=αt(j)βt(j)∑j=1N‍αt(j)βt(j)·cjkNxt,μjk,Ujk∑m=1M‍cjmNxt,μjm,Ujm.
By using the former quantity, the parameters cjk, μjk, and Ujk are reestimated through the following formulas:
(44)c¯jk=∑t=1T‍γt(j,k)∑t=1T‍∑m=1M‍γt(j,m),μ¯jk=∑t=1T‍γtj,k·xt∑t=1T‍γt(j,k),U¯jk=∑t=1T‍γt(j,k)·(xt-μjk)xt-μjkT∑t=1T‍γt(j,k),
where superscript T denotes vector transpose.For discrete observations, the reestimation formula for the observation matrixbj(l) is
(45)b¯jl=∑t=1,withxt=XlT‍γtj∑t=1T‍γtj,
where the quantity γt(j), which takes into account the duration dependent forward variable αt(j), is calculated through (17).The reader is referred to Rabiner’s work [13] for the interpretation on the observation parameters reestimation formulas.


## 2.2.1. The Forward-Backward Algorithm

Given a generic sequence of observationsx=x1x2⋯xT, the goal is to calculate the model likelihood, that is, P(x∣λ). This quantity is useful for the training procedure, where the parameters that locally maximize the model likelihood are chosen, as well as for classification problems. The latter is the case in which the observation sequence x has to be mapped to one of a finite set of C classes, represented by a set of HSMM parameters L={λ1,…,λC}. The class of x is chosen such that λ(x)=argmax⁡λ∈LP(X∣λ).To calculate the model likelihood, we first define the forward variable at each timet as
(12)αt(i)=Px1x2⋯xt,st=Si∣λ1≤i≤N.Contrarily to HMMs, for HSMMs the state duration must be taken into account in the the forward variable calculation. Consequently, Yu [26] proposed the following inductive formula:
(13)αtj,d=∑d′=1D‍∑i=1N‍αt-d′i,d′aij0pjjd′∏k=t-d+1tbjxkiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤j≤N,1≤t≤T
that is, the sum of the probabilities of being in the current state Sj (at time t) for the past d′ time units, (with 1≤d′≤D and D the maximum allowed duration for each state), coming from all the possible previous states Si, 1≤i≤N, and i≠j.The disadvantage of the above formulation is that, as discussed in Introduction, the specification of the maximum durationD represents a limitation to the modeling generalization. Moreover, from (13), it is clear that the computation and memory complexities drastically increase with D, which can be very large in many applications, in particular for online failure prediction.To alleviate this problem, Azimi et al. [30–32] introduced a new forward algorithm for HSMMs that, by keeping track of the estimated average state duration at each iteration, has a computational complexity comparable to the forward algorithm for HMMs [13]. However, the average state duration represents an approximation. Consequently, the forward algorithm of Azimi, compared with (13), pays the price of a lower precision in favor of a (indispensable) better computational efficiency.To calculate the forward variableαt(j) using Azimi’s approach, the duration-dependent transition matrix, defined in (8), is taken in consideration in the induction formula of (13), which becomes [30]
(14)αtj=∑i=1N‍αt-1iaijd^t-1bjxt.
To calculate the above formula, the average state duration of (2) must be estimated, for each time t, by means of the variable d^t=[d^t(i)], defined as
(15)d^ti=Edti∣x1x2⋯xt,st=Si,λ1≤i≤N,
where E denotes the expected value. To calculate the above quantity, Azimi et al. [30–32] use the following formula:
(16)d^t=γt-1⊙d^t-1+1,
where ⊙ represents the element by element product between two matrices/vectors and the vector γt=γt(i) (the probability of being in state Si at time t given the observation sequence and the model parameters) with dimensions N×1 is calculated in terms of αt(i) as
(17)γti=Pst=Si∣x1x2⋯xt,λ=αti∑j=1N‍αtjiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.Equation (16) is based on the following induction formula [30–32] that rules the dynamics of the duration vector when the system’s state is known:
(18)dt(i)=st-1(i)·dt-1(i)+1,
where, for each t, st(i) is 1 if st=Si, 0 otherwise.A simple example shows that (18) is incorrect: assuming an HSMM with three states and considering the state sequence (S1,S1,S2,…) the correct sequence of the duration vector is d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,1,1]T, where the superscript T denotes vector transpose. If we apply (18), we obtain d1=[1,1,1]T, d2=[2,1,1]T, and d3=[1,2,1]T, which is in contradiction with the definition of the state duration vector given in (2).To calculate the average state duration variabled^t(i) we propose a new induction formula that estimates, for each time t, the time spent in the ith state prior to t as
(19)d^t(i)=Pst-1=Si∣st=Si,x1,…,xt·d^t-1i+1(20)=aiid^t-1·αt-1(i)·bixtαti·d^t-1i+1,iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i≤N.
The derivation of (20) is given in Appendix. The intuition behind (19) is that the current average duration is the previous average duration plus one, weighted with the “amount” of the current state that was already in state Si in the previous step.Using the proposed (20), the forward algorithm can be specified as follows:(1)
initialization, with1≤i≤N:
(21)α1i=πibix1,d^1(i)=1,Ad^1=Pd^1+I-Pd^1A0,whereP(d^i) is estimated using (6);(2)
induction, with1≤j≤N and 1≤t≤T-1:
(22)αt+1j=∑i=1N‍αtiaijd^tbjxt+1,(23)d^t+1i=aiid^t·αti·bixt+1αt+1i·d^ti+1,(24)Ad^t+1=Pd^t+1+I-Pd^t+1A0,whereaij(d^t) are the coefficients of the matrix Ad^t;(3)
termination:(25)Px∣λ=∑i=1N‍αTi.Similar considerations as the forward procedure can be made for the backward algorithm, which is implemented by defining the variableβt(i) as
(26)βt(i)=Pxt+1xt+2⋯xT∣st=Si,λ1≤i≤N.Having estimated the dynamic transition matrixAd^t for each 1≤t≤T using (24), the backward variable can be calculated inductively as follows.(1)
Initialization:(27)βT(i)=1,1≤i≤N.(2)
Induction:(28)βti=∑j=1N‍aijd^tbjxt+1βt+1j,it=T-1,T-2,…,1,1≤i≤N.Although the variableβt(i) is not necessary for the calculation of the model likelihood, it will be useful in the parameter reestimation procedure, as it will be explained in Section 2.2.3.


## 2.2.2. The Viterbi Algorithm

The Viterbi algorithm [46, 47] (also known asdecoding) allows determining the best state sequence corresponding to a given observation sequence.Formally, given a sequence of observationx=x1x2⋯xT, the best state sequence S*=s1*s2*⋯sT* corresponding to x is calculated by defining the variable δt(i) as
(29)δt(i)=max⁡s1,s2,…,st-1Ps1s2,…,st=Si,x1x2⋯xt∣λ.The procedure to recursively calculate the variableδt(i) and to retrieve the target state sequence (i.e., the arguments which maximize the δt(i)’s) for the proposed HSMM is a straightforward extension of the Viterbi algorithm for HMMs [13]. The only change is the usage, in the recursive calculation of δt(i), of the dynamic transition matrix Ad^t=[aij(d^t)], calculated through (24). The Viterbi algorithm for the introduced parametric HSMMs can be summarized as follows:(1)
initialization, with1≤i≤N:
(30)δ1i=πibix1,ψ1i=0,(2)
recursion, with1≤j≤N and 2≤t≤T:
(31)δtj=max⁡1≤i≤Nδt-1iaijd^tbjxt,(32)ψtj=argmax1≤i≤N⁡δt-1(i)aij(d^t),(3)
termination:(33)P*=max⁡1≤i≤NδT(i),(34)sT*=argmax1≤i≤N⁡δT(i),where we keep track of the argument maximizing (31) using the vector ψt, which, tracked back, gives the desired best state sequence:
(35)st*=ψt+1st+1*,t=T-1,T-2,…,1.


## 2.2.3. The Training Algorithm

The training algorithm consists of estimating the model parameters from the observation data. As discussed in Section2.1, a parametric HSMM is defined by λ=(A0,Θ,C,μ,U,π) if the observations are continuous, or λ=(A0,Θ,B,π) if the observations are discrete. Given a generic observation sequence x=x1x2⋯xT, referred to astraining set in the following, the training procedure consists of finding the model parameter set λ* which locally maximizes the model likelihood Px∣λ.We use the modified Baum-Welch algorithm of Azimi et al. [30–32]. However, in our implementation we do not make assumption on the density function used to model the state duration, and we consider both continuous and discrete observations.Being a variant of the more general Expectation-Maximization (EM) algorithm, Baum-Welch is an iterative procedure which consists of two steps: (i) the expectation step, in which the forward and backward variables are calculated and the model likelihood is estimated and (ii) the maximization step, in which the model parameters are updated and used in the next iteration. This process usually starts from a random guess of the model parametersλ0 and it is iterated until the likelihood function does not improve between two consecutive iterations.Similarly to HMMs, the reestimation formulas are derived by firstly introducing the variableξt(i,j), which represents the probability of being in state Si at time t, and in state Sj at time t+1, given the model and the observation sequence, as
(36)ξti,j=Pst=Si,st+1=Sj∣x,λ.
However, in the HSMM case, the variable ξt(i,j) considers the duration estimation performed in the forward algorithm (see Equation (24)). Formulated in terms of the forward and backward variables, it is given by
(37)ξti,j=Pst=Si,st+1=Sj∣x,λ=Pst=Si,st+1=Sj,x∣λPx∣λ=αtiaijd^tbjxt+1βt+1jPx∣λ=αtiaijd^tbjxt+1βt+1j∑i=1N‍∑j=1N‍αtiaijd^tbjxt+1βt+1j.Fromξt(i,j) we can derive the quantity γt(i) (already defined in (17)) representing the probability of being in state Si at time t given the observation sequence and the model parameters:
(38)γt(i)=∑j=1N‍ξt(i,j).Finally, the the reestimation formulas for the parametersπ and A0 are given by
(39)π¯i=γ1i,(40)a¯ij0=∑t=1T-1‍ξt(i,j)⊙G∑j=1N‍∑t=1T-1‍ξt(i,j)⊙G,
where G=[gij] is a square matrix of dimensions N×N where gij=0 for i=j and gij=1 for i≠j, ⊙ represents the element by element product between two matrices, ∑t=1T-1‍γt(i) is the expected number of transitions from state Si, and ∑t=1T-1‍ξt(i,j) is the expected number of transitions from state Si to state Sj.Equation (39) represents the expected number of times that the model starts in state Si, while (40) represents the expected number of transitions from state Si to state Sj with i≠j over the total expected number of transitions from state Si to any other state different from Si.For the matrixA¯0, being normalized, the stochastic constraints are satisfied at each iteration, that is, ∑j=1N‍a¯ij0=1 for each 1≤i≤N, while the estimation of the prior probability π¯i inherently sums up to 1 at each iteration, since it represents the expected frequency in state Si at time t=1 for each 1≤i≤N.With respect to the reestimation of the state duration parametersΘ, firstly, we estimate the mean μi,d and the variance σi,d2 of the ith state duration for each 1≤i≤N, from the forward and backward variables and the estimation of the state duration variable
(41)μi,d=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1jd^ti∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j,(42)σi,d2=∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j∑j=1,j≠iN‍·d^t(i)-μi,d2·∑t=1T-1‍αti∑j=1,j≠iN‍aijd^tibjxt+1βt+1j-1,
where (41) can be interpreted as the probability of transition from state Si to Sj with i≠j at time t weighted by the duration of state Si at t, giving the desired expected value, while in (42) the same quantity is weighted by the squared distance of the duration at time t from its mean, giving the estimation of the variance.Then, the parameters of the desired duration distribution can be estimated fromμi,d and σi,d2. For example, if a Gamma distribution with shape parameter ν and scale parameter η is chosen to model the state duration, the parameters νi and ηi for each 1≤i≤N can be calculated as νi=μi,d2/σi,d2 and ηi=σi,d2/μi,d.Concerning the observation parameters, once the modified forward and backward variables, accounting for the state duration, are defined as in (22) and (28), the reestimation formulas are the same as for Hidden Markov Models [13].In particular, for continuous observations, the parameters of the Gaussians’ mixture defined in (9) are reestimated by firstly defining the probability of being in state Sj at time t with the probability of the observation vector xt evaluated by the kth mixture component, as
(43)γtj,k=αt(j)βt(j)∑j=1N‍αt(j)βt(j)·cjkNxt,μjk,Ujk∑m=1M‍cjmNxt,μjm,Ujm.
By using the former quantity, the parameters cjk, μjk, and Ujk are reestimated through the following formulas:
(44)c¯jk=∑t=1T‍γt(j,k)∑t=1T‍∑m=1M‍γt(j,m),μ¯jk=∑t=1T‍γtj,k·xt∑t=1T‍γt(j,k),U¯jk=∑t=1T‍γt(j,k)·(xt-μjk)xt-μjkT∑t=1T‍γt(j,k),
where superscript T denotes vector transpose.For discrete observations, the reestimation formula for the observation matrixbj(l) is
(45)b¯jl=∑t=1,withxt=XlT‍γtj∑t=1T‍γtj,
where the quantity γt(j), which takes into account the duration dependent forward variable αt(j), is calculated through (17).The reader is referred to Rabiner’s work [13] for the interpretation on the observation parameters reestimation formulas.


## 3. AIC-Based Model Selection

In the framework of the proposed parametric HSMMs, the model selection procedure aims to select the optimal number of hidden statesN, the right duration distribution family, and, in the case of mixture observation modeling, the number of Gaussian mixtures M to be used. In this work, we make use of the Akaike Information Criterion (AIC). Indeed, it has been seen that in case of complex models and in presence of a limited number of training observations, AIC represents a satisfactory methodology for model selection, outperforming other approaches like Bayesian Information Criterion.In general, information criteria are represented as a two-term structure. They account for a compromise between a measure of model fitness, which is based on the likelihood of the model, and a penalty term which takes into account the model complexity. Usually, the model complexity is measured in terms of the number of parameters that have to be estimated and in terms of the number of observations.The Akaike Information Criterion is an estimate of the asymptotic value of the expected distance between the unknown true likelihood function of the data and the fitted likelihood function of the model. In particular, the AIC can be expressed as(46)AIC=-log⁡L(λ^)+pT,
where L(λ^) is the likelihood of the model with the estimated parameters λ^, as defined in (25), p is the number of model parameters, and T is the length of the observed sequence. The best model is the one minimizing equation (46).Concerningp, the number of parameters to be estimated for a parametric HSMM with N states is p=ph+po, where ph are the parameters of the hidden states layer, while po are those of the observation layer.In particularph=(N-1)+(N-1)·N+z·N where(i)
N-1 accounts for the prior probabilities π;(ii)
(N-1)·N accounts for the nonrecurrent transition matrix A0;(iii)
z·N accounts for the duration probability, being z the number of parameters θ of the duration distribution.Concerningpo, a distinction must be made between discrete and continuous observations:(i)
in the case of discrete observations withL possible observable values, po=(L-1)·N, which accounts for the elements of the observation matrix B;(ii)
if the observations are continuous and a multivariate mixture ofM Gaussians with O variates is used as observation model, po=O·N·M+O·O·N·M+(M-1)·N where each term accounts, respectively, for the mean vector μ, the covariance matrix U, and the mixture coefficients C.


## 4. Remaining Useful Lifetime Estimation

One of the most important advantages of the time modeling of HSMMs is the possibility to effectively face the prediction problem. The knowledge of the state duration distributions allows the estimation of the remaining time in a certain state and, in general, the prediction of the expected timeD~ before entering in a determinate state.As already mentioned, an interesting application of the prediction problem is the Remaining Useful Lifetime (RUL) estimation of industrial equipments. Indeed, if each state of an HSMM is mapped to a different condition of an industrial machine and if the stateSk that represents thefailure condition is identified, at each moment, the RUL can be defined as the expected time D~ to reach the failure state Sk. If we assume that the time to failure is a random variable D following a determinate probability density, we define the RUL at the current time t as
(47)RULt=D~=ED:st+D~=Sk,st+D~-1=Siiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii1≤i,k≤N,i≠k,
where E denotes the expected value.Having fixed thefailure state, the estimation of the RUL is performed, in two steps, every time a new observation is acquired (online):(1)
estimation of the current state;(2)
projection of the future state transitions until the failure state is reached and estimation of the expected sojourn time.The estimation of the current state is performed via the Viterbi path, that is, the variableδt=[δt(i)]1≤i≤N defined in (29). To correctly model the uncertainty of the current state estimation, we use the normalized variable δ¯t(i) obtained as
(48)δ¯ti=max⁡s1,s2,…,st-1Pst=Si∣s1s2⋯st-1,x1x2⋯xt,λ=δt(i)∑j=1N‍δt(j),1≤i≤N
that is, an estimate of the probability of being in state Si at time t.Together with the normalized variableδ¯t(i), the maximum a posteriori estimate of the current state s^t* is taken into account according to (34). If s^t* coincides with the failure state, the desired event is detected by the model and the time to this event is obviously zero. Otherwise, an estimation of the average remaining time in the current state d~avg(s^t*) is calculated as
(49)d~avgs^t*=∑i=1N‍μdi-d^ti⊙δ¯ti,
where with μdi we denote the expected value of the duration variable in state Si according to the duration distribution specified by the parameters θi. Equation (49) thus estimates the remaining time in the current state by subtracting the estimated states duration, d^t(i), at time t from the expected sojourn time of state Si, and weighting the result using the uncertainty about the current state, δ¯t(i), and, finally, by summing up all the contributions from each state.In addition to the average remaining time, a lower and an upper bound value can be calculated based on the standard deviation,σdi, of the duration distribution for state Si:
(50)d~lows^t*=∑i=1N‍μdi-σdi-d^ti⊙δ¯ti,(51)d~ups^t*=∑i=1N‍μdi+σdi-d^ti⊙δ¯ti.Once the remaining time in the current state is estimated, the probability of the next state is calculated by multiplying the transpose of the nonrecurrent transition matrix by the current state probability estimation as follows:(52)δ¯next=δ¯t+d~i1≤i≤N=A0T·δ¯t,
while the maximum a posteriori estimate of the next state s^next* is calculated as
(53)s^next*=s^t+d~*=argmax1≤i≤N⁡δ¯t+d~(i).Again, ifs^t+d~* coincides with the failure state, the failure will happen after the remaining time at the current state is over and the average estimation of the failure time is D~avg=d~avg(s^t*) calculated at the previous step, with the bound values D~low=d~low(s^t*) and D~up=d~up(s^t*). Otherwise the estimation of the sojourn time of the next state is calculated as follows:
(54)d~avgs^t+d~*=∑i=1N‍μdi⊙δ¯t+d~i,(55)d~lows^t+d~*=∑i=1N‍μdi-σdi⊙δ¯t+d~i,(56)d~ups^t+d~*=∑i=1N‍μdi+σdi⊙δ¯t+d~i.This procedure is repeated until the failure state is encountered in the prediction of the next state. The calculation of the RUL is then simply obtained by summing all the estimated remaining time in each intermediate state before encountering the failure state:(57)D~avg=∑‍d~avg,(58)D~low=∑‍d~low,(59)D~up=∑‍d~up.Finally, Algorithm1 details the above described RUL estimation procedure.Algorithm 1:Remaining Useful Lifetime estimation (Pseudo-Code).
(1)function  RulEstimation(xt, Sk)       ⊳ xt: The last observation acquired
(2)                         ⊳ Sk: The failure state
(3)Initialization:
(4)   D~avg←0
(5)   D~low←0
(6)   D~up←0
(7)Current state estimation:
(8)Calculate δ¯t                    ⊳ Using (48)
(9)Calculate s^t*                    ⊳ Using (34)
(10)Calculate d^t                    ⊳ Using (20)
(11)S^←s^t*
(12)Loop:
(13)while  S^≠Sk  do
(14)Calculate d~avg                 ⊳ Using (49) or (54)
(15) Calculate d~low                 ⊳ Using (50) or (55)
(16) Calculate d~up                  ⊳ Using (51) or (56)
(17) D~avg←D~avg+d~avg
(18) D~low←D~low+d~low
(19) D~up←D~up+d~up
(20)Calculate δ¯next                    ⊳ Using (52)
(21)Calculate s^next*                   ⊳ Using (53)
(22)S^←s^next*
end
(23) return  D~avg, D~low, D~up


## 5. Experimental Results

To demonstrate the effectiveness of the proposed HSMM models, we make use of a series of experiments, performed both on simulated and real data.The simulated data were generated by considering a left-right HSMM and adapting the parameters of the artificial example reported in the work of Lee et al. [15]. The real case data are monitoring data related to the entire operational life of bearings, made available for the IEEE PHM 2012 data challenge (http://www.femto-st.fr/en/Research-departments/AS2M/Research-groups/PHM/IEEE-PHM-2012-Data-challenge.php).


### 5.1. Simulated Experiment

Data have been generated with the idea of simulating the behavior of an industrial machine that, during its functioning, experiences several degradation modalities until a failure state is reached at the end of its lifetime. The generated data are used to test the performances of our methodology for (i) automatic model selection, (ii) online condition monitoring, and (iii) Remaining Useful Lifetime estimation, considering both continuous and discrete observations.


#### 5.1.1. Data Generation

The industrial machine subject of these experiments has been modeled as a left-right parametric HSMM, withN=5 states, having state S5 as absorbing (failure) state. The choice of a left-right setting has been made for simplicity reasons, since the primary goal of this work is to demonstrate that the proposed model specification coupled with the Akaike Information Criterion is effective to solve automatic model selection, online condition monitoring, and prediction problems. At this purpose, we divided the experiments in two cases according to the nature of the observation, being continuous or discrete.For each of the continuous and the discrete cases, three data sets have been generated by considering the following duration models: Gaussian, Gamma, and Weibull densities. For each of the three data sets, 30 series of data are used as training set and 10 series as testing set. Each time series containsT=650 observations. The parameters used to generate the data are taken from the work of Lee et al. [15] and are adapted to obtain an equivalent left-right parametric HSMM as follows:
(60)π=10000,A0=0100000100000100000100001,ΘN=θ1=100;20,θ2=90;15,θ3=100;20,θ4=80;25,θ5=200;1,ΘG=θ1=500;0.2,θ2=540;0.1667,θ3=500;0.2,θ4=256;0.3125,θ5=800;0.005,ΘW=θ1=102;28,θ2=92;29,θ3=102;28,θ4=82;20,θ5=200;256,
where ΘN, ΘG, and ΘW are the different distribution parameters used for the Gaussian, Gamma, and Weibull duration models, respectively. More precisely, they represent the values of the mean μd and the the variance σd2 of the Gaussian distribution, the shape νd and the scale ηd of the Gamma distribution, and the scale ad and the shape bd of the Weibull distribution. It must be noticed that, as explained in Section 2.1, for state S5, being the absorbing state, the duration parameters θ5 have no influence on the data, since, once the state S5 is reached, the system will remain there forever.Concerning the continuous observation modeling, a bivariate Gaussian distribution has been used with the following parameters [15]:
(61)μ1=2020,μ2=2035,μ3=3535,μ5=2828,U1=200020,U2=150015,U3=15-2-215,U4=5005,U5=103310
while for the discrete case, L=7 distinct observation symbols have been taken into consideration with the following observation probability distribution
(62)B=0.80.2000000.10.80.1000000.10.80.1000000.10.70.10.100000.20.60.10.1.
An example of simulated data both for the continuous and the discrete cases is shown in Figure 2, where a Gaussian duration model has been used.The data generated with the parameter described in Section5.1.1, both for the continuous and the discrete case.
(a)
Example of simulated data for the continuous case
(b)
Example of simulated data for the discrete case


#### 5.1.2. Training and Model Selection

The goal of this experimental phase is to test the effectiveness of the AIC in solving the automatic model selection. For this purpose, the training sets of the 6 data sets (continuous/discrte observation with Gaussian, Gamma, and Weibull duration models) have been taken individually and, for each one of them, a series of learning procedure has been run, each one with a variable HSMM structure. In particular we took into account all the combinations of the duration distribution families (Gaussian, Gamma, and Weibull), an increasing number of states,N, from 2 to 8 and, for the continuous observation cases, an increasing number of Gaussian mixtures, M, in the observation distribution from 1 to 4.As accurate parameter initialization is crucial for obtaining a good model fitting [14], a series of 40 learning procedures corresponding to 40 random initializations of the initial parameters λ0 have been executed for each of the considered HSMM structures. For each model structure, the AIC value, as defined in (46), has been evaluated. The final trained set of parameters λ* corresponding to the minimum AIC value has been retained, resulting in 7 HSMMs with a number of states from 2 to 8.The obtained results are shown in Figure3, for both, the continuous and discrete observation data. As it can be noticed, for all the 6 test cases of Figure 3 the AIC values do not improve much for a number of states higher than 5, meaning that adding more states does not add considerable information to the HSMM modeling power. Hence 5 states can be considered as an optimal number of states. Moreover, as shown in the zoomed sections of Figure 3, for the HSMMs with 5 states, the minimum AIC values are obtained for the duration distributions corresponding to the ones used to generate the data. As a consequence AIC can be considered as an effective approach to perform model selection for HSMM, as well as selecting the appropriate parametric distribution family for the state duration modeling.Akaike Information Criterion (AIC) applied to continuous and discrete observations data. AIC is effective for automatic model selection, since its minimum value provides the same number of states and duration model used to generate the data.
(a)
AIC values for continuous data and Gaussian duration distribution
(b)
AIC values for continuous data and Gamma duration distribution
(c)
AIC values for continuous data and Weibull duration distribution
(d)
AIC values for discrete data and Gaussian duration distribution
(e)
AIC values for discrete data and Gamma duration distribution
(f)
AIC values for discrete data and Weibull duration distribution


#### 5.1.3. Condition Monitoring

The optimal parametersλ* obtained in the previous phase have been used to perform the online condition monitoring experiment on the 10 testing cases for all the 6 considered HSMM configurations. In this experiment, we simulate online monitoring by considering all the testing observations up to the current time, that is, {x1x2⋯xt}. Each time a new data point is acquired, the Viterbi algorithm is used to estimate the current state st*=argmax1≤i≤Nδt(i) as specified in (34).An example of execution of the condition monitoring experiment is shown in Figure4, for both, continuous and discrete observations, respectively. In Figure 4(a) the state duration has been modeled with a Gamma distribution, while in Figure 4(b) with a Gaussian distribution. In Figures 4(a) and 4(b), the first display represents the true state of the HSMM and the second display represents the estimated state from the Viterbi algorithm, while the third display represents the observed time series.Condition monitoring using the Viterbi path. HSMMs can be effective to solve condition monitoring problems in time-dependent applications due to their high accuracy in hidden state recognition.
(a)
State estimation with Viterbi path for continuous data and Gamma duration distribution
(b)
State estimation with Viterbi path for discrete data and Gaussian duration distributionKnowing the true state sequence we calculated the accuracy, defined as the percentage of correctly estimated states over the total length of the state sequence, for each of the testing cases. The results are summarized in Table1(a) for the continuous observations and in Table 1(b) for the discrete observations. The high percentage of correct classified states shows that HSMMs can be effectively used to solve condition monitoring problems for applications in which the system shows a strong time and state duration dependency.Table 1
State recognition accuracy.
(a)
Continuous observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
99.4%
98.5%
99.2%
2
99.7%
98.6%
99.5%
3
99.4%
99.2%
99.7%
4
98.9%
98.9%
99.7%
5
98.2%
98.9%
100%
6
99.1%
98.8%
99.7%
7
98.5%
99.4%
99.7%
8
99.2%
99.1%
99.5%
9
99.2%
98.6%
99.7%
10
99.2%
99.1%
99.5%
Average
99.1%
98.9%
99.6%
(b)
Discrete observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
97.4%
96.7%
97.4%
2
97.2%
97.6%
96.5%
3
99.4%
95.8%
96.6%
4
98.2%
95.3%
97.7%
5
99.1%
97.4%
97.5%
6
97.8%
97.7%
97.8%
7
95.8%
97.2%
96.6%
8
97.7%
96.4%
97.2%
9
98.9%
97.2%
98.5%
10
99.2%
95.6%
96.9%
Average
98.1%
96.7%
97.3%


#### 5.1.4. Remaining Useful Lifetime Estimation

In this experimental phase we considered the stateS5 as the failure state and the trained parameters λ* of Section 5.1.2 for each HSMM configuration. As the online RUL estimation procedure is intended to be used in real time, we simulated condition monitoring experiment where we progressively consider the observations {x1x2⋯xt} up to time t. When a new observation is acquired, after the current state probability δ¯t(i) is estimated (Equation (48)), the calculation of the average, upper, and lower RUL ((57), (58), and (59)) is performed.Examples of RUL estimation are illustrated in Figure5. In particular Figure 5(a) represents the case of continuous observations and duration modeled by a Weibull distribution, while Figure 5(b) shows the case of discrete observations and duration modeled by a Gamma distribution. From the figures one can notice that the average, as well as the lower and the upper bound estimations, converges to the real RUL as the failure time approaches. Moreover, as expected, the uncertainty about the estimation decreases with the time, since predictions performed at an early stage are more imprecise. As a consequence, the upper and the lower bound become more narrow as the failure state approaches, and the estimation becomes more precise until it converges to the actual RUL value, with the prediction error tending to zero at the end of the evaluation.HSMMs effectively solve RUL estimation problems. The prediction converges to the actual RUL value and its uncertainty decreases as the real failure time approaches.
(a)
Remaining Useful Lifetime estimation for continuous data and Weibull duration distribution
(b)
Remaining Useful Lifetime estimation for discrete data and Gamma duration distributionTo quantitatively estimate the performance of our methodology for the RUL estimation, we considered at each timet the absolute prediction error (APE) between the real RUL and the predicted value defined as
(63)APEt=RULrealt-RULt,
where RULreal(t) is the (known) value of the RUL at time t, while RUL(t) is RUL predicted by the model. To evaluate the overall performance of our methodology, we considered the average absolute prediction error of the RUL estimation, defined as
(64)APE¯=∑t=1T‍APE(t)T,
where T is the length of the testing signal. APE¯ being a prediction error average, values of (64) close to zero correspond to good predictive performances.The results for each of the 10 testing cases and the different HSMM configurations are reported in Tables2(a) and 2(b), for the continuous and the discrete observation cases, respectively. As it can be noticed, the online Remaining Useful Lifetime estimation and in general the online prediction of the time to a certain event can be effectively faced with HSMMs, which achieve a reliable estimation power with a small prediction error.Table 2
Average absolute prediction error (APE¯) of the RUL estimation using the proposed state duration estimator of (20).
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
5.1
17.0
6.7
14.0
29.0
0.91
4.5
17.0
8.1
2
7.6
19.0
5.0
6.1
21.0
8.5
6.6
19.0
6.1
3
7.7
5.4
19.0
2.9
12.0
17.0
16.0
29.0
3.0
4
9.0
21.0
2.9
7.5
22.0
6.8
6.0
19.0
6.7
5
7.3
19.0
4.7
2.2
14.0
14.0
3.9
17.0
8.7
6
6.5
18.0
5.6
5.1
18.0
10.0
14.0
27.0
2.7
7
4.7
16.0
7.5
4.8
17.0
11.0
1.2
13.0
12.0
8
10.0
22.0
2.9
5.2
18.0
10.0
9.2
22.0
3.9
9
3.1
9.2
14.0
2.0
16.0
13.0
8.2
21.0
4.9
10
6.4
18.0
5.6
7.5
22.0
6.9
3.3
12.0
13.0
Average
6.8
17.0
7.4
5.7
19.0
9.9
7.3
20.0
7.0
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
2.1
11.0
14.0
3.1
8.8
14.0
2.4
12.0
13.0
2
2.1
11.0
13.0
11.0
22.0
3.3
19.0
32.0
7.1
3
5.1
17.0
7.6
6.6
18.0
5.1
2.3
14.0
11.0
4
5.9
6.5
18.0
5.2
17.0
6.7
4.2
16.0
9.0
5
3.2
14.0
10.0
8.3
19.0
3.4
12.0
24.0
2.9
6
12.0
24.0
2.7
6.2
18.0
5.2
4.1
8.4
16.0
7
2.9
15.0
9.7
9.3
21.0
2.3
19.0
31.0
6.6
8
15.0
27.0
7.0
7.4
18.0
4.3
4.3
17.0
9.4
9
5.9
18.0
7.7
11.0
23.0
5.5
3.9
16.0
8.8
10
3.5
11.0
14.0
5.5
6.0
16.0
5.2
17.0
7.1
Average
5.7
15.0
10.0
7.4
17.0
6.6
7.7
19.0
9.0Finally, we tested our RUL estimation methodology using the state duration estimation of (16) introduced by Azimi et al. [30–32]. The results are shown in Tables 3(a) and 3(b), in which, respectively, the prediction errors obtained for continuous and discrete observations are reported.Table 3
Average absolute prediction error (APE¯) of the RUL estimation using the state duration estimator of (16) introduced by Azimi et al. [30–32].
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
57.8
51.0
66.8
26.2
9.7
52.7
25.9
28.4
64.6
2
50.2
44.4
57.7
21.3
17.0
46.9
29.0
19.2
70.8
3
50.3
44.7
57.3
27.1
8.7
56.5
34.5
13.9
73.4
4
51.8
46.0
60.4
21.3
14.3
45.9
34.9
17.1
78.7
5
59.4
53.7
66.2
29.0
9.5
55.4
33.4
15.6
74.9
6
58.0
51.7
67.1
25.8
8.3
54.1
23.1
25.8
66.5
7
59.4
53.6
66.9
18.2
12.5
47.7
36.0
17.1
74.4
8
63.4
55.6
72.3
19.4
15.7
44.1
34.8
17.8
77.0
9
49.1
43.5
57.0
14.5
17.1
43.2
25.1
26.7
67.0
10
54.4
48.4
62.8
23.2
7.9
52.7
24.1
24.5
67.4
Average
55.4
49.3
63.5
22.6
12.1
49.9
30.1
20.6
71.5
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
51.4
41.0
62.4
42.4
31.8
53.0
32.6
26.4
73.6
2
49.6
39.9
60.4
59.5
48.3
70.8
31.3
27.6
69.3
3
50.2
38.6
62.3
46.5
35.7
57.4
32.4
25.7
70.2
4
42.2
31.5
53.8
50.1
40.5
60.6
23.7
36.1
60.3
5
44.3
33.9
55.8
47.8
37.4
59.1
36.0
25.6
76.5
6
52.2
43.2
62.7
55.2
44.3
66.9
27.2
31.6
64.3
7
55.0
43.9
66.8
56.0
45.7
67.0
34.7
23.2
74.4
8
50.3
39.0
62.0
60.4
50.5
71.0
35.1
26.4
72.4
9
55.5
47.4
64.0
48.0
37.2
59.5
31.8
22.2
73.6
10
49.0
38.2
60.7
52.1
41.2
63.1
29.4
28.9
68.7
Average
50.0
39.7
61.1
51.8
41.3
62.9
31.4
27.4
70.4Comparing Table2 and Table 3, one can notice that the proposed RUL method outperforms the one of Azimi. This is mainly due to the proposed average state duration of (20), compared to the one of Azimi given by (16).


### 5.2. Real Data

In this section we apply the proposed HSMM-based approach for RUL estimation to a real case study using bearing monitoring data recorded during experiments carried out on the Pronostia experimental platform and made available for the IEEE Prognostics and Health Management (PHM) 2012 challenge [49]. The data correspond to normally degraded bearings, leading to cases which closely correspond to the industrial reality.The choice of testing the proposed methodology on bearings derives from two facts: (i) bearings are the most critical components related to failures of rotating machines [50] and (ii) their monotonically increasing degradation pattern justifies the usage of left-right HSMM models.


#### 5.2.1. Data Description

Pronostia is an experimental platform designed and realized at the Automatic Control and Micro-Mechatronic Systems (AS2M) Department of Franche-Comté Electronics, Mechanics, Thermal Processing, Optics-Science and Technology (FEMTO-ST) Institute (http://www.femto-st.fr/) (Besançon, France), with the aim of collecting real data related to accelerated degradation of bearings. Such data are used to validate methods for bearing condition assessment, diagnostic and prognostic [19, 51–59].The Pronostia platform allows to perform run-to-failure tests under constant or variable operating conditions. The operating conditions are determined by two factors that can be controlled online: (i) the rotation speed and (ii) the radial force load. During each experiment, temperature and vibration monitoring measurements are gathered online, through two type of sensors placed in the bearing housing: a temperature probe and two accelerometers (one on the vertical and one on the horizontal axis).The platform is composed of three main parts: a rotating part, a load profile generation part, and a measurement part, as illustrated in Figure6.Figure 6
Global overview of the Pronostia experimental platform [19].The rotating part is composed of an asynchronous motor which develops a power equal to 250 W, two shafts, and a gearbox, which allows the motor to reach its rated speed of 2830 rpm. The motor’s rotation speed and the direction are set through a human machine interface.The load profiles part issues a radial force on the external ring of the test bearing through a pneumatic jack connected to a lever arm, which indirectly transmits the load through a clamping ring. The goal of the applied radial force is to accelerate the bearing’s degradation process.The measurement part consists of a data acquisition card connected to the monitoring sensors, which provides the user with the measured temperature and vibration data. The vibration measurements are provided in snapshots of 0.1 s collected each 10 seconds at a sampling frequency of 25.6 kHz (2560 samples per each snapshot), while the temperature has been continuously recorded at a sampling frequency of 10 Hz (600 samples collected each minute).Further details on the Pronostia test rig can be found on the data presentation paper [49] and on the web page of the data challenge (http://www.femto-st.fr/en/Research-departments/AS2M/Research-groups/PHM/IEEE-PHM-2012-Data-challenge.php).Regarding the data provided for the PHM 2012 challenge, 3 different operating conditions were considered:(i)
first operating conditions: speed of 1800 rpm and load of 4000 Newton;(ii)
second operating conditions: speed of 1650 rpm and load of 4200 Newton;(iii)
third operating conditions: speed of 1500 rpm and load of 5000 Newton.Under the above operating conditions, a total of 17 accelerated life tests were realized on bearings of type NSK 6804 DD, which can operate at a maximum speed of 13000 rpm and a load limit of 4000 N. The tests were stopped when the amplitude of the vibration signal was higher than 20 g; thus this moment was defined as the bearing failure time. An example of bearing before and after the experiment is shown in Figure7 together with the corresponding vibration signal collected during the whole test.Figure 7
A tested bearing before and after the experiment with its recorded vibration signal [49].Table4 reports how the 17 tested bearings were separated into the three operating conditions. Moreover, the duration of each experiment, being the RUL to be predicted for each bearing, is also given. We performed two sets of experiments by considering, respectively, the bearings relative to the first and the second operating condition (i.e.,Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Table 4
Lifetime duration (in seconds) and operating conditions of the bearings tested in the IEEE PHM 2012 Prognostic Challenge dataset [49].
Condition 1
Condition 2
Condition 3
1800 rpm and 4000 N
1650 rpm and 4200 N
1500 rpm and 5000 N
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing1_1
28030
Bearing2_1
9110
Bearing3_1
5150
Bearing1_2
8710
Bearing2_2
7970
Bearing3_2
16370
Bearing1_3
23750
Bearing2_3
19550
Bearing3_3
4340
Bearing1_4
14280
Bearing2_4
7510
Bearing1_5
24630
Bearing2_5
23110
Bearing1_6
24480
Bearing2_6
7010
Bearing1_7
22590
Bearing2_7
2300As already mentioned, the available data correspond to normally degraded bearings, meaning that the defects were not initially induced and that each degraded bearing contains almost all the types of defects (balls, rings, and cage), resembling faithfully a common real industrial situation. Moreover, no assumption about the type of failure to be occurred is provided with the data and, since the variability in experiment durations is high (from 1 h to 7 h), performing good estimates of the RUL is a difficult task [49].In our experiments we considered, as input to our model, the horizontal channel of the accelerometer. We preprocessed the raw signals by extracting two time-domain features, that is, root mean square (RMS) and kurtosis, within windows of the same length as the given snapshots (L=2560). Let rw(t) be the raw signal of the wth window; for each w we estimate RMS as xwRMS=1/L∑t=1L‍rw2(t) and kurtosis as xwKURT=1/L∑t=1L‍rw(t)-r¯w4/1/L∑t=1L‍rw(t)-r¯w22, where r¯w is the mean of rw. An example of feature extraction forBearing1_1 is shown in Figure 8.Figure 8
Raw vibration data (a) versus RMS and kurtosis features (b) forBearing1_1.To assess the performance of the proposed HSMM, after the model selection procedure, we implemented a leave-one-out cross validation scheme: by considering separately conditions 1 and 2, we performed the online RUL estimation for each of the 7 bearings, using an HSMM trained with the remaining 6 bearing histories. Similarly to the simulated case, we considered the average absolute prediction error, defined in (64), to quantitatively evaluate our method.


#### 5.2.2. Bearings RUL Estimation

We performed our experiments in two steps: firstly we applied model selection in order to determine an optimal model structure, and secondly, we estimated the RUL of the bearings. The full procedure is detailed in the following.(A) HSMM Structure. To determine an appropriate HSMM structure for effectively modeling the considered data, we considered several HSMM structures, characterized by (i) the duration distribution family (being Gaussian, Gamma, or Weibull), (ii) an increasing number of states, N, from 2 to 6, and (iii) an increasing number of Gaussian mixtures, M, in the observation density from 1 to 4. For each model structure, obtained by systematically considering all the combinations of (i) to (iii), we run 120 parameter learnings, corresponding to 120 random initializations, λ0, on the data sets (Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Similar to Section5.1.2, at the end of each learning we evaluated the AIC value (Equation (46)) as reported in Figures 9(a) and 9(b) for conditions 1 and 2, respectively. In both cases, the global minimum AIC value corresponds to an HSMM with N=4 states, a Weibull duration model, and a M=1 Gaussians mixture for the observation density.In both cases the minimum AIC value corresponds to an HSMM withN=4 states, a Weibull duration model, and M=1 mixture in the observation density.
(a)
AIC values for Condition 1
(b)
AIC values for Condition 2(B) RUL Estimation. Using the above obtained optimal HSMM structure, we trained it via a leave-one-out cross validation scheme by using for condition 1, at each iteration,Bearing1_i, 1≤i≤7, as the testing bearing, while the remaining six bearings were used for training. Once the trained parameters λi* were estimated for the ith testing bearing, we progressively collected the observations of the testedBearing1_i to calculate, at each time t, the average, lower, and upper RUL, as specified in (57), (58), and (59), respectively, considering the state S4 as the failure state. The same procedure has been performed for the bearings in condition 2.Examples of RUL estimation forBearing1_7 andBearing2_6 are shown in Figures 10(a) and 10(b), respectively, where the black solid line represents the real RUL which goes to zero as the time goes on. As it can be seen, the average as well as the lower and the upper bound estimations converge to the real RUL as the real failure time approaches and the uncertainty about the estimation decreases with time.By obtaining a low average absolute prediction error, the proposed parametric HSMM is effective for estimating the Remaining Useful Lifetime of bearings.
(a)
RUL estimation for Bearing1_7
(b)
RUL estimation for Bearing2_6Concerning the quantitative estimation of the predictive performances, we report in Table5 the average absolute prediction error of the RUL estimation (see Equation (64)), expressed in seconds. As it can be noticed, the average absolute prediction error of the average RUL is, respectively, 1 hour and 15 minutes for condition 1, and 1 hour and 5 minutes for condition 2, which are good values, considering the high variability in the training set durations and the fact that the performance metric takes into account also the less accurate predictions performed at an early stage of the bearings life. Moreover, for both conditions, the average prediction errors of 5 tests out of 7 are below the average, while the best average error of the mean RUL is only 23 minutes for condition 1 while it further decreases to 14 minutes for condition 2.Table 5
Average absolute prediction error (APE¯) of the RUL estimation, expressed in seconds.
(a)
Condition 1
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing1_1
10571.6
12723.0
9414.6
Bearing1_2
4331.2
3815.6
3821.3
Bearing1_3
2997.0
9730.9
6091.2
Bearing1_4
6336.3
2876.6
14871.9
Bearing1_5
1968.9
7448.4
10411.5
Bearing1_6
4253.0
9896.4
9793.7
Bearing1_7
1388.0
7494.3
10088.1
Average
4549.4
7712.2
9213.2
(b)
Condition 2
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing2_1
2475.9
5006.5
7287.5
Bearing2_2
1647.3
4497.2
8288.6
Bearing2_3
8877.1
9508.3
7962.1
Bearing2_4
1769.8
4248.6
4982.5
Bearing2_5
8663.1
10490.0
10730.0
Bearing2_6
877.1
3504.7
6687.0
Bearing2_7
3012.5
3866.4
6651.9
Average
3903.3
5874.5
7512.8


## 5.1. Simulated Experiment

Data have been generated with the idea of simulating the behavior of an industrial machine that, during its functioning, experiences several degradation modalities until a failure state is reached at the end of its lifetime. The generated data are used to test the performances of our methodology for (i) automatic model selection, (ii) online condition monitoring, and (iii) Remaining Useful Lifetime estimation, considering both continuous and discrete observations.


### 5.1.1. Data Generation

The industrial machine subject of these experiments has been modeled as a left-right parametric HSMM, withN=5 states, having state S5 as absorbing (failure) state. The choice of a left-right setting has been made for simplicity reasons, since the primary goal of this work is to demonstrate that the proposed model specification coupled with the Akaike Information Criterion is effective to solve automatic model selection, online condition monitoring, and prediction problems. At this purpose, we divided the experiments in two cases according to the nature of the observation, being continuous or discrete.For each of the continuous and the discrete cases, three data sets have been generated by considering the following duration models: Gaussian, Gamma, and Weibull densities. For each of the three data sets, 30 series of data are used as training set and 10 series as testing set. Each time series containsT=650 observations. The parameters used to generate the data are taken from the work of Lee et al. [15] and are adapted to obtain an equivalent left-right parametric HSMM as follows:
(60)π=10000,A0=0100000100000100000100001,ΘN=θ1=100;20,θ2=90;15,θ3=100;20,θ4=80;25,θ5=200;1,ΘG=θ1=500;0.2,θ2=540;0.1667,θ3=500;0.2,θ4=256;0.3125,θ5=800;0.005,ΘW=θ1=102;28,θ2=92;29,θ3=102;28,θ4=82;20,θ5=200;256,
where ΘN, ΘG, and ΘW are the different distribution parameters used for the Gaussian, Gamma, and Weibull duration models, respectively. More precisely, they represent the values of the mean μd and the the variance σd2 of the Gaussian distribution, the shape νd and the scale ηd of the Gamma distribution, and the scale ad and the shape bd of the Weibull distribution. It must be noticed that, as explained in Section 2.1, for state S5, being the absorbing state, the duration parameters θ5 have no influence on the data, since, once the state S5 is reached, the system will remain there forever.Concerning the continuous observation modeling, a bivariate Gaussian distribution has been used with the following parameters [15]:
(61)μ1=2020,μ2=2035,μ3=3535,μ5=2828,U1=200020,U2=150015,U3=15-2-215,U4=5005,U5=103310
while for the discrete case, L=7 distinct observation symbols have been taken into consideration with the following observation probability distribution
(62)B=0.80.2000000.10.80.1000000.10.80.1000000.10.70.10.100000.20.60.10.1.
An example of simulated data both for the continuous and the discrete cases is shown in Figure 2, where a Gaussian duration model has been used.The data generated with the parameter described in Section5.1.1, both for the continuous and the discrete case.
(a)
Example of simulated data for the continuous case
(b)
Example of simulated data for the discrete case


### 5.1.2. Training and Model Selection

The goal of this experimental phase is to test the effectiveness of the AIC in solving the automatic model selection. For this purpose, the training sets of the 6 data sets (continuous/discrte observation with Gaussian, Gamma, and Weibull duration models) have been taken individually and, for each one of them, a series of learning procedure has been run, each one with a variable HSMM structure. In particular we took into account all the combinations of the duration distribution families (Gaussian, Gamma, and Weibull), an increasing number of states,N, from 2 to 8 and, for the continuous observation cases, an increasing number of Gaussian mixtures, M, in the observation distribution from 1 to 4.As accurate parameter initialization is crucial for obtaining a good model fitting [14], a series of 40 learning procedures corresponding to 40 random initializations of the initial parameters λ0 have been executed for each of the considered HSMM structures. For each model structure, the AIC value, as defined in (46), has been evaluated. The final trained set of parameters λ* corresponding to the minimum AIC value has been retained, resulting in 7 HSMMs with a number of states from 2 to 8.The obtained results are shown in Figure3, for both, the continuous and discrete observation data. As it can be noticed, for all the 6 test cases of Figure 3 the AIC values do not improve much for a number of states higher than 5, meaning that adding more states does not add considerable information to the HSMM modeling power. Hence 5 states can be considered as an optimal number of states. Moreover, as shown in the zoomed sections of Figure 3, for the HSMMs with 5 states, the minimum AIC values are obtained for the duration distributions corresponding to the ones used to generate the data. As a consequence AIC can be considered as an effective approach to perform model selection for HSMM, as well as selecting the appropriate parametric distribution family for the state duration modeling.Akaike Information Criterion (AIC) applied to continuous and discrete observations data. AIC is effective for automatic model selection, since its minimum value provides the same number of states and duration model used to generate the data.
(a)
AIC values for continuous data and Gaussian duration distribution
(b)
AIC values for continuous data and Gamma duration distribution
(c)
AIC values for continuous data and Weibull duration distribution
(d)
AIC values for discrete data and Gaussian duration distribution
(e)
AIC values for discrete data and Gamma duration distribution
(f)
AIC values for discrete data and Weibull duration distribution


### 5.1.3. Condition Monitoring

The optimal parametersλ* obtained in the previous phase have been used to perform the online condition monitoring experiment on the 10 testing cases for all the 6 considered HSMM configurations. In this experiment, we simulate online monitoring by considering all the testing observations up to the current time, that is, {x1x2⋯xt}. Each time a new data point is acquired, the Viterbi algorithm is used to estimate the current state st*=argmax1≤i≤Nδt(i) as specified in (34).An example of execution of the condition monitoring experiment is shown in Figure4, for both, continuous and discrete observations, respectively. In Figure 4(a) the state duration has been modeled with a Gamma distribution, while in Figure 4(b) with a Gaussian distribution. In Figures 4(a) and 4(b), the first display represents the true state of the HSMM and the second display represents the estimated state from the Viterbi algorithm, while the third display represents the observed time series.Condition monitoring using the Viterbi path. HSMMs can be effective to solve condition monitoring problems in time-dependent applications due to their high accuracy in hidden state recognition.
(a)
State estimation with Viterbi path for continuous data and Gamma duration distribution
(b)
State estimation with Viterbi path for discrete data and Gaussian duration distributionKnowing the true state sequence we calculated the accuracy, defined as the percentage of correctly estimated states over the total length of the state sequence, for each of the testing cases. The results are summarized in Table1(a) for the continuous observations and in Table 1(b) for the discrete observations. The high percentage of correct classified states shows that HSMMs can be effectively used to solve condition monitoring problems for applications in which the system shows a strong time and state duration dependency.Table 1
State recognition accuracy.
(a)
Continuous observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
99.4%
98.5%
99.2%
2
99.7%
98.6%
99.5%
3
99.4%
99.2%
99.7%
4
98.9%
98.9%
99.7%
5
98.2%
98.9%
100%
6
99.1%
98.8%
99.7%
7
98.5%
99.4%
99.7%
8
99.2%
99.1%
99.5%
9
99.2%
98.6%
99.7%
10
99.2%
99.1%
99.5%
Average
99.1%
98.9%
99.6%
(b)
Discrete observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
97.4%
96.7%
97.4%
2
97.2%
97.6%
96.5%
3
99.4%
95.8%
96.6%
4
98.2%
95.3%
97.7%
5
99.1%
97.4%
97.5%
6
97.8%
97.7%
97.8%
7
95.8%
97.2%
96.6%
8
97.7%
96.4%
97.2%
9
98.9%
97.2%
98.5%
10
99.2%
95.6%
96.9%
Average
98.1%
96.7%
97.3%


### 5.1.4. Remaining Useful Lifetime Estimation

In this experimental phase we considered the stateS5 as the failure state and the trained parameters λ* of Section 5.1.2 for each HSMM configuration. As the online RUL estimation procedure is intended to be used in real time, we simulated condition monitoring experiment where we progressively consider the observations {x1x2⋯xt} up to time t. When a new observation is acquired, after the current state probability δ¯t(i) is estimated (Equation (48)), the calculation of the average, upper, and lower RUL ((57), (58), and (59)) is performed.Examples of RUL estimation are illustrated in Figure5. In particular Figure 5(a) represents the case of continuous observations and duration modeled by a Weibull distribution, while Figure 5(b) shows the case of discrete observations and duration modeled by a Gamma distribution. From the figures one can notice that the average, as well as the lower and the upper bound estimations, converges to the real RUL as the failure time approaches. Moreover, as expected, the uncertainty about the estimation decreases with the time, since predictions performed at an early stage are more imprecise. As a consequence, the upper and the lower bound become more narrow as the failure state approaches, and the estimation becomes more precise until it converges to the actual RUL value, with the prediction error tending to zero at the end of the evaluation.HSMMs effectively solve RUL estimation problems. The prediction converges to the actual RUL value and its uncertainty decreases as the real failure time approaches.
(a)
Remaining Useful Lifetime estimation for continuous data and Weibull duration distribution
(b)
Remaining Useful Lifetime estimation for discrete data and Gamma duration distributionTo quantitatively estimate the performance of our methodology for the RUL estimation, we considered at each timet the absolute prediction error (APE) between the real RUL and the predicted value defined as
(63)APEt=RULrealt-RULt,
where RULreal(t) is the (known) value of the RUL at time t, while RUL(t) is RUL predicted by the model. To evaluate the overall performance of our methodology, we considered the average absolute prediction error of the RUL estimation, defined as
(64)APE¯=∑t=1T‍APE(t)T,
where T is the length of the testing signal. APE¯ being a prediction error average, values of (64) close to zero correspond to good predictive performances.The results for each of the 10 testing cases and the different HSMM configurations are reported in Tables2(a) and 2(b), for the continuous and the discrete observation cases, respectively. As it can be noticed, the online Remaining Useful Lifetime estimation and in general the online prediction of the time to a certain event can be effectively faced with HSMMs, which achieve a reliable estimation power with a small prediction error.Table 2
Average absolute prediction error (APE¯) of the RUL estimation using the proposed state duration estimator of (20).
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
5.1
17.0
6.7
14.0
29.0
0.91
4.5
17.0
8.1
2
7.6
19.0
5.0
6.1
21.0
8.5
6.6
19.0
6.1
3
7.7
5.4
19.0
2.9
12.0
17.0
16.0
29.0
3.0
4
9.0
21.0
2.9
7.5
22.0
6.8
6.0
19.0
6.7
5
7.3
19.0
4.7
2.2
14.0
14.0
3.9
17.0
8.7
6
6.5
18.0
5.6
5.1
18.0
10.0
14.0
27.0
2.7
7
4.7
16.0
7.5
4.8
17.0
11.0
1.2
13.0
12.0
8
10.0
22.0
2.9
5.2
18.0
10.0
9.2
22.0
3.9
9
3.1
9.2
14.0
2.0
16.0
13.0
8.2
21.0
4.9
10
6.4
18.0
5.6
7.5
22.0
6.9
3.3
12.0
13.0
Average
6.8
17.0
7.4
5.7
19.0
9.9
7.3
20.0
7.0
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
2.1
11.0
14.0
3.1
8.8
14.0
2.4
12.0
13.0
2
2.1
11.0
13.0
11.0
22.0
3.3
19.0
32.0
7.1
3
5.1
17.0
7.6
6.6
18.0
5.1
2.3
14.0
11.0
4
5.9
6.5
18.0
5.2
17.0
6.7
4.2
16.0
9.0
5
3.2
14.0
10.0
8.3
19.0
3.4
12.0
24.0
2.9
6
12.0
24.0
2.7
6.2
18.0
5.2
4.1
8.4
16.0
7
2.9
15.0
9.7
9.3
21.0
2.3
19.0
31.0
6.6
8
15.0
27.0
7.0
7.4
18.0
4.3
4.3
17.0
9.4
9
5.9
18.0
7.7
11.0
23.0
5.5
3.9
16.0
8.8
10
3.5
11.0
14.0
5.5
6.0
16.0
5.2
17.0
7.1
Average
5.7
15.0
10.0
7.4
17.0
6.6
7.7
19.0
9.0Finally, we tested our RUL estimation methodology using the state duration estimation of (16) introduced by Azimi et al. [30–32]. The results are shown in Tables 3(a) and 3(b), in which, respectively, the prediction errors obtained for continuous and discrete observations are reported.Table 3
Average absolute prediction error (APE¯) of the RUL estimation using the state duration estimator of (16) introduced by Azimi et al. [30–32].
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
57.8
51.0
66.8
26.2
9.7
52.7
25.9
28.4
64.6
2
50.2
44.4
57.7
21.3
17.0
46.9
29.0
19.2
70.8
3
50.3
44.7
57.3
27.1
8.7
56.5
34.5
13.9
73.4
4
51.8
46.0
60.4
21.3
14.3
45.9
34.9
17.1
78.7
5
59.4
53.7
66.2
29.0
9.5
55.4
33.4
15.6
74.9
6
58.0
51.7
67.1
25.8
8.3
54.1
23.1
25.8
66.5
7
59.4
53.6
66.9
18.2
12.5
47.7
36.0
17.1
74.4
8
63.4
55.6
72.3
19.4
15.7
44.1
34.8
17.8
77.0
9
49.1
43.5
57.0
14.5
17.1
43.2
25.1
26.7
67.0
10
54.4
48.4
62.8
23.2
7.9
52.7
24.1
24.5
67.4
Average
55.4
49.3
63.5
22.6
12.1
49.9
30.1
20.6
71.5
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
51.4
41.0
62.4
42.4
31.8
53.0
32.6
26.4
73.6
2
49.6
39.9
60.4
59.5
48.3
70.8
31.3
27.6
69.3
3
50.2
38.6
62.3
46.5
35.7
57.4
32.4
25.7
70.2
4
42.2
31.5
53.8
50.1
40.5
60.6
23.7
36.1
60.3
5
44.3
33.9
55.8
47.8
37.4
59.1
36.0
25.6
76.5
6
52.2
43.2
62.7
55.2
44.3
66.9
27.2
31.6
64.3
7
55.0
43.9
66.8
56.0
45.7
67.0
34.7
23.2
74.4
8
50.3
39.0
62.0
60.4
50.5
71.0
35.1
26.4
72.4
9
55.5
47.4
64.0
48.0
37.2
59.5
31.8
22.2
73.6
10
49.0
38.2
60.7
52.1
41.2
63.1
29.4
28.9
68.7
Average
50.0
39.7
61.1
51.8
41.3
62.9
31.4
27.4
70.4Comparing Table2 and Table 3, one can notice that the proposed RUL method outperforms the one of Azimi. This is mainly due to the proposed average state duration of (20), compared to the one of Azimi given by (16).


## 5.1.1. Data Generation

The industrial machine subject of these experiments has been modeled as a left-right parametric HSMM, withN=5 states, having state S5 as absorbing (failure) state. The choice of a left-right setting has been made for simplicity reasons, since the primary goal of this work is to demonstrate that the proposed model specification coupled with the Akaike Information Criterion is effective to solve automatic model selection, online condition monitoring, and prediction problems. At this purpose, we divided the experiments in two cases according to the nature of the observation, being continuous or discrete.For each of the continuous and the discrete cases, three data sets have been generated by considering the following duration models: Gaussian, Gamma, and Weibull densities. For each of the three data sets, 30 series of data are used as training set and 10 series as testing set. Each time series containsT=650 observations. The parameters used to generate the data are taken from the work of Lee et al. [15] and are adapted to obtain an equivalent left-right parametric HSMM as follows:
(60)π=10000,A0=0100000100000100000100001,ΘN=θ1=100;20,θ2=90;15,θ3=100;20,θ4=80;25,θ5=200;1,ΘG=θ1=500;0.2,θ2=540;0.1667,θ3=500;0.2,θ4=256;0.3125,θ5=800;0.005,ΘW=θ1=102;28,θ2=92;29,θ3=102;28,θ4=82;20,θ5=200;256,
where ΘN, ΘG, and ΘW are the different distribution parameters used for the Gaussian, Gamma, and Weibull duration models, respectively. More precisely, they represent the values of the mean μd and the the variance σd2 of the Gaussian distribution, the shape νd and the scale ηd of the Gamma distribution, and the scale ad and the shape bd of the Weibull distribution. It must be noticed that, as explained in Section 2.1, for state S5, being the absorbing state, the duration parameters θ5 have no influence on the data, since, once the state S5 is reached, the system will remain there forever.Concerning the continuous observation modeling, a bivariate Gaussian distribution has been used with the following parameters [15]:
(61)μ1=2020,μ2=2035,μ3=3535,μ5=2828,U1=200020,U2=150015,U3=15-2-215,U4=5005,U5=103310
while for the discrete case, L=7 distinct observation symbols have been taken into consideration with the following observation probability distribution
(62)B=0.80.2000000.10.80.1000000.10.80.1000000.10.70.10.100000.20.60.10.1.
An example of simulated data both for the continuous and the discrete cases is shown in Figure 2, where a Gaussian duration model has been used.The data generated with the parameter described in Section5.1.1, both for the continuous and the discrete case.
(a)
Example of simulated data for the continuous case
(b)
Example of simulated data for the discrete case


## 5.1.2. Training and Model Selection

The goal of this experimental phase is to test the effectiveness of the AIC in solving the automatic model selection. For this purpose, the training sets of the 6 data sets (continuous/discrte observation with Gaussian, Gamma, and Weibull duration models) have been taken individually and, for each one of them, a series of learning procedure has been run, each one with a variable HSMM structure. In particular we took into account all the combinations of the duration distribution families (Gaussian, Gamma, and Weibull), an increasing number of states,N, from 2 to 8 and, for the continuous observation cases, an increasing number of Gaussian mixtures, M, in the observation distribution from 1 to 4.As accurate parameter initialization is crucial for obtaining a good model fitting [14], a series of 40 learning procedures corresponding to 40 random initializations of the initial parameters λ0 have been executed for each of the considered HSMM structures. For each model structure, the AIC value, as defined in (46), has been evaluated. The final trained set of parameters λ* corresponding to the minimum AIC value has been retained, resulting in 7 HSMMs with a number of states from 2 to 8.The obtained results are shown in Figure3, for both, the continuous and discrete observation data. As it can be noticed, for all the 6 test cases of Figure 3 the AIC values do not improve much for a number of states higher than 5, meaning that adding more states does not add considerable information to the HSMM modeling power. Hence 5 states can be considered as an optimal number of states. Moreover, as shown in the zoomed sections of Figure 3, for the HSMMs with 5 states, the minimum AIC values are obtained for the duration distributions corresponding to the ones used to generate the data. As a consequence AIC can be considered as an effective approach to perform model selection for HSMM, as well as selecting the appropriate parametric distribution family for the state duration modeling.Akaike Information Criterion (AIC) applied to continuous and discrete observations data. AIC is effective for automatic model selection, since its minimum value provides the same number of states and duration model used to generate the data.
(a)
AIC values for continuous data and Gaussian duration distribution
(b)
AIC values for continuous data and Gamma duration distribution
(c)
AIC values for continuous data and Weibull duration distribution
(d)
AIC values for discrete data and Gaussian duration distribution
(e)
AIC values for discrete data and Gamma duration distribution
(f)
AIC values for discrete data and Weibull duration distribution


## 5.1.3. Condition Monitoring

The optimal parametersλ* obtained in the previous phase have been used to perform the online condition monitoring experiment on the 10 testing cases for all the 6 considered HSMM configurations. In this experiment, we simulate online monitoring by considering all the testing observations up to the current time, that is, {x1x2⋯xt}. Each time a new data point is acquired, the Viterbi algorithm is used to estimate the current state st*=argmax1≤i≤Nδt(i) as specified in (34).An example of execution of the condition monitoring experiment is shown in Figure4, for both, continuous and discrete observations, respectively. In Figure 4(a) the state duration has been modeled with a Gamma distribution, while in Figure 4(b) with a Gaussian distribution. In Figures 4(a) and 4(b), the first display represents the true state of the HSMM and the second display represents the estimated state from the Viterbi algorithm, while the third display represents the observed time series.Condition monitoring using the Viterbi path. HSMMs can be effective to solve condition monitoring problems in time-dependent applications due to their high accuracy in hidden state recognition.
(a)
State estimation with Viterbi path for continuous data and Gamma duration distribution
(b)
State estimation with Viterbi path for discrete data and Gaussian duration distributionKnowing the true state sequence we calculated the accuracy, defined as the percentage of correctly estimated states over the total length of the state sequence, for each of the testing cases. The results are summarized in Table1(a) for the continuous observations and in Table 1(b) for the discrete observations. The high percentage of correct classified states shows that HSMMs can be effectively used to solve condition monitoring problems for applications in which the system shows a strong time and state duration dependency.Table 1
State recognition accuracy.
(a)
Continuous observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
99.4%
98.5%
99.2%
2
99.7%
98.6%
99.5%
3
99.4%
99.2%
99.7%
4
98.9%
98.9%
99.7%
5
98.2%
98.9%
100%
6
99.1%
98.8%
99.7%
7
98.5%
99.4%
99.7%
8
99.2%
99.1%
99.5%
9
99.2%
98.6%
99.7%
10
99.2%
99.1%
99.5%
Average
99.1%
98.9%
99.6%
(b)
Discrete observations
Test case
Duration distribution
Gaussian
Gamma
Weibull
1
97.4%
96.7%
97.4%
2
97.2%
97.6%
96.5%
3
99.4%
95.8%
96.6%
4
98.2%
95.3%
97.7%
5
99.1%
97.4%
97.5%
6
97.8%
97.7%
97.8%
7
95.8%
97.2%
96.6%
8
97.7%
96.4%
97.2%
9
98.9%
97.2%
98.5%
10
99.2%
95.6%
96.9%
Average
98.1%
96.7%
97.3%


## 5.1.4. Remaining Useful Lifetime Estimation

In this experimental phase we considered the stateS5 as the failure state and the trained parameters λ* of Section 5.1.2 for each HSMM configuration. As the online RUL estimation procedure is intended to be used in real time, we simulated condition monitoring experiment where we progressively consider the observations {x1x2⋯xt} up to time t. When a new observation is acquired, after the current state probability δ¯t(i) is estimated (Equation (48)), the calculation of the average, upper, and lower RUL ((57), (58), and (59)) is performed.Examples of RUL estimation are illustrated in Figure5. In particular Figure 5(a) represents the case of continuous observations and duration modeled by a Weibull distribution, while Figure 5(b) shows the case of discrete observations and duration modeled by a Gamma distribution. From the figures one can notice that the average, as well as the lower and the upper bound estimations, converges to the real RUL as the failure time approaches. Moreover, as expected, the uncertainty about the estimation decreases with the time, since predictions performed at an early stage are more imprecise. As a consequence, the upper and the lower bound become more narrow as the failure state approaches, and the estimation becomes more precise until it converges to the actual RUL value, with the prediction error tending to zero at the end of the evaluation.HSMMs effectively solve RUL estimation problems. The prediction converges to the actual RUL value and its uncertainty decreases as the real failure time approaches.
(a)
Remaining Useful Lifetime estimation for continuous data and Weibull duration distribution
(b)
Remaining Useful Lifetime estimation for discrete data and Gamma duration distributionTo quantitatively estimate the performance of our methodology for the RUL estimation, we considered at each timet the absolute prediction error (APE) between the real RUL and the predicted value defined as
(63)APEt=RULrealt-RULt,
where RULreal(t) is the (known) value of the RUL at time t, while RUL(t) is RUL predicted by the model. To evaluate the overall performance of our methodology, we considered the average absolute prediction error of the RUL estimation, defined as
(64)APE¯=∑t=1T‍APE(t)T,
where T is the length of the testing signal. APE¯ being a prediction error average, values of (64) close to zero correspond to good predictive performances.The results for each of the 10 testing cases and the different HSMM configurations are reported in Tables2(a) and 2(b), for the continuous and the discrete observation cases, respectively. As it can be noticed, the online Remaining Useful Lifetime estimation and in general the online prediction of the time to a certain event can be effectively faced with HSMMs, which achieve a reliable estimation power with a small prediction error.Table 2
Average absolute prediction error (APE¯) of the RUL estimation using the proposed state duration estimator of (20).
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
5.1
17.0
6.7
14.0
29.0
0.91
4.5
17.0
8.1
2
7.6
19.0
5.0
6.1
21.0
8.5
6.6
19.0
6.1
3
7.7
5.4
19.0
2.9
12.0
17.0
16.0
29.0
3.0
4
9.0
21.0
2.9
7.5
22.0
6.8
6.0
19.0
6.7
5
7.3
19.0
4.7
2.2
14.0
14.0
3.9
17.0
8.7
6
6.5
18.0
5.6
5.1
18.0
10.0
14.0
27.0
2.7
7
4.7
16.0
7.5
4.8
17.0
11.0
1.2
13.0
12.0
8
10.0
22.0
2.9
5.2
18.0
10.0
9.2
22.0
3.9
9
3.1
9.2
14.0
2.0
16.0
13.0
8.2
21.0
4.9
10
6.4
18.0
5.6
7.5
22.0
6.9
3.3
12.0
13.0
Average
6.8
17.0
7.4
5.7
19.0
9.9
7.3
20.0
7.0
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
2.1
11.0
14.0
3.1
8.8
14.0
2.4
12.0
13.0
2
2.1
11.0
13.0
11.0
22.0
3.3
19.0
32.0
7.1
3
5.1
17.0
7.6
6.6
18.0
5.1
2.3
14.0
11.0
4
5.9
6.5
18.0
5.2
17.0
6.7
4.2
16.0
9.0
5
3.2
14.0
10.0
8.3
19.0
3.4
12.0
24.0
2.9
6
12.0
24.0
2.7
6.2
18.0
5.2
4.1
8.4
16.0
7
2.9
15.0
9.7
9.3
21.0
2.3
19.0
31.0
6.6
8
15.0
27.0
7.0
7.4
18.0
4.3
4.3
17.0
9.4
9
5.9
18.0
7.7
11.0
23.0
5.5
3.9
16.0
8.8
10
3.5
11.0
14.0
5.5
6.0
16.0
5.2
17.0
7.1
Average
5.7
15.0
10.0
7.4
17.0
6.6
7.7
19.0
9.0Finally, we tested our RUL estimation methodology using the state duration estimation of (16) introduced by Azimi et al. [30–32]. The results are shown in Tables 3(a) and 3(b), in which, respectively, the prediction errors obtained for continuous and discrete observations are reported.Table 3
Average absolute prediction error (APE¯) of the RUL estimation using the state duration estimator of (16) introduced by Azimi et al. [30–32].
(a)
APE
¯ of the RUL estimation for the continuous observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
57.8
51.0
66.8
26.2
9.7
52.7
25.9
28.4
64.6
2
50.2
44.4
57.7
21.3
17.0
46.9
29.0
19.2
70.8
3
50.3
44.7
57.3
27.1
8.7
56.5
34.5
13.9
73.4
4
51.8
46.0
60.4
21.3
14.3
45.9
34.9
17.1
78.7
5
59.4
53.7
66.2
29.0
9.5
55.4
33.4
15.6
74.9
6
58.0
51.7
67.1
25.8
8.3
54.1
23.1
25.8
66.5
7
59.4
53.6
66.9
18.2
12.5
47.7
36.0
17.1
74.4
8
63.4
55.6
72.3
19.4
15.7
44.1
34.8
17.8
77.0
9
49.1
43.5
57.0
14.5
17.1
43.2
25.1
26.7
67.0
10
54.4
48.4
62.8
23.2
7.9
52.7
24.1
24.5
67.4
Average
55.4
49.3
63.5
22.6
12.1
49.9
30.1
20.6
71.5
(b)
APE
¯ of the RUL estimation for the discrete observation test cases
Test case
Duration distribution
Gaussian
Gamma
Weibull
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
APE
¯
avg
up
low
avg
up
low
avg
up
low
1
51.4
41.0
62.4
42.4
31.8
53.0
32.6
26.4
73.6
2
49.6
39.9
60.4
59.5
48.3
70.8
31.3
27.6
69.3
3
50.2
38.6
62.3
46.5
35.7
57.4
32.4
25.7
70.2
4
42.2
31.5
53.8
50.1
40.5
60.6
23.7
36.1
60.3
5
44.3
33.9
55.8
47.8
37.4
59.1
36.0
25.6
76.5
6
52.2
43.2
62.7
55.2
44.3
66.9
27.2
31.6
64.3
7
55.0
43.9
66.8
56.0
45.7
67.0
34.7
23.2
74.4
8
50.3
39.0
62.0
60.4
50.5
71.0
35.1
26.4
72.4
9
55.5
47.4
64.0
48.0
37.2
59.5
31.8
22.2
73.6
10
49.0
38.2
60.7
52.1
41.2
63.1
29.4
28.9
68.7
Average
50.0
39.7
61.1
51.8
41.3
62.9
31.4
27.4
70.4Comparing Table2 and Table 3, one can notice that the proposed RUL method outperforms the one of Azimi. This is mainly due to the proposed average state duration of (20), compared to the one of Azimi given by (16).


## 5.2. Real Data

In this section we apply the proposed HSMM-based approach for RUL estimation to a real case study using bearing monitoring data recorded during experiments carried out on the Pronostia experimental platform and made available for the IEEE Prognostics and Health Management (PHM) 2012 challenge [49]. The data correspond to normally degraded bearings, leading to cases which closely correspond to the industrial reality.The choice of testing the proposed methodology on bearings derives from two facts: (i) bearings are the most critical components related to failures of rotating machines [50] and (ii) their monotonically increasing degradation pattern justifies the usage of left-right HSMM models.


### 5.2.1. Data Description

Pronostia is an experimental platform designed and realized at the Automatic Control and Micro-Mechatronic Systems (AS2M) Department of Franche-Comté Electronics, Mechanics, Thermal Processing, Optics-Science and Technology (FEMTO-ST) Institute (http://www.femto-st.fr/) (Besançon, France), with the aim of collecting real data related to accelerated degradation of bearings. Such data are used to validate methods for bearing condition assessment, diagnostic and prognostic [19, 51–59].The Pronostia platform allows to perform run-to-failure tests under constant or variable operating conditions. The operating conditions are determined by two factors that can be controlled online: (i) the rotation speed and (ii) the radial force load. During each experiment, temperature and vibration monitoring measurements are gathered online, through two type of sensors placed in the bearing housing: a temperature probe and two accelerometers (one on the vertical and one on the horizontal axis).The platform is composed of three main parts: a rotating part, a load profile generation part, and a measurement part, as illustrated in Figure6.Figure 6
Global overview of the Pronostia experimental platform [19].The rotating part is composed of an asynchronous motor which develops a power equal to 250 W, two shafts, and a gearbox, which allows the motor to reach its rated speed of 2830 rpm. The motor’s rotation speed and the direction are set through a human machine interface.The load profiles part issues a radial force on the external ring of the test bearing through a pneumatic jack connected to a lever arm, which indirectly transmits the load through a clamping ring. The goal of the applied radial force is to accelerate the bearing’s degradation process.The measurement part consists of a data acquisition card connected to the monitoring sensors, which provides the user with the measured temperature and vibration data. The vibration measurements are provided in snapshots of 0.1 s collected each 10 seconds at a sampling frequency of 25.6 kHz (2560 samples per each snapshot), while the temperature has been continuously recorded at a sampling frequency of 10 Hz (600 samples collected each minute).Further details on the Pronostia test rig can be found on the data presentation paper [49] and on the web page of the data challenge (http://www.femto-st.fr/en/Research-departments/AS2M/Research-groups/PHM/IEEE-PHM-2012-Data-challenge.php).Regarding the data provided for the PHM 2012 challenge, 3 different operating conditions were considered:(i)
first operating conditions: speed of 1800 rpm and load of 4000 Newton;(ii)
second operating conditions: speed of 1650 rpm and load of 4200 Newton;(iii)
third operating conditions: speed of 1500 rpm and load of 5000 Newton.Under the above operating conditions, a total of 17 accelerated life tests were realized on bearings of type NSK 6804 DD, which can operate at a maximum speed of 13000 rpm and a load limit of 4000 N. The tests were stopped when the amplitude of the vibration signal was higher than 20 g; thus this moment was defined as the bearing failure time. An example of bearing before and after the experiment is shown in Figure7 together with the corresponding vibration signal collected during the whole test.Figure 7
A tested bearing before and after the experiment with its recorded vibration signal [49].Table4 reports how the 17 tested bearings were separated into the three operating conditions. Moreover, the duration of each experiment, being the RUL to be predicted for each bearing, is also given. We performed two sets of experiments by considering, respectively, the bearings relative to the first and the second operating condition (i.e.,Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Table 4
Lifetime duration (in seconds) and operating conditions of the bearings tested in the IEEE PHM 2012 Prognostic Challenge dataset [49].
Condition 1
Condition 2
Condition 3
1800 rpm and 4000 N
1650 rpm and 4200 N
1500 rpm and 5000 N
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing1_1
28030
Bearing2_1
9110
Bearing3_1
5150
Bearing1_2
8710
Bearing2_2
7970
Bearing3_2
16370
Bearing1_3
23750
Bearing2_3
19550
Bearing3_3
4340
Bearing1_4
14280
Bearing2_4
7510
Bearing1_5
24630
Bearing2_5
23110
Bearing1_6
24480
Bearing2_6
7010
Bearing1_7
22590
Bearing2_7
2300As already mentioned, the available data correspond to normally degraded bearings, meaning that the defects were not initially induced and that each degraded bearing contains almost all the types of defects (balls, rings, and cage), resembling faithfully a common real industrial situation. Moreover, no assumption about the type of failure to be occurred is provided with the data and, since the variability in experiment durations is high (from 1 h to 7 h), performing good estimates of the RUL is a difficult task [49].In our experiments we considered, as input to our model, the horizontal channel of the accelerometer. We preprocessed the raw signals by extracting two time-domain features, that is, root mean square (RMS) and kurtosis, within windows of the same length as the given snapshots (L=2560). Let rw(t) be the raw signal of the wth window; for each w we estimate RMS as xwRMS=1/L∑t=1L‍rw2(t) and kurtosis as xwKURT=1/L∑t=1L‍rw(t)-r¯w4/1/L∑t=1L‍rw(t)-r¯w22, where r¯w is the mean of rw. An example of feature extraction forBearing1_1 is shown in Figure 8.Figure 8
Raw vibration data (a) versus RMS and kurtosis features (b) forBearing1_1.To assess the performance of the proposed HSMM, after the model selection procedure, we implemented a leave-one-out cross validation scheme: by considering separately conditions 1 and 2, we performed the online RUL estimation for each of the 7 bearings, using an HSMM trained with the remaining 6 bearing histories. Similarly to the simulated case, we considered the average absolute prediction error, defined in (64), to quantitatively evaluate our method.


### 5.2.2. Bearings RUL Estimation

We performed our experiments in two steps: firstly we applied model selection in order to determine an optimal model structure, and secondly, we estimated the RUL of the bearings. The full procedure is detailed in the following.(A) HSMM Structure. To determine an appropriate HSMM structure for effectively modeling the considered data, we considered several HSMM structures, characterized by (i) the duration distribution family (being Gaussian, Gamma, or Weibull), (ii) an increasing number of states, N, from 2 to 6, and (iii) an increasing number of Gaussian mixtures, M, in the observation density from 1 to 4. For each model structure, obtained by systematically considering all the combinations of (i) to (iii), we run 120 parameter learnings, corresponding to 120 random initializations, λ0, on the data sets (Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Similar to Section5.1.2, at the end of each learning we evaluated the AIC value (Equation (46)) as reported in Figures 9(a) and 9(b) for conditions 1 and 2, respectively. In both cases, the global minimum AIC value corresponds to an HSMM with N=4 states, a Weibull duration model, and a M=1 Gaussians mixture for the observation density.In both cases the minimum AIC value corresponds to an HSMM withN=4 states, a Weibull duration model, and M=1 mixture in the observation density.
(a)
AIC values for Condition 1
(b)
AIC values for Condition 2(B) RUL Estimation. Using the above obtained optimal HSMM structure, we trained it via a leave-one-out cross validation scheme by using for condition 1, at each iteration,Bearing1_i, 1≤i≤7, as the testing bearing, while the remaining six bearings were used for training. Once the trained parameters λi* were estimated for the ith testing bearing, we progressively collected the observations of the testedBearing1_i to calculate, at each time t, the average, lower, and upper RUL, as specified in (57), (58), and (59), respectively, considering the state S4 as the failure state. The same procedure has been performed for the bearings in condition 2.Examples of RUL estimation forBearing1_7 andBearing2_6 are shown in Figures 10(a) and 10(b), respectively, where the black solid line represents the real RUL which goes to zero as the time goes on. As it can be seen, the average as well as the lower and the upper bound estimations converge to the real RUL as the real failure time approaches and the uncertainty about the estimation decreases with time.By obtaining a low average absolute prediction error, the proposed parametric HSMM is effective for estimating the Remaining Useful Lifetime of bearings.
(a)
RUL estimation for Bearing1_7
(b)
RUL estimation for Bearing2_6Concerning the quantitative estimation of the predictive performances, we report in Table5 the average absolute prediction error of the RUL estimation (see Equation (64)), expressed in seconds. As it can be noticed, the average absolute prediction error of the average RUL is, respectively, 1 hour and 15 minutes for condition 1, and 1 hour and 5 minutes for condition 2, which are good values, considering the high variability in the training set durations and the fact that the performance metric takes into account also the less accurate predictions performed at an early stage of the bearings life. Moreover, for both conditions, the average prediction errors of 5 tests out of 7 are below the average, while the best average error of the mean RUL is only 23 minutes for condition 1 while it further decreases to 14 minutes for condition 2.Table 5
Average absolute prediction error (APE¯) of the RUL estimation, expressed in seconds.
(a)
Condition 1
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing1_1
10571.6
12723.0
9414.6
Bearing1_2
4331.2
3815.6
3821.3
Bearing1_3
2997.0
9730.9
6091.2
Bearing1_4
6336.3
2876.6
14871.9
Bearing1_5
1968.9
7448.4
10411.5
Bearing1_6
4253.0
9896.4
9793.7
Bearing1_7
1388.0
7494.3
10088.1
Average
4549.4
7712.2
9213.2
(b)
Condition 2
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing2_1
2475.9
5006.5
7287.5
Bearing2_2
1647.3
4497.2
8288.6
Bearing2_3
8877.1
9508.3
7962.1
Bearing2_4
1769.8
4248.6
4982.5
Bearing2_5
8663.1
10490.0
10730.0
Bearing2_6
877.1
3504.7
6687.0
Bearing2_7
3012.5
3866.4
6651.9
Average
3903.3
5874.5
7512.8


## 5.2.1. Data Description

Pronostia is an experimental platform designed and realized at the Automatic Control and Micro-Mechatronic Systems (AS2M) Department of Franche-Comté Electronics, Mechanics, Thermal Processing, Optics-Science and Technology (FEMTO-ST) Institute (http://www.femto-st.fr/) (Besançon, France), with the aim of collecting real data related to accelerated degradation of bearings. Such data are used to validate methods for bearing condition assessment, diagnostic and prognostic [19, 51–59].The Pronostia platform allows to perform run-to-failure tests under constant or variable operating conditions. The operating conditions are determined by two factors that can be controlled online: (i) the rotation speed and (ii) the radial force load. During each experiment, temperature and vibration monitoring measurements are gathered online, through two type of sensors placed in the bearing housing: a temperature probe and two accelerometers (one on the vertical and one on the horizontal axis).The platform is composed of three main parts: a rotating part, a load profile generation part, and a measurement part, as illustrated in Figure6.Figure 6
Global overview of the Pronostia experimental platform [19].The rotating part is composed of an asynchronous motor which develops a power equal to 250 W, two shafts, and a gearbox, which allows the motor to reach its rated speed of 2830 rpm. The motor’s rotation speed and the direction are set through a human machine interface.The load profiles part issues a radial force on the external ring of the test bearing through a pneumatic jack connected to a lever arm, which indirectly transmits the load through a clamping ring. The goal of the applied radial force is to accelerate the bearing’s degradation process.The measurement part consists of a data acquisition card connected to the monitoring sensors, which provides the user with the measured temperature and vibration data. The vibration measurements are provided in snapshots of 0.1 s collected each 10 seconds at a sampling frequency of 25.6 kHz (2560 samples per each snapshot), while the temperature has been continuously recorded at a sampling frequency of 10 Hz (600 samples collected each minute).Further details on the Pronostia test rig can be found on the data presentation paper [49] and on the web page of the data challenge (http://www.femto-st.fr/en/Research-departments/AS2M/Research-groups/PHM/IEEE-PHM-2012-Data-challenge.php).Regarding the data provided for the PHM 2012 challenge, 3 different operating conditions were considered:(i)
first operating conditions: speed of 1800 rpm and load of 4000 Newton;(ii)
second operating conditions: speed of 1650 rpm and load of 4200 Newton;(iii)
third operating conditions: speed of 1500 rpm and load of 5000 Newton.Under the above operating conditions, a total of 17 accelerated life tests were realized on bearings of type NSK 6804 DD, which can operate at a maximum speed of 13000 rpm and a load limit of 4000 N. The tests were stopped when the amplitude of the vibration signal was higher than 20 g; thus this moment was defined as the bearing failure time. An example of bearing before and after the experiment is shown in Figure7 together with the corresponding vibration signal collected during the whole test.Figure 7
A tested bearing before and after the experiment with its recorded vibration signal [49].Table4 reports how the 17 tested bearings were separated into the three operating conditions. Moreover, the duration of each experiment, being the RUL to be predicted for each bearing, is also given. We performed two sets of experiments by considering, respectively, the bearings relative to the first and the second operating condition (i.e.,Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Table 4
Lifetime duration (in seconds) and operating conditions of the bearings tested in the IEEE PHM 2012 Prognostic Challenge dataset [49].
Condition 1
Condition 2
Condition 3
1800 rpm and 4000 N
1650 rpm and 4200 N
1500 rpm and 5000 N
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing
Lifetime [s]
Bearing1_1
28030
Bearing2_1
9110
Bearing3_1
5150
Bearing1_2
8710
Bearing2_2
7970
Bearing3_2
16370
Bearing1_3
23750
Bearing2_3
19550
Bearing3_3
4340
Bearing1_4
14280
Bearing2_4
7510
Bearing1_5
24630
Bearing2_5
23110
Bearing1_6
24480
Bearing2_6
7010
Bearing1_7
22590
Bearing2_7
2300As already mentioned, the available data correspond to normally degraded bearings, meaning that the defects were not initially induced and that each degraded bearing contains almost all the types of defects (balls, rings, and cage), resembling faithfully a common real industrial situation. Moreover, no assumption about the type of failure to be occurred is provided with the data and, since the variability in experiment durations is high (from 1 h to 7 h), performing good estimates of the RUL is a difficult task [49].In our experiments we considered, as input to our model, the horizontal channel of the accelerometer. We preprocessed the raw signals by extracting two time-domain features, that is, root mean square (RMS) and kurtosis, within windows of the same length as the given snapshots (L=2560). Let rw(t) be the raw signal of the wth window; for each w we estimate RMS as xwRMS=1/L∑t=1L‍rw2(t) and kurtosis as xwKURT=1/L∑t=1L‍rw(t)-r¯w4/1/L∑t=1L‍rw(t)-r¯w22, where r¯w is the mean of rw. An example of feature extraction forBearing1_1 is shown in Figure 8.Figure 8
Raw vibration data (a) versus RMS and kurtosis features (b) forBearing1_1.To assess the performance of the proposed HSMM, after the model selection procedure, we implemented a leave-one-out cross validation scheme: by considering separately conditions 1 and 2, we performed the online RUL estimation for each of the 7 bearings, using an HSMM trained with the remaining 6 bearing histories. Similarly to the simulated case, we considered the average absolute prediction error, defined in (64), to quantitatively evaluate our method.


## 5.2.2. Bearings RUL Estimation

We performed our experiments in two steps: firstly we applied model selection in order to determine an optimal model structure, and secondly, we estimated the RUL of the bearings. The full procedure is detailed in the following.(A) HSMM Structure. To determine an appropriate HSMM structure for effectively modeling the considered data, we considered several HSMM structures, characterized by (i) the duration distribution family (being Gaussian, Gamma, or Weibull), (ii) an increasing number of states, N, from 2 to 6, and (iii) an increasing number of Gaussian mixtures, M, in the observation density from 1 to 4. For each model structure, obtained by systematically considering all the combinations of (i) to (iii), we run 120 parameter learnings, corresponding to 120 random initializations, λ0, on the data sets (Bearing1_1, Bearing1_2, …, Bearing1_7 andBearing2_1, Bearing2_2, …, Bearing2_7).Similar to Section5.1.2, at the end of each learning we evaluated the AIC value (Equation (46)) as reported in Figures 9(a) and 9(b) for conditions 1 and 2, respectively. In both cases, the global minimum AIC value corresponds to an HSMM with N=4 states, a Weibull duration model, and a M=1 Gaussians mixture for the observation density.In both cases the minimum AIC value corresponds to an HSMM withN=4 states, a Weibull duration model, and M=1 mixture in the observation density.
(a)
AIC values for Condition 1
(b)
AIC values for Condition 2(B) RUL Estimation. Using the above obtained optimal HSMM structure, we trained it via a leave-one-out cross validation scheme by using for condition 1, at each iteration,Bearing1_i, 1≤i≤7, as the testing bearing, while the remaining six bearings were used for training. Once the trained parameters λi* were estimated for the ith testing bearing, we progressively collected the observations of the testedBearing1_i to calculate, at each time t, the average, lower, and upper RUL, as specified in (57), (58), and (59), respectively, considering the state S4 as the failure state. The same procedure has been performed for the bearings in condition 2.Examples of RUL estimation forBearing1_7 andBearing2_6 are shown in Figures 10(a) and 10(b), respectively, where the black solid line represents the real RUL which goes to zero as the time goes on. As it can be seen, the average as well as the lower and the upper bound estimations converge to the real RUL as the real failure time approaches and the uncertainty about the estimation decreases with time.By obtaining a low average absolute prediction error, the proposed parametric HSMM is effective for estimating the Remaining Useful Lifetime of bearings.
(a)
RUL estimation for Bearing1_7
(b)
RUL estimation for Bearing2_6Concerning the quantitative estimation of the predictive performances, we report in Table5 the average absolute prediction error of the RUL estimation (see Equation (64)), expressed in seconds. As it can be noticed, the average absolute prediction error of the average RUL is, respectively, 1 hour and 15 minutes for condition 1, and 1 hour and 5 minutes for condition 2, which are good values, considering the high variability in the training set durations and the fact that the performance metric takes into account also the less accurate predictions performed at an early stage of the bearings life. Moreover, for both conditions, the average prediction errors of 5 tests out of 7 are below the average, while the best average error of the mean RUL is only 23 minutes for condition 1 while it further decreases to 14 minutes for condition 2.Table 5
Average absolute prediction error (APE¯) of the RUL estimation, expressed in seconds.
(a)
Condition 1
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing1_1
10571.6
12723.0
9414.6
Bearing1_2
4331.2
3815.6
3821.3
Bearing1_3
2997.0
9730.9
6091.2
Bearing1_4
6336.3
2876.6
14871.9
Bearing1_5
1968.9
7448.4
10411.5
Bearing1_6
4253.0
9896.4
9793.7
Bearing1_7
1388.0
7494.3
10088.1
Average
4549.4
7712.2
9213.2
(b)
Condition 2
Test Bearings
APE
¯
avg
APE
¯
low
APE
¯
up
Bearing2_1
2475.9
5006.5
7287.5
Bearing2_2
1647.3
4497.2
8288.6
Bearing2_3
8877.1
9508.3
7962.1
Bearing2_4
1769.8
4248.6
4982.5
Bearing2_5
8663.1
10490.0
10730.0
Bearing2_6
877.1
3504.7
6687.0
Bearing2_7
3012.5
3866.4
6651.9
Average
3903.3
5874.5
7512.8


## 6. Conclusion and Future Work

In this paper, we introduced an approach based on Hidden Semi-Markov Models (HSMM) and Akaike Information Criteria (AIC) to perform (i) automatic model selection, (ii) online condition monitoring, and (iii) online time to event estimation.The proposed HSMM models the state duration distribution with a parametric density, allowing a less computationally expensive learning procedure due to the few required parameters to estimate. Together with the provided general model specification, the modified learning, inference, and prediction algorithms allow the usage of any parametric distribution to model the state duration, as well as continuous or discrete observations. As a consequence, a wide class of different applications can be modeled with the proposed methodology.This paper highlights, through experiments performed on simulated data, that the proposed approach is effective in (i) automatically selecting the correct configuration of the HSMM in terms of number of states and correct duration distribution family, (ii) performing online state estimation, and (iii) correctly predict the time to a determinate event, identified as the entrance of the model in a target state. As a consequence, the proposed parametric HSMM combined with AIC can be used in practice for condition monitoring and Remaining Useful Lifetime applications.As the targeted application of the proposed methodology is failure prognosis in industrial machines, combining the proposed HSMM model with online learning procedure, capable of adapting the model parameter to new conditions, would be considered in a future extension.



---

*Source: 278120-2015-02-22.xml*