# Optimal Cost-Effective Maintenance Policy for a Helicopter Gearbox Early Fault Detection under Varying Load

Authors: Xin Li; Jing Cai; Hongfu Zuo; Huaiyuan Li
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Most of the existing fault detection methods rarely consider the cost-optimal maintenance policy. A novel multivariate Bayesian control approach is proposed, which enables the implementation of early fault detection for a helicopter gearbox with cost minimization maintenance policy under varying load. A continuous time hidden semi-Markov model (HSMM) is employed to describe the stochastic relationship between the unobservable states and observable observations of the gear system. Explicit expressions for the remaining useful life prediction are derived using HSMM. Considering the maintenance cost in fault detection, the multivariate Bayesian control scheme based on HSMM is developed; the objective is to minimize the long-run expected average cost per unit time. An effective computational algorithm in the semi-Markov decision process (SMDP) framework is designed to obtain the optimal control limit. A comparison with the multivariate Bayesian control chart based on hidden Markov model (HMM) and the traditional age-based replacement policy is given, which illustrates the effectiveness of the proposed approach.


---

## Body

## 1. Introduction

In a helicopter system, the mechanical drive system is the most efficient and compact device to transmit torque and change angular velocity. As the key mechanical transmission parts, gears are extensively used in aerospace, shipbuilding industry, and wind power industry, such as helicopters, cruisers, and wind turbines. Gear failure will cause machine halts in the whole mechanical system, resulting in great economic losses and even human casualties.Condition monitoring and fault detection technique can significantly improve the reliability of the gear transmission system and reduce the occurrence of failure. In condition-based maintenance (CBM) and prognostics and health management (PHM), the vibration monitoring data obtained from the accelerometers carry a large amount of information, which were used in the fault detection of the gears for decades and were very helpful to prevent unnecessary machine halts [1, 2]. The accurate early fault detection can effectively prevent the occurrence of secondary damage. Appropriate preventive maintenance can significantly improve equipment availability and save cost for gearbox users. Therefore, the early gear fault detection and maintenance scheme optimization have become research hotspots in recent years.In recent research, advanced nonparametric methods were widely applied to gearbox early fault detection, such as the continuous wavelet transform (CWT) [3], the multiscale chirplet path pursuit (MSCPP) [4], and multiresolution Fourier transform [5]. Some researchers applied parametric time series model to gearbox fault detection, most of whom assumed that gearbox ran under constant load conditions. For example, Wang and Wong [6] indicated that the autoregressive (AR) model with a high diagnostic confidence level can detect the gear crack much earlier than the conventional method. Rofe [7] used the autoregressive moving-average (ARMA) model under the condition of load fluctuation for gear fault diagnosis and selected the variance and kurtosis of the residuals of ARMA model as the indicators of failure. It should be noted that the fault diagnosis methods in the literature mostly assume the gears run under constant speed or constant load conditions, meaning that all the signal responses are due to component degradation or failures. This is not always consistent with the reality that gears often run under varying load, such as the wind turbine gearboxes and the helicopter planetary gearboxes [8, 9]. Under the assumption of constant load, it is difficult to separate the interference signals caused by load changes, which greatly affects the characteristics of the signals. In most cases, the change of the load variation in the feature level is similar to that caused by fault. This increases the difficulty of feature extraction and thereby also increases the difficulty of effective monitoring, diagnosis, and prediction. The gear motion residual (GMR) signal is an effective algorithm and highly insensitive to the varying load [10]. The healthy portion of the GMR signals is used to fit the VAR model and the residuals of the whole data sets are calculated to process the early fault detection scheme in this paper.The operational states (healthy or unhealthy) cannot be observed directly when the gear system is running. Several approaches have been used to model the degradation process, such as proportional hazards model (PH) [11], stochastic filtering model (SFM) [12], and hidden Markov model (HMM) [13]. In recent years, HMM has been widely applied to speech recognition [14], system fault diagnosis [13, 15], and life prediction [16]. However, in HMM, the sojourn time in each operational state of the system is assumed to follow exponential distribution, that is, the “memoryless” property. As an extension of the HMM, the sojourn time in each operational state is not restricted by exponential distribution in a hidden semi-Markov model (HSMM), and its degradation process is much closer to the real degradation situations. Liu et al. [17, 18] proposed a diagnosis and prediction method of equipment based on the HSMM and the sojourn time follows normal distribution. Jiang et al. [19] carried out a small extension to HMM, considering a HSMM with the sojourn time obeys Erlang (2,λ) distribution. A general Erlang HSMM is considered in this paper to model the deterioration process of the gear system. The expectation-maximization (EM) algorithm is employed to estimate the unknown state and observation parameters of the model. After the model parameters are estimated, the posterior probability that the system is in each state is used to derive the expressions for the reliability distribution function of remaining useful life (RUL), the probability density function (PDF), and the mean residual life (MRL) function.When using the HSMM to model the deterioration of the gear system with the consideration of maintenance policy, it is necessary to consider the best stopping time. The Bayesian methods are used extensively to determine the stopping time. Girshick and Rubin [20] introduced the Bayesian method for the production process control for the first time in 1952, showing that repair was needed when the posterior probability was beyond the control limit. Considering the minimum maintenance cost, Kim et al. [13] put forward the optimal Bayesian forecasting model using the actual oil spectral analysis data. Most of the stochastic models assumed that the system degradation process was modeled by a continuous time HMM and each sojourn time of hidden states was exponentially distributed, which is not always realistic. Relaxing the assumption of Markovian process deterioration, Panagiotidou and Tagaras [21] proposed a model that integrated process control and maintenance for two operational states. However, this method is only applicable to univariate monitoring. Such a drawback of the existing optimization model motivates us to set up an optimal maintenance model under real circumstances. In order to find the optimal stopping time, based on the HSMM, a Bayesian control scheme is established with the objective of minimizing the long-run expected average cost per unit time to determine the optimal stopping time and when to perform the preventive maintenance actions.To the best of our knowledge, this is the first paper applying HSMM with general Erlang distribution of sojourn time in two hidden states and optimal Bayesian control scheme to find the cost-effective strategy for early fault detection of the helicopter gearbox. With the condition monitoring information, the proposed approach can not only update the remaining useful life at each sampling epoch, but also process early fault detection of the helicopter gear system and determine the optimal stopping time with minimum cost simultaneously.The remainder of the paper is organized as follows. In Section2, the healthy part of GMR signal of the helicopter gearbox life test vibration data is used to fit a VAR model, and the residuals of the whole data sets are obtained. In Section 3, the HSMM is applied, in which each operational state obeys the general Erlang distribution, to describe the degradation process of the gear system, and the EM algorithm is employed to estimate the unknown state and observation parameters of the model. The conditional reliability distribution functions of the RUL and the MRL function are derived. In Section 4, on the basis of hidden semi-Markov degradation modeling, an optimal Bayesian control chart is developed, and the optimal control limit is solved in the semi-Markov decision process framework, followed by a comparison with other methods. The conclusions are presented in Section 5.


## 2. Residuals Computation for a Helicopter Gearbox

### 2.1. Experimental Scheme

The overall process of the optimal Bayesian control scheme for the helicopter gearbox early fault detection is illustrated schematically in Figure1. Firstly, we need to obtain the residuals of the multivariate condition monitoring data. The test data were obtained from the Mechanical Diagnostic Test Bed (MDTB) in Pennsylvania State University, as shown in Figure 2.Figure 1
The flow chart of Bayesian control scheme.Figure 2
Mechanical diagnostic test bed.Gearbox used in this test includes a pinion gear with 21 teeth and a drive gear with 70 teeth. The gearbox was driven by induction motors (22.38 kW, 1750 rpm) and the torques were provided by alternating current absorption motors (55.95 kW, 1750 rpm). The power of the gearbox is 3.73~14.92 kW and the speed ratio range is 1.2 : 1~6. #2 and #3 accelerometers were installed on the gearbox body, respectively, in the axial and horizontal directions, as shown in Figure3.Figure 3
Locations of #2 and #3 accelerometers.
(a)
Side view
(b)
Front viewSamples were taken every 8 minutes (0.1333 h) and stored in a new file with a sampling frequency of 20 kHz. Each sampling time is 10 s. The 16-bit resolution AD converter was used to ensure the accelerometers have sufficient precision. Driving motor speed V01 and gearbox output torque V05 were also collected with sampling frequency of 1 kHz and sampling time of 10 s. The test is divided into two stages: at a certain rotating speed, gearbox firstly was run at 100% rated output torque for 96 hours (555 in-lbs); it was then run under varying load until failure. In the second stage, the output torque was periodically increased from 50% to 300% rated torque (50%, 100%, 150%, 200%, 250%, and 300%) and then dropped from 300% to 50% rated torque (300%, 250%, 200%, 150%, 100%, and 50%) gradually. The output torque of the gearbox is shown in Figure4. Figure 5 shows the mean value of the input shaft speed. The input shaft speed fluctuation error in Figure 5 is within 0.06%; therefore the speed can be considered as a constant.Figure 4
Output torque of the gearbox.Figure 5
Input speed of the driving motor.The gearbox was run normally in the first stage for 96 hours; thus only the vibration data collected in the second stage (varying load) were analyzed. Sampling points collected by #2 and #3 sensors were, respectively, saved in files A02 and A03, and each of them contained 145 data files under varying load. The original vibration data is shown in Figure6. The results of shutdown fault inspection were shown in Figure 7 that the pinion gear was normal while five fully broken and two partially broken teeth were found in the drive gear.Figure 6
Original data of gearbox under varying load.Figure 7
Test results of the drive gear.
(a)
Broken drive gear
(b)
Five fully broken and two partially broken teeth


### 2.2. Gear Motion Residual (GMR) Signal

TSA algorithm, which extracts the meshing frequency components from the gear vibration data and then synchronously adds and averages them with the rotation of the gear shaft, is widely used in fault diagnosis of gears. The number of sampling points in a single rotational period can be calculated using the following formula:(1)K=fsfm/Ng,where fs is sampling frequency, fm is the fundamental meshing frequency of the gear, Ng is the number of teeth of the gear, and · is the ceil function.Suppose thatV(k), k=1,2,…,n data points were contained in each sampling data file, the rotation period number M of the gear is given by(2)M=nK,where · represents the floor function.In a complete revolution of the gear of interest, the TSA signal is given by(3)Vnk=1M∑i=0M-1Vk+iK,k=1,2,…,K.The TSA signal of the helicopter gearbox is shown in Figure8.Figure 8
TSA signal of the original data.It is efficacious to use the TSA method to reduce the influence of the background noise and the nongear vibration source. GMR signal, which can be obtained by removing the gear meshing frequency and its harmonic from the TSA signals, is highly insensitive to the varying load conditions [10]. The GMR signal can be expressed as follows:(4)Rn=Vn-En,where En is the signal composed of the eliminated components.The GMR signal of the gear of interest is shown in Figure9. It can be seen from the comparison between Figure 9 and the TSA signal that the GMR signal effectively removes the noise, representing more obvious signal change features. Further, we consider a VAR model for the healthy part of the GMR signal to obtain the residual.Figure 9
GMR signal of the original data.


### 2.3. GMR Residuals Computation Using VAR Model

In this section, the healthy portion of the GMR signal is used to fit a VAR model. Two-dimensional data can be expressed asZ1i,Z2i,…,Ztii, i=1,2. Assume that the healthy portion of the data obeys a stationary VAR process:(5)Zn=μ+∑r=1pΦrZn-r-μ+εn,n∈Z,where μ∈R2 is the mean value; p∈N is the model order; Φr∈R2×2 is the autocorrelation matrices; εn are independent identically distributed(i.i.d) and obey N20,Σ; Σ∈R2×2 is the covariance. Let δ=μ-∑r=1pΦrμ, so (5) can be expressed with the following form:(6)W=ΦA+B,where(7)W=Zt22,Zt2-12,…,Zp+22,Zp+12,Zt11,Zt1-11,…,Zp+21,Zp+11′,A=δ,Φ1,Φ2,…,Φp′,B=εt22,εt2-12,…,εp+22,εp+12,εt11,εt1-11,…,εp+21,εp+11′,Φ=1111Zt2-12Zp2Zt1-11Zp1Zt2-22Zp-12Zt1-21Zp-11⋮⋮⋮⋮Zt2-p2Z12Zt1-p1Z11′.Reinsel [22] indicated that A and Σ can be obtained using least squares estimates, and the AIC criterion is selected in this paper to determine the model order p. Thus the estimated VAR model parameters γ^=p^,δ^,Φ^1,Φ^2,…,Φ^p can be utilized to calculate the residual Yn of the TSA signals:(8)Yn=Zn-Eγ^Zn∣Z⃑n-1,where Z⃑n-1=Z1,Z2,…,Zn-1.Forn>p^,(9)Yn=Zn-δ^+∑r=1p^Φ^rZn-r.Forn≤p^, the residual Yn can be computed recursively by Kalman filter, for which the details can be found in the Appendix.Using the above method, we obtain the GMR residual signal of the gear, which is shown in Figure10. To compress the huge amount of data, we select the standard deviation, which is commonly used characteristic parameter in the vibration signal processing, to process the residual of the gear, as shown in Figure 11.Figure 10
Residuals for data files 1–145.Figure 11
Standard deviation indicator of residuals.Figure11 illustrates that the residual standard deviation values of A02 and A03 are relatively stable from files 1 to 100; therefore, we assume the system is in the healthy state in this portion, while from files 101 to 145 the system operates in the unhealthy state as the standard deviation values of this part show an obvious increase compared with those of the healthy portion. Normality and independence tests are taken for these two parts of data, the results of which are shown in Table 1.Table 1
p value of the normality and independence tests.
RMS of residuals
Healthy portion (files 1–100)
Unhealthy portion (files 101–145)
Normality
0.1641
0.4032
Independence
0.3321
0.1861From Table1, both the standard deviation values for A02 and A03 passed the normality and independence tests, indicating that standard deviation values of the residual are independent and obey the multivariate normal distribution, which is consistent with the assumption that the observation of the HSMM described in Section 3 obeys multivariate normal distribution.


## 2.1. Experimental Scheme

The overall process of the optimal Bayesian control scheme for the helicopter gearbox early fault detection is illustrated schematically in Figure1. Firstly, we need to obtain the residuals of the multivariate condition monitoring data. The test data were obtained from the Mechanical Diagnostic Test Bed (MDTB) in Pennsylvania State University, as shown in Figure 2.Figure 1
The flow chart of Bayesian control scheme.Figure 2
Mechanical diagnostic test bed.Gearbox used in this test includes a pinion gear with 21 teeth and a drive gear with 70 teeth. The gearbox was driven by induction motors (22.38 kW, 1750 rpm) and the torques were provided by alternating current absorption motors (55.95 kW, 1750 rpm). The power of the gearbox is 3.73~14.92 kW and the speed ratio range is 1.2 : 1~6. #2 and #3 accelerometers were installed on the gearbox body, respectively, in the axial and horizontal directions, as shown in Figure3.Figure 3
Locations of #2 and #3 accelerometers.
(a)
Side view
(b)
Front viewSamples were taken every 8 minutes (0.1333 h) and stored in a new file with a sampling frequency of 20 kHz. Each sampling time is 10 s. The 16-bit resolution AD converter was used to ensure the accelerometers have sufficient precision. Driving motor speed V01 and gearbox output torque V05 were also collected with sampling frequency of 1 kHz and sampling time of 10 s. The test is divided into two stages: at a certain rotating speed, gearbox firstly was run at 100% rated output torque for 96 hours (555 in-lbs); it was then run under varying load until failure. In the second stage, the output torque was periodically increased from 50% to 300% rated torque (50%, 100%, 150%, 200%, 250%, and 300%) and then dropped from 300% to 50% rated torque (300%, 250%, 200%, 150%, 100%, and 50%) gradually. The output torque of the gearbox is shown in Figure4. Figure 5 shows the mean value of the input shaft speed. The input shaft speed fluctuation error in Figure 5 is within 0.06%; therefore the speed can be considered as a constant.Figure 4
Output torque of the gearbox.Figure 5
Input speed of the driving motor.The gearbox was run normally in the first stage for 96 hours; thus only the vibration data collected in the second stage (varying load) were analyzed. Sampling points collected by #2 and #3 sensors were, respectively, saved in files A02 and A03, and each of them contained 145 data files under varying load. The original vibration data is shown in Figure6. The results of shutdown fault inspection were shown in Figure 7 that the pinion gear was normal while five fully broken and two partially broken teeth were found in the drive gear.Figure 6
Original data of gearbox under varying load.Figure 7
Test results of the drive gear.
(a)
Broken drive gear
(b)
Five fully broken and two partially broken teeth


## 2.2. Gear Motion Residual (GMR) Signal

TSA algorithm, which extracts the meshing frequency components from the gear vibration data and then synchronously adds and averages them with the rotation of the gear shaft, is widely used in fault diagnosis of gears. The number of sampling points in a single rotational period can be calculated using the following formula:(1)K=fsfm/Ng,where fs is sampling frequency, fm is the fundamental meshing frequency of the gear, Ng is the number of teeth of the gear, and · is the ceil function.Suppose thatV(k), k=1,2,…,n data points were contained in each sampling data file, the rotation period number M of the gear is given by(2)M=nK,where · represents the floor function.In a complete revolution of the gear of interest, the TSA signal is given by(3)Vnk=1M∑i=0M-1Vk+iK,k=1,2,…,K.The TSA signal of the helicopter gearbox is shown in Figure8.Figure 8
TSA signal of the original data.It is efficacious to use the TSA method to reduce the influence of the background noise and the nongear vibration source. GMR signal, which can be obtained by removing the gear meshing frequency and its harmonic from the TSA signals, is highly insensitive to the varying load conditions [10]. The GMR signal can be expressed as follows:(4)Rn=Vn-En,where En is the signal composed of the eliminated components.The GMR signal of the gear of interest is shown in Figure9. It can be seen from the comparison between Figure 9 and the TSA signal that the GMR signal effectively removes the noise, representing more obvious signal change features. Further, we consider a VAR model for the healthy part of the GMR signal to obtain the residual.Figure 9
GMR signal of the original data.


## 2.3. GMR Residuals Computation Using VAR Model

In this section, the healthy portion of the GMR signal is used to fit a VAR model. Two-dimensional data can be expressed asZ1i,Z2i,…,Ztii, i=1,2. Assume that the healthy portion of the data obeys a stationary VAR process:(5)Zn=μ+∑r=1pΦrZn-r-μ+εn,n∈Z,where μ∈R2 is the mean value; p∈N is the model order; Φr∈R2×2 is the autocorrelation matrices; εn are independent identically distributed(i.i.d) and obey N20,Σ; Σ∈R2×2 is the covariance. Let δ=μ-∑r=1pΦrμ, so (5) can be expressed with the following form:(6)W=ΦA+B,where(7)W=Zt22,Zt2-12,…,Zp+22,Zp+12,Zt11,Zt1-11,…,Zp+21,Zp+11′,A=δ,Φ1,Φ2,…,Φp′,B=εt22,εt2-12,…,εp+22,εp+12,εt11,εt1-11,…,εp+21,εp+11′,Φ=1111Zt2-12Zp2Zt1-11Zp1Zt2-22Zp-12Zt1-21Zp-11⋮⋮⋮⋮Zt2-p2Z12Zt1-p1Z11′.Reinsel [22] indicated that A and Σ can be obtained using least squares estimates, and the AIC criterion is selected in this paper to determine the model order p. Thus the estimated VAR model parameters γ^=p^,δ^,Φ^1,Φ^2,…,Φ^p can be utilized to calculate the residual Yn of the TSA signals:(8)Yn=Zn-Eγ^Zn∣Z⃑n-1,where Z⃑n-1=Z1,Z2,…,Zn-1.Forn>p^,(9)Yn=Zn-δ^+∑r=1p^Φ^rZn-r.Forn≤p^, the residual Yn can be computed recursively by Kalman filter, for which the details can be found in the Appendix.Using the above method, we obtain the GMR residual signal of the gear, which is shown in Figure10. To compress the huge amount of data, we select the standard deviation, which is commonly used characteristic parameter in the vibration signal processing, to process the residual of the gear, as shown in Figure 11.Figure 10
Residuals for data files 1–145.Figure 11
Standard deviation indicator of residuals.Figure11 illustrates that the residual standard deviation values of A02 and A03 are relatively stable from files 1 to 100; therefore, we assume the system is in the healthy state in this portion, while from files 101 to 145 the system operates in the unhealthy state as the standard deviation values of this part show an obvious increase compared with those of the healthy portion. Normality and independence tests are taken for these two parts of data, the results of which are shown in Table 1.Table 1
p value of the normality and independence tests.
RMS of residuals
Healthy portion (files 1–100)
Unhealthy portion (files 101–145)
Normality
0.1641
0.4032
Independence
0.3321
0.1861From Table1, both the standard deviation values for A02 and A03 passed the normality and independence tests, indicating that standard deviation values of the residual are independent and obey the multivariate normal distribution, which is consistent with the assumption that the observation of the HSMM described in Section 3 obeys multivariate normal distribution.


## 3. Residual Life Prediction Based on Hidden Semi-Markov Model

### 3.1. Hidden Semi-Markov Model

It has been demonstrated from in [13, 19, 23] that a 3-state Markov model is adequate to model the deterioration process of the system. We assume that the degradation process of the system is described by a 3-state HSMM, the state space of which is S=0,1,2. States 0 and 1 are unobservable and represent the good and warning system states, respectively. Only state 2 is assumed to be observable and represents the failure state. Suppose the system starts in the healthy state, that is, PX0=0=1. While the sampling interval is fixed, the observation of the system Y=Yk:k∈N is independent given the state of the gear system. yΔ,y2Δ,…,ykΔ∈Rd is used to denote the d-dimension observation vector at each sampling epoch; thus ykΔ given the state x obeys d-dimension multivariate normal distribution with Ndμx,Σx, x=0,1, and the probability density function is given by(10)fykΔ∣XnΔykΔ∣x=12πddet⁡Σxexp⁡-12ykΔ-μx′Σx-1ykΔ-μx,where μ0,μ1∈R2, Σ0,Σ1∈R2×2 are unknown observation parameters.The state transition of the gear system is illustrated in Figure12. The transition probability matrix is P=Piji,j∈S, where Pij represents the probability that the system transits to state j after it leaves from state i.Figure 12
Schematic for system state transition of the 3-state HSMM.Generally, the transition probability matrix of the system is assumed as follows:(11)P=p00p01p02p10p11p12p20p21p22=0p01p02001001.Eq. (11) indicates that the system starts in the healthy state, where p01+p02=1.Suppose that in statei for i=0,1, the sojourn time obeys the general Erlang distribution, and the probability density function is given by(12)fit=λkitki-1e-λtki-1!fort≥0,where ki∈N+ is the unknown shape parameter and λ>0 is unknown rate parameter. Erlang distribution can also be modeled as a series of a given number of exponential phases running one after another until the end of the sojourn time. In order to record the phases in the Markov model, we need to enlarge the state space. Let Θ=K1,K2,K3 be the new state space, where K1=1,…,k1 represents the set of healthy states, K2=k1+1,…,k1+k2 represents the set of unhealthy states, and K3=k1+k2+1 represents the failure state.The system’s failure time is denoted byξ=inf⁡t≥0:Xt=2. Let O denote the data histories and Λ=p01,p02,k1,k2,λ and Ψ=μ0,μ1,Σ0,Σ1 denote the sets of state and observation parameters to be estimated, respectively. Since the sample path (Xt:t≥0) of the state process in HSMM is unobservable, it is difficult to determine the analytical expression of the maximum likelihood function. The expectation-maximization (EM) algorithm can solve it by iteratively maximizing the pseudo likelihood function. Let Λ0=p^01,p^02,k^1,k^2,λ^ and Ψ0=μ^0,μ^1,Σ^0,Σ^1 be the initial values; the EM algorithm is shown in Figure 13.Figure 13
Flow chart of the EM algorithm.Updated parametersΛ∗,Ψ∗ in each step are then used as new initial values to the E-step, leading the iteration in E-step and M-step until the Euclidean norm satisfies Λ∗,Ψ∗-Λ^,Ψ^<ε, for small given ε>0. Detailed expressions of likelihood function and pseudo likelihood function can be found in Khaleghei and Makis [23]. The state parameters p01,p02 can be updated by calculating the unique stationary point from (13)∂QΛ,Ψ∣Λ^,Ψ^∂p01=∂QΛ,Ψ∣Λ^,Ψ^∂p02=0.The explicit expressions of updatedp01∗ in each step are given by(14)p01∗=∑i=1Nα^1i+∑j=1Mβ^1j∑i=1Nα^1+ω^i+∑j=1Mβ^1j.The updatedk1∗,k2∗,λ∗ in each step can be obtained by (15)maxk1,k2,λ⁡∑i=1Nα^1+ω^i+∑j=1Mβ^1+γ^1jk1ln⁡λ+∑i=1Nα^1ik2ln⁡λ+∑i=1Nα^3ik2-1+∑j=1Ma4,g^jf^j+∑i=1Nα^2+ω^ln⁡ti+∑j=1Mβ^3+γ^3jk1-1-∑i=1Nα^1i+∑j=1Mβ^1jln⁡k2-1!-∑i=1Nα^1+ω^i+∑j=1Mβ^1+γ^1jln⁡k1-1!-∑i=1Ntα^1+ω^i+∑j=1Mβ^2+γ^2jλ.For observation parameters, the unique stationary point can be obtained by(16)∂QΛ,Ψ∣Λ^,Ψ^∂μ0=0,∂QΛ,Ψ∣Λ^,Ψ^∂μ1=0,∂QΛ,Ψ∣Λ^,Ψ^∂Σ0-1=0,∂QΛ,Ψ∣Λ^,Ψ^∂Σ1-1=0.Thus the explicit expressions of the updatedΨ∗=μ0∗,μ1∗,Σ0∗,Σ1∗ in each step are given by(17)μ0∗=∑i=1Nn1,C^+∑j=1Mn1,D^∑i=1Ns1,C^+∑j=1Ms1,D^,Σ0∗=∑i=1Nn3,C^+∑j=1Mn3,D^∑i=1Ns1,C^+∑j=1Ms1,D^,μ1∗=∑i=1Nn2,C^+∑j=1Mn2,D^∑i=1Ns2,C^+∑j=1Ms2,D^,Σ1∗=∑i=1Nn4,C^+∑j=1Mn4,D^∑i=1Ns2,C^+∑j=1Ms2,D^,where the form a,b=a′b, n1=0,y1,∑n=12yn,…,∑n=1Tyn′, s1=0,1,…,T, n2=∑n=1Tyn,∑n=2Tyn,…,yT,0′, s2=T,…,1,0′, n3=0,y1-μ^1′y1-μ^1,…,∑n=1Ty1-μ^1′y1-μ^1′, and n4=∑n=1Tyn-μ^2′yn-μ^2,∑n=2Tyn-μ^2′yn-μ^2,…,yT-μ^2′yT-μ^2,0′.The initial values toΛ^ and Ψ^ were assigned and EM algorithm was employed to solve the parameters to be estimated, the convergence criterion for which is Λ∗,Ψ∗-Λ^,Ψ^<1e-6. We assigned different initial values and the results were very similar, which are shown in Tables 2 and 3.Table 2
Model parameters estimation with EM algorithm.
Parameters
Initial values
First iteration
Second iteration
Last iteration
p
^
01
0.8
0.6897
0.9187
1.0000
p
^
02
0.2
0.3103
0.0813
7.6502
e
-
05
k
^
1
3
2
2
2
k
^
2
3
1
2
2
λ
^
0.5
0.4063
0.2918
0.2069
μ
^
0
10
15
13.6567
19.3175
13.9976
19.3198
15.9207
19.4560
μ
^
1
20
30
25.0970
33.5655
27.0138
36.3763
29.9528
38.8550
Σ
^
0
5
8
8
10
5.2117
5.9879
5.9879
9.7890
5.2980
6.3219
6.3219
9.7924
5.3835
6.6110
6.6110
9.2166
Σ
^
1
50
100
100
100
89.6879
138.6759
138.6759
200.8768
102.8798
139.0890
139.0890
203.3276
123.1435
159.1407
159.1407
212.8811
Q
×
1
0
3
—
−2.48
−1.89
−1.45Table 3
Model parameters estimation with EM algorithm using different initial values.
Parameters
Initial values
First iteration
Second iteration
Last iteration
p
^
01
0.5
0.3176
0.8920
1.0000
p
^
02
0.5
0.6824
0.1080
2.2018
e
-
08
k
^
1
5
3
3
2
k
^
2
5
2
2
2
λ
^
0.3
0.2663
0.2199
0.2070
μ
^
0
15
15
13.8769
16.7687
14.8768
18.9879
15.9214
19.4561
μ
^
1
25
35
28.1879
36.9768
29.0180
37.9671
29.9731
38.8583
Σ
^
0
5
10
10
5
3.9112
7.0879
7.0879
6.0788
5.1939
6.9078
6.9078
8.6892
5.4165
6.6329
6.6329
9.1875
Σ
^
1
100
150
150
200
113.6877
152.3243
152.3243
209.7609
119.0451
157.9880
157.9880
210.6791
122.8959
159.2142
159.2142
212.7090
Q
×
1
0
3
—
−1.98
−1.53
−1.32


### 3.2. Remaining Useful Life (RUL) Prediction

LetΠki be the posterior probability that given the observation data yΔ,y2Δ,…,ykΔ∈R2 the system is in state i1≤i≤k1+k2 at sampling time kΔ:(18)Πki=PΘkΔ=i∣ξ>kΔ,yΔ,y2Δ,…,ykΔ,Π→k-1,where ∑i=1k1+k2Πki=1, Π→k=Πk1,Πk2,…,Πkk1+k2 denotes the posterior probability vector. Suppose the system always starts in the first phase of the healthy state, that is, Π0(1)=1. According to the Bayes’ theorem, the posterior probability Πki can be updated iteratively at each sampling epoch by the following formula:(19)Πki=gykΔ∣ΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1×PΘkΔ=i∣ξ>kΔ,yk-1Δ,Π→k-1∑j=1k1+k2gykΔ∣ΘkΔ=j,ξ>kΔ,yk-1Δ,Π→k-1×PΘkΔ=j∣ξ>kΔ,yk-1Δ,Π→k-1,where PΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1=∑1≤m≤iPmiΔΠk-1m/∑m=1k1+k2∑j=mk1+k2PmjΔΠk-1m.It can be obtained from (10) that ykΔ∣ΘkΔ~Ndμx,Σx for x=1,2; thus(20)gykΔ∣ΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1=fykΔ∣μ1,Σ1ifi∈K1fykΔ∣μ2,Σ2ifi∈K2,where fykΔ∣μx,Σx=1/2πddet⁡Σxexp⁡-1/2ykΔ-μx′Σx-1ykΔ-μx,x=1,2.For the gearbox degeneration system, the estimated shape parameters of Erlang distribution in the hidden states arek1=k2=2 (see Table 2); therefore the transition probability matrix PijΔ can be obtained by solving Kolmogorov backward differential equations:(21)Pijt=e-λtλte-λtp01λt22!e-λtp01λt33!e-λt1-e-λt1+λt+p01λt22!+p01λt33!0e-λtp01λte-λtp01λt22!e-λt1-e-λt1+p01λt+p01λt22!00e-λtλte-λt1-e-λt1+λt000e-λt1-e-λt00001.Suppose at sampling timekΔ, the gear system has not failed, that is, ξ>kΔ; then for t≥0 and k1=k2=2, the conditional reliability function of RUL of the gear system is(22)Rt∣Π→k=Pξ-kΔ>t∣ξ>kΔ,Π→k=PΘkΔ+t∉K3∣ξ>kΔ,Π→k=Πk1e-λt1+λt+p01λt22!+p01λt33!+Πk2e-λt1+p01λt+p01λt22!+Πk3e-λt1+λt+Πk4e-λt.Thus the PDF is:(23)ft∣Π→k=d1-Rt∣Π→kdt.Given the model parametersΛ=p01,p02,k1,k2,λ and Ψ=μ0,μ1,Σ0,Σ1, (22) and (23) can be used to update the conditional reliability function and the PDF using the condition monitoring data obtained from each sampling epoch. Taking file 88 to file 93 for analysis, the corresponding conditional reliability functions are shown in Figure 14 and the probability density functions in Figure 15, where “∗” represents the actual residual values corresponding to each sampling file. As can be seen from Figure 15, from file 89, the RUL distribution of the gear system is highly concentrated; the main cause is that with the steady accumulation of condition monitoring data, the uncertainty of the remaining distribution lifetime gradually decreases while the accuracy of the RUL prediction continually improves.Figure 14
Conditional reliability function in files 88~93.Figure 15
Probability density function in files 88~93.The MRL of the gear system is given by:(24)μkΔ=Eξ-kΔ>t∣ξ>kΔ,Π→k=1λ2Πk11+p01+Πk21+2p01+2Πk3+Πk4.The data from file 1 to file 95 are selected for MRL prediction. After each sampling completed, the life prediction is performed when the new observations are available, and then the predicted value is compared with the values that are obtained using HMM. Meanwhile, as illustrated in Table4, the relative error (RE) between actual remaining useful lives and the prediction results using different models are calculated. As can be seen in Table 4, with the increase of collected data, the RE in the prediction results using HSMM are smaller than those using HMM, which is closer to the actual values.Table 4
Comparison of RUL prediction and RE using different models.
File number
ARL
RUL
RE
HSMM
HMM
HSMM
HMM
1
19.34
19.801
19.800
0.024
0.024
2
19.206
19.013
19.942
0.010
0.038
3
19.071
18.203
19.801
0.046
0.038
4
18.937
18.019
19.689
0.049
0.040
5
18.803
16.981
19.332
0.097
0.028
6
18.669
15.002
19.340
0.196
0.036
7
18.535
13.200
19.185
0.288
0.035
⋮
⋮
⋮
⋮
⋮
⋮
88
7.662
8.013
18.343
0.045
1.394
89
7.527
8.003
18.271
0.063
1.427
90
7.393
7.830
18.001
0.059
1.435
91
7.259
7.482
7.701
0.031
0.061
92
7.125
7.450
7.645
0.045
0.073
93
6.990
7.308
7.632
0.046
0.091
94
6.856
7.225
7.545
0.054
0.100
95
6.722
7.031
7.416
0.046
0.103


## 3.1. Hidden Semi-Markov Model

It has been demonstrated from in [13, 19, 23] that a 3-state Markov model is adequate to model the deterioration process of the system. We assume that the degradation process of the system is described by a 3-state HSMM, the state space of which is S=0,1,2. States 0 and 1 are unobservable and represent the good and warning system states, respectively. Only state 2 is assumed to be observable and represents the failure state. Suppose the system starts in the healthy state, that is, PX0=0=1. While the sampling interval is fixed, the observation of the system Y=Yk:k∈N is independent given the state of the gear system. yΔ,y2Δ,…,ykΔ∈Rd is used to denote the d-dimension observation vector at each sampling epoch; thus ykΔ given the state x obeys d-dimension multivariate normal distribution with Ndμx,Σx, x=0,1, and the probability density function is given by(10)fykΔ∣XnΔykΔ∣x=12πddet⁡Σxexp⁡-12ykΔ-μx′Σx-1ykΔ-μx,where μ0,μ1∈R2, Σ0,Σ1∈R2×2 are unknown observation parameters.The state transition of the gear system is illustrated in Figure12. The transition probability matrix is P=Piji,j∈S, where Pij represents the probability that the system transits to state j after it leaves from state i.Figure 12
Schematic for system state transition of the 3-state HSMM.Generally, the transition probability matrix of the system is assumed as follows:(11)P=p00p01p02p10p11p12p20p21p22=0p01p02001001.Eq. (11) indicates that the system starts in the healthy state, where p01+p02=1.Suppose that in statei for i=0,1, the sojourn time obeys the general Erlang distribution, and the probability density function is given by(12)fit=λkitki-1e-λtki-1!fort≥0,where ki∈N+ is the unknown shape parameter and λ>0 is unknown rate parameter. Erlang distribution can also be modeled as a series of a given number of exponential phases running one after another until the end of the sojourn time. In order to record the phases in the Markov model, we need to enlarge the state space. Let Θ=K1,K2,K3 be the new state space, where K1=1,…,k1 represents the set of healthy states, K2=k1+1,…,k1+k2 represents the set of unhealthy states, and K3=k1+k2+1 represents the failure state.The system’s failure time is denoted byξ=inf⁡t≥0:Xt=2. Let O denote the data histories and Λ=p01,p02,k1,k2,λ and Ψ=μ0,μ1,Σ0,Σ1 denote the sets of state and observation parameters to be estimated, respectively. Since the sample path (Xt:t≥0) of the state process in HSMM is unobservable, it is difficult to determine the analytical expression of the maximum likelihood function. The expectation-maximization (EM) algorithm can solve it by iteratively maximizing the pseudo likelihood function. Let Λ0=p^01,p^02,k^1,k^2,λ^ and Ψ0=μ^0,μ^1,Σ^0,Σ^1 be the initial values; the EM algorithm is shown in Figure 13.Figure 13
Flow chart of the EM algorithm.Updated parametersΛ∗,Ψ∗ in each step are then used as new initial values to the E-step, leading the iteration in E-step and M-step until the Euclidean norm satisfies Λ∗,Ψ∗-Λ^,Ψ^<ε, for small given ε>0. Detailed expressions of likelihood function and pseudo likelihood function can be found in Khaleghei and Makis [23]. The state parameters p01,p02 can be updated by calculating the unique stationary point from (13)∂QΛ,Ψ∣Λ^,Ψ^∂p01=∂QΛ,Ψ∣Λ^,Ψ^∂p02=0.The explicit expressions of updatedp01∗ in each step are given by(14)p01∗=∑i=1Nα^1i+∑j=1Mβ^1j∑i=1Nα^1+ω^i+∑j=1Mβ^1j.The updatedk1∗,k2∗,λ∗ in each step can be obtained by (15)maxk1,k2,λ⁡∑i=1Nα^1+ω^i+∑j=1Mβ^1+γ^1jk1ln⁡λ+∑i=1Nα^1ik2ln⁡λ+∑i=1Nα^3ik2-1+∑j=1Ma4,g^jf^j+∑i=1Nα^2+ω^ln⁡ti+∑j=1Mβ^3+γ^3jk1-1-∑i=1Nα^1i+∑j=1Mβ^1jln⁡k2-1!-∑i=1Nα^1+ω^i+∑j=1Mβ^1+γ^1jln⁡k1-1!-∑i=1Ntα^1+ω^i+∑j=1Mβ^2+γ^2jλ.For observation parameters, the unique stationary point can be obtained by(16)∂QΛ,Ψ∣Λ^,Ψ^∂μ0=0,∂QΛ,Ψ∣Λ^,Ψ^∂μ1=0,∂QΛ,Ψ∣Λ^,Ψ^∂Σ0-1=0,∂QΛ,Ψ∣Λ^,Ψ^∂Σ1-1=0.Thus the explicit expressions of the updatedΨ∗=μ0∗,μ1∗,Σ0∗,Σ1∗ in each step are given by(17)μ0∗=∑i=1Nn1,C^+∑j=1Mn1,D^∑i=1Ns1,C^+∑j=1Ms1,D^,Σ0∗=∑i=1Nn3,C^+∑j=1Mn3,D^∑i=1Ns1,C^+∑j=1Ms1,D^,μ1∗=∑i=1Nn2,C^+∑j=1Mn2,D^∑i=1Ns2,C^+∑j=1Ms2,D^,Σ1∗=∑i=1Nn4,C^+∑j=1Mn4,D^∑i=1Ns2,C^+∑j=1Ms2,D^,where the form a,b=a′b, n1=0,y1,∑n=12yn,…,∑n=1Tyn′, s1=0,1,…,T, n2=∑n=1Tyn,∑n=2Tyn,…,yT,0′, s2=T,…,1,0′, n3=0,y1-μ^1′y1-μ^1,…,∑n=1Ty1-μ^1′y1-μ^1′, and n4=∑n=1Tyn-μ^2′yn-μ^2,∑n=2Tyn-μ^2′yn-μ^2,…,yT-μ^2′yT-μ^2,0′.The initial values toΛ^ and Ψ^ were assigned and EM algorithm was employed to solve the parameters to be estimated, the convergence criterion for which is Λ∗,Ψ∗-Λ^,Ψ^<1e-6. We assigned different initial values and the results were very similar, which are shown in Tables 2 and 3.Table 2
Model parameters estimation with EM algorithm.
Parameters
Initial values
First iteration
Second iteration
Last iteration
p
^
01
0.8
0.6897
0.9187
1.0000
p
^
02
0.2
0.3103
0.0813
7.6502
e
-
05
k
^
1
3
2
2
2
k
^
2
3
1
2
2
λ
^
0.5
0.4063
0.2918
0.2069
μ
^
0
10
15
13.6567
19.3175
13.9976
19.3198
15.9207
19.4560
μ
^
1
20
30
25.0970
33.5655
27.0138
36.3763
29.9528
38.8550
Σ
^
0
5
8
8
10
5.2117
5.9879
5.9879
9.7890
5.2980
6.3219
6.3219
9.7924
5.3835
6.6110
6.6110
9.2166
Σ
^
1
50
100
100
100
89.6879
138.6759
138.6759
200.8768
102.8798
139.0890
139.0890
203.3276
123.1435
159.1407
159.1407
212.8811
Q
×
1
0
3
—
−2.48
−1.89
−1.45Table 3
Model parameters estimation with EM algorithm using different initial values.
Parameters
Initial values
First iteration
Second iteration
Last iteration
p
^
01
0.5
0.3176
0.8920
1.0000
p
^
02
0.5
0.6824
0.1080
2.2018
e
-
08
k
^
1
5
3
3
2
k
^
2
5
2
2
2
λ
^
0.3
0.2663
0.2199
0.2070
μ
^
0
15
15
13.8769
16.7687
14.8768
18.9879
15.9214
19.4561
μ
^
1
25
35
28.1879
36.9768
29.0180
37.9671
29.9731
38.8583
Σ
^
0
5
10
10
5
3.9112
7.0879
7.0879
6.0788
5.1939
6.9078
6.9078
8.6892
5.4165
6.6329
6.6329
9.1875
Σ
^
1
100
150
150
200
113.6877
152.3243
152.3243
209.7609
119.0451
157.9880
157.9880
210.6791
122.8959
159.2142
159.2142
212.7090
Q
×
1
0
3
—
−1.98
−1.53
−1.32


## 3.2. Remaining Useful Life (RUL) Prediction

LetΠki be the posterior probability that given the observation data yΔ,y2Δ,…,ykΔ∈R2 the system is in state i1≤i≤k1+k2 at sampling time kΔ:(18)Πki=PΘkΔ=i∣ξ>kΔ,yΔ,y2Δ,…,ykΔ,Π→k-1,where ∑i=1k1+k2Πki=1, Π→k=Πk1,Πk2,…,Πkk1+k2 denotes the posterior probability vector. Suppose the system always starts in the first phase of the healthy state, that is, Π0(1)=1. According to the Bayes’ theorem, the posterior probability Πki can be updated iteratively at each sampling epoch by the following formula:(19)Πki=gykΔ∣ΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1×PΘkΔ=i∣ξ>kΔ,yk-1Δ,Π→k-1∑j=1k1+k2gykΔ∣ΘkΔ=j,ξ>kΔ,yk-1Δ,Π→k-1×PΘkΔ=j∣ξ>kΔ,yk-1Δ,Π→k-1,where PΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1=∑1≤m≤iPmiΔΠk-1m/∑m=1k1+k2∑j=mk1+k2PmjΔΠk-1m.It can be obtained from (10) that ykΔ∣ΘkΔ~Ndμx,Σx for x=1,2; thus(20)gykΔ∣ΘkΔ=i,ξ>kΔ,yk-1Δ,Π→k-1=fykΔ∣μ1,Σ1ifi∈K1fykΔ∣μ2,Σ2ifi∈K2,where fykΔ∣μx,Σx=1/2πddet⁡Σxexp⁡-1/2ykΔ-μx′Σx-1ykΔ-μx,x=1,2.For the gearbox degeneration system, the estimated shape parameters of Erlang distribution in the hidden states arek1=k2=2 (see Table 2); therefore the transition probability matrix PijΔ can be obtained by solving Kolmogorov backward differential equations:(21)Pijt=e-λtλte-λtp01λt22!e-λtp01λt33!e-λt1-e-λt1+λt+p01λt22!+p01λt33!0e-λtp01λte-λtp01λt22!e-λt1-e-λt1+p01λt+p01λt22!00e-λtλte-λt1-e-λt1+λt000e-λt1-e-λt00001.Suppose at sampling timekΔ, the gear system has not failed, that is, ξ>kΔ; then for t≥0 and k1=k2=2, the conditional reliability function of RUL of the gear system is(22)Rt∣Π→k=Pξ-kΔ>t∣ξ>kΔ,Π→k=PΘkΔ+t∉K3∣ξ>kΔ,Π→k=Πk1e-λt1+λt+p01λt22!+p01λt33!+Πk2e-λt1+p01λt+p01λt22!+Πk3e-λt1+λt+Πk4e-λt.Thus the PDF is:(23)ft∣Π→k=d1-Rt∣Π→kdt.Given the model parametersΛ=p01,p02,k1,k2,λ and Ψ=μ0,μ1,Σ0,Σ1, (22) and (23) can be used to update the conditional reliability function and the PDF using the condition monitoring data obtained from each sampling epoch. Taking file 88 to file 93 for analysis, the corresponding conditional reliability functions are shown in Figure 14 and the probability density functions in Figure 15, where “∗” represents the actual residual values corresponding to each sampling file. As can be seen from Figure 15, from file 89, the RUL distribution of the gear system is highly concentrated; the main cause is that with the steady accumulation of condition monitoring data, the uncertainty of the remaining distribution lifetime gradually decreases while the accuracy of the RUL prediction continually improves.Figure 14
Conditional reliability function in files 88~93.Figure 15
Probability density function in files 88~93.The MRL of the gear system is given by:(24)μkΔ=Eξ-kΔ>t∣ξ>kΔ,Π→k=1λ2Πk11+p01+Πk21+2p01+2Πk3+Πk4.The data from file 1 to file 95 are selected for MRL prediction. After each sampling completed, the life prediction is performed when the new observations are available, and then the predicted value is compared with the values that are obtained using HMM. Meanwhile, as illustrated in Table4, the relative error (RE) between actual remaining useful lives and the prediction results using different models are calculated. As can be seen in Table 4, with the increase of collected data, the RE in the prediction results using HSMM are smaller than those using HMM, which is closer to the actual values.Table 4
Comparison of RUL prediction and RE using different models.
File number
ARL
RUL
RE
HSMM
HMM
HSMM
HMM
1
19.34
19.801
19.800
0.024
0.024
2
19.206
19.013
19.942
0.010
0.038
3
19.071
18.203
19.801
0.046
0.038
4
18.937
18.019
19.689
0.049
0.040
5
18.803
16.981
19.332
0.097
0.028
6
18.669
15.002
19.340
0.196
0.036
7
18.535
13.200
19.185
0.288
0.035
⋮
⋮
⋮
⋮
⋮
⋮
88
7.662
8.013
18.343
0.045
1.394
89
7.527
8.003
18.271
0.063
1.427
90
7.393
7.830
18.001
0.059
1.435
91
7.259
7.482
7.701
0.031
0.061
92
7.125
7.450
7.645
0.045
0.073
93
6.990
7.308
7.632
0.046
0.091
94
6.856
7.225
7.545
0.054
0.100
95
6.722
7.031
7.416
0.046
0.103


## 4. Optimal Bayesian Control Scheme

### 4.1. Optimal Bayesian Control Chart Based on HSMM

In this section, we design a multivariate Bayesian control chart to detect the early fault of the helicopter gearbox. The optimal control limit is used to determine the stopping time, of which the objective is to minimize the long-run expected average cost per unit time. In partially observable Markov decision process (POMDP), it is well known that the posterior probability that the system is in the unhealthy state provides enough information for maintenance decision making [13]. Let Πk1 be the posterior probability that the helicopter gear system is in the unhealthy state:(25)Πk1=∑i∈K2Πki=Πk3+Πk4,where the initial value Π01=0.By renewal theory, the cost minimization problem with fixed sampling interval is equivalent to seeking an optimal value ofΠ∗∈0,1, such that(26)gΠ∗=infΠ∈0,1⁡EΠ∗CCEΠ∗CL,where CC and CL denote, respectively, the cycle cost and cycle length.The Bayesian control scheme is illustrated in Figure16. For the fixed control limit Π-∈0,1, at each decision epoch, the posterior probability vector Π→k is updated by Bayes’ theorem. At the sampling time kΔ, if Πk1<Π-, the system will continue to run, of which the cost for each sampling is CS. If Πk1≥Π-, system should be stopped and full inspection is performed to determine whether the system is in state 0 or not, for which the inspection cost rate and inspection time are, respectively, CI and TI. The gear system has probability Πk1 to be in state 1 and probability 1-Πk1 to be in state 0. If the system is found to be in state 0 after a full inspection, then no maintenance measures are taken and the system continues operating in the healthy condition. If the system is found to be in state 1, preventive maintenance measures will be taken with the maintenance cost rate CPM and the maintenance time TPM. If the posterior probability Πk1 is not beyond the fixed control limit Π- but still fails, that is, a random failure occurs, we need to immediately take failure replacement measures with corresponding maintenance cost rate and replacement time CF and TF, respectively. When the system takes full inspection, preventive maintenance, or failure replacement, the lost production cost rate CLP will be incurred. When the system operates in the warning state, additional operating cost rate CAO and maintenance cost rate CAM will be generated. We further assume that after the inspection, the preventive maintenance, or failure replacement, the system will begin from a new cycle.Figure 16
Schematic of Bayesian control approach.We design an effective algorithm in the semi-Markov decision process (SMDP) framework to obtain the optimal control limit. In order to discretize the interval0,1, it is sufficient to choose the number of subintervals L=30 to provide effective accuracy. For L=30, we define the SMDP in state i1,i2,…,ik1+k2-1; that is, if the posterior probability in state 1 lies in the interval i1-1/L,i1/L, the posterior probability in state 2 lies in the interval i2-1/L,i2/L,…, and the posterior probability in state k1+k2-1 lies in the interval ik1+k2-1-1/L,ik1+k2-1/L. If Πk1≥Π-, the preventive maintenance is performed, then the SMDP is defined in state PM. If the system fails, the SMDP is defined in state F. Let S-=i1,i2,…,ik1+k2-1:∑m=1k1+k2-1im≤L,im∈N+, then the state space for the SMDP is defined as S=1,0,…,0∪S-∪PM∪F, where 1,0,…,0 presents the gear system is in state one at the initial moment.For the fixed control limitΠ-, gΠ- can be obtained by solving the following linear equations:(27)vx=cx-gΠ-τx+∑x′∈SPx,x′vx′forx,x′∈S,(28)v0=0,where vx is the relative value until the next decision epoch given the current state x∈S. cx is the expected cost incurred until the next decision epoch given the current state x∈S. τx is the expected sojourn time until the next decision epoch given the current state x∈S. Px,x′ is the probability that at the next decision epoch the system will be in state x′∈S given the current state is x′∈S.Next, we derive the closed form expressions of each quantity in (27). For each i1,i2,…,ik1+k2-1, i1′,i2′,…,ik1+k2-1′∈S-, and ja-=ia′-1/L, ja+=ia′/L for a=1,…,k1+k2-1, the transition probabilities of SMDP are calculated by the following formula:(29)Pi1,i2,…,ik1+k2-1i1′,i2′,…,ik1+k2-1′=Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+,ξ>kΔ∣ξ>k-1Δ,Π→k-1=∑m=1k1+k2Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+∣ξ>kΔ,Π→k-1,ΘkΔ=m×PΘkΔ=m∣Π→k-1×RΔ∣Π→k-1.From (10) we can further obtain:(30)fykΔ∣μ1,Σ1fykΔ∣μ0,Σ0=Σ0·Σ1-11/2exp⁡-12Vk+C,where Vk=YkΔ-B′AYkΔ-B, A=Σ1-1-Σ0-1, B=Σ1-1-Σ0-1-1Σ1-1μ1-Σ0-1μ0, C=μ1TΣ1-1μ1-μ0TΣ0-1μ0-BTΣ1-1μ1-Σ0-1μ0.Thus (18) can be expressed as follows:(31)Πkm=cΠk-1m∑∀j∈K1cΠk-1j+fykΔ∣μ1,Σ1/fykΔ∣μ0,Σ0∑∀j∈K2cΠk-1jifm∈K1cΠn-1mfykΔ∣μ0,Σ0/fykΔ∣μ1,Σ1∑∀j∈K1cΠk-1j+∑∀j∈K2cΠk-1jifm∈K2,where cΠk-1m=∑n=1mΠk-1nPnmΔ.In (29), for m∈K1, Pjm-≤Πk(m)<jm+∣ξ>kΔ,Π→k-1 is computed as follows:(32)Pjm-≤Πkm<jm+∣ξ>kΔ,Π→k-1=∑m=1k1+k2P-2ln⁡cΠk-1m-jm-∑j∈K1cΠk-1jΣ0-1·Σ11/2jm-∑j∈K2cΠk-1j-C≤Vk<-2ln⁡cΠk-1m-jm+∑j∈K1cΠk-1jΣ0-1·Σ11/2jm+∑j∈K2cΠk-1j-C∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j=∑m=1k1+k2PTjm-≤Vk<Tjm+∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j.Form∈K2, the probability Pjm-≤Πk(m)<jm+∣ξ>kΔ,Π→k-1 can be derived using the similar method in (32) and is omitted. So the first item in (29) can be written as follows:(33)Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+∣ξ>kΔ,Π→k-1=∑m=1k1+k2Pmax⁡Tj1-,…,Tjk1+k2-1-≤Vk<min⁡Tj1+,…,Tjk1+k2-1+∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j.SinceYkΔ-B∣ΘkΔ∈K1~Nμ0-B,Σ0, YkΔ-B∣ΘkΔ∈K2~Nμ1-B,Σ1, the probability in (33) can be calculated using Theorem  3.1 proposed by Provost and Rudiuk [24]. Next, the other transition probabilities in SMDP are derived.IfΠk-11≥Π∗, after full inspection, the probability of system in state 1, that is, true alarm occurs, is given by(34)Pi1,i2,…,ik1+k2-1,PM=Πk-11=Πk3+Πk4.Also, after full inspection, the gear system may be found in healthy state; the probability is given by(35)Pi1,i2,…,ik1+k2-1,1,0,…,0=1-Πk-11=1-Πk3-Πk4.IfΠk-11<Π∗, the probability of failure for the gear system is given by(36)Pi1,i2,…,ik1+k2-1,F=1-RΔ∣Π→k-1=1-Πk-11e-λΔ1+λΔ+p01λΔ22!+p01λΔ33!-Πk-12e-λΔ1+p01λΔ+p01λΔ22!-Πk-13e-λΔ1+λΔ-Πk-14e-λΔ.The remainder transition probabilities are given by(37)PF,1,0,…,0=1,PPM,1,0,…,0=1.IfΠk-11≥Π∗, the corresponding expected cost and expected sojourn time are given by(38)ci1,i2,…,ik1+k2-1=CI+CLP·TI,τi1,i2,…,ik1+k2-1=TI.IfΠk-11<Π∗, the expected cost and expected sojourn time can be computed as follows:(39)ci1,i2,…,ik1+k2-1=Cs∑i1′,i2′,…,ik1+k2-1′∈S-Pi1,i2,…,ik1+k2-1,i1′,i2′,…,ik1+k2-1′+CF+CLP·TF·Pi1,i2,…,ik1+k2-1,F+CAO+CAM∫0Δ∑j∈K2∑i=1jΠk-1iPijtdt,τi1,i2,…,ik1+k2-1=Δ∑i1′,i2′,…,ik1+k2-1′∈S-Pi1,i2,…,ik1+k2-1i1′,i2′,…,ik1+k2-1′+∫0ΔTF+tft∣Π→k-1dt.If the gear system is inPM or F state, the corresponding expected costs and expected sojourn times are given by(40)cPM=CP+CLP·TPM,cF=CF+CLP·TF,τPM=TPM,τF=TF.By setting the maintenance time parametersTI=3 h, TPM=15 h, and TF=30 h and maintenance cost parameters CS=20, CI=100, CPM=300, CF=1000, CLP=300, CAO=100, and CAM=200, we coded the computational algorithm in MATLAB R2015b to calculate the different expected average costs under different control limits Π-. The optimal control limit was Π∗=0.71 and the corresponding average cost was equal to 113.59. The Bayesian control chart of the helicopter gear system is shown in Figure 17. It is observed that based on HSMM the posterior probability Πk1 firstly exceeds the optimal control limit at file 89, which indicates that the system need to be stopped to take full inspection at this moment.Figure 17
Optimal multivariate Bayesian control chart of the gearbox.


### 4.2. Comparison with Other Maintenance Policies

Using the same maintenance time and cost parameters, the optimal control limitΠ∗=0.28 was obtained in the Bayesian control chart based on HMM, in which the sojourn time in each hidden states follows exponential distribution. The corresponding average cost was equal to 121.72. It should be noted that the average cost concluded by HSMM was less than that by HMM. As shown in Figure 17, while based on HMM, the posterior probability Πk1 firstly exceeds the optimal control limit at file 91. The multivariate Bayesian control chart indicates that the proposed Bayesian control scheme can detect the fault of the gear of interest much earlier than the method based on HMM.Next, we considered the traditional age-based replacement policy, which is well known that it does not take into account condition monitoring information. For the helicopter gearbox, the expected cycle cost (CC) under age-based replacement policy is the sum of failure cost, additional operating and maintenance costs while the system operates in unhealthy state and the preventive maintenance cost:(41)EτCC=CF+CLPTF1-Rτ∣Π⃑0+CAO+CAM∫0τP13t+P14tdt+CP+CLPTPRτ∣Π⃑0,where τ is the stopping time.The expected cycle length (CL) is the sum of the expected time under normal operation, replacement time caused by failure, and preventive maintenance time:(42)EτCL=∫0τRt∣Π⃑0dt+TF1-Rτ∣Π⃑0+TPRτ∣Π⃑0.Then the long-run expected average cost per unit time under age-based replacement policy is equivalent to finding the optimal stopping timeτ∗, such that(43)τ∗=arg⁡inf⁡EτCCEτCL.Table5 shows the optimal values obtained under different policies, which are, respectively, the multivariate Bayesian control chart based on HSMM and HMM and the age-based replacement policy. It is interesting to find that the expected average cost using the Bayesian control chart based on the HSMM is considerably lower than the other policies, while the cost obtained by the age-based replacement policy is the highest due to the irrespective of condition monitoring information. This is reasonable because, for the helicopter gear system, the cost incurred by failure is always higher than the cost incurred by inspection and preventive maintenance measures.Table 5
The optimal maintenance cost under different maintenance policies.
Policy
Optimal parameters
Expected average cost
Π
∗
τ
∗
Bayesian (HSMM)
0.71
—
113.59
Bayesian (HMM)
0.28
—
121.72
Aged-based replacement
—
17.24
140.23


## 4.1. Optimal Bayesian Control Chart Based on HSMM

In this section, we design a multivariate Bayesian control chart to detect the early fault of the helicopter gearbox. The optimal control limit is used to determine the stopping time, of which the objective is to minimize the long-run expected average cost per unit time. In partially observable Markov decision process (POMDP), it is well known that the posterior probability that the system is in the unhealthy state provides enough information for maintenance decision making [13]. Let Πk1 be the posterior probability that the helicopter gear system is in the unhealthy state:(25)Πk1=∑i∈K2Πki=Πk3+Πk4,where the initial value Π01=0.By renewal theory, the cost minimization problem with fixed sampling interval is equivalent to seeking an optimal value ofΠ∗∈0,1, such that(26)gΠ∗=infΠ∈0,1⁡EΠ∗CCEΠ∗CL,where CC and CL denote, respectively, the cycle cost and cycle length.The Bayesian control scheme is illustrated in Figure16. For the fixed control limit Π-∈0,1, at each decision epoch, the posterior probability vector Π→k is updated by Bayes’ theorem. At the sampling time kΔ, if Πk1<Π-, the system will continue to run, of which the cost for each sampling is CS. If Πk1≥Π-, system should be stopped and full inspection is performed to determine whether the system is in state 0 or not, for which the inspection cost rate and inspection time are, respectively, CI and TI. The gear system has probability Πk1 to be in state 1 and probability 1-Πk1 to be in state 0. If the system is found to be in state 0 after a full inspection, then no maintenance measures are taken and the system continues operating in the healthy condition. If the system is found to be in state 1, preventive maintenance measures will be taken with the maintenance cost rate CPM and the maintenance time TPM. If the posterior probability Πk1 is not beyond the fixed control limit Π- but still fails, that is, a random failure occurs, we need to immediately take failure replacement measures with corresponding maintenance cost rate and replacement time CF and TF, respectively. When the system takes full inspection, preventive maintenance, or failure replacement, the lost production cost rate CLP will be incurred. When the system operates in the warning state, additional operating cost rate CAO and maintenance cost rate CAM will be generated. We further assume that after the inspection, the preventive maintenance, or failure replacement, the system will begin from a new cycle.Figure 16
Schematic of Bayesian control approach.We design an effective algorithm in the semi-Markov decision process (SMDP) framework to obtain the optimal control limit. In order to discretize the interval0,1, it is sufficient to choose the number of subintervals L=30 to provide effective accuracy. For L=30, we define the SMDP in state i1,i2,…,ik1+k2-1; that is, if the posterior probability in state 1 lies in the interval i1-1/L,i1/L, the posterior probability in state 2 lies in the interval i2-1/L,i2/L,…, and the posterior probability in state k1+k2-1 lies in the interval ik1+k2-1-1/L,ik1+k2-1/L. If Πk1≥Π-, the preventive maintenance is performed, then the SMDP is defined in state PM. If the system fails, the SMDP is defined in state F. Let S-=i1,i2,…,ik1+k2-1:∑m=1k1+k2-1im≤L,im∈N+, then the state space for the SMDP is defined as S=1,0,…,0∪S-∪PM∪F, where 1,0,…,0 presents the gear system is in state one at the initial moment.For the fixed control limitΠ-, gΠ- can be obtained by solving the following linear equations:(27)vx=cx-gΠ-τx+∑x′∈SPx,x′vx′forx,x′∈S,(28)v0=0,where vx is the relative value until the next decision epoch given the current state x∈S. cx is the expected cost incurred until the next decision epoch given the current state x∈S. τx is the expected sojourn time until the next decision epoch given the current state x∈S. Px,x′ is the probability that at the next decision epoch the system will be in state x′∈S given the current state is x′∈S.Next, we derive the closed form expressions of each quantity in (27). For each i1,i2,…,ik1+k2-1, i1′,i2′,…,ik1+k2-1′∈S-, and ja-=ia′-1/L, ja+=ia′/L for a=1,…,k1+k2-1, the transition probabilities of SMDP are calculated by the following formula:(29)Pi1,i2,…,ik1+k2-1i1′,i2′,…,ik1+k2-1′=Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+,ξ>kΔ∣ξ>k-1Δ,Π→k-1=∑m=1k1+k2Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+∣ξ>kΔ,Π→k-1,ΘkΔ=m×PΘkΔ=m∣Π→k-1×RΔ∣Π→k-1.From (10) we can further obtain:(30)fykΔ∣μ1,Σ1fykΔ∣μ0,Σ0=Σ0·Σ1-11/2exp⁡-12Vk+C,where Vk=YkΔ-B′AYkΔ-B, A=Σ1-1-Σ0-1, B=Σ1-1-Σ0-1-1Σ1-1μ1-Σ0-1μ0, C=μ1TΣ1-1μ1-μ0TΣ0-1μ0-BTΣ1-1μ1-Σ0-1μ0.Thus (18) can be expressed as follows:(31)Πkm=cΠk-1m∑∀j∈K1cΠk-1j+fykΔ∣μ1,Σ1/fykΔ∣μ0,Σ0∑∀j∈K2cΠk-1jifm∈K1cΠn-1mfykΔ∣μ0,Σ0/fykΔ∣μ1,Σ1∑∀j∈K1cΠk-1j+∑∀j∈K2cΠk-1jifm∈K2,where cΠk-1m=∑n=1mΠk-1nPnmΔ.In (29), for m∈K1, Pjm-≤Πk(m)<jm+∣ξ>kΔ,Π→k-1 is computed as follows:(32)Pjm-≤Πkm<jm+∣ξ>kΔ,Π→k-1=∑m=1k1+k2P-2ln⁡cΠk-1m-jm-∑j∈K1cΠk-1jΣ0-1·Σ11/2jm-∑j∈K2cΠk-1j-C≤Vk<-2ln⁡cΠk-1m-jm+∑j∈K1cΠk-1jΣ0-1·Σ11/2jm+∑j∈K2cΠk-1j-C∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j=∑m=1k1+k2PTjm-≤Vk<Tjm+∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j.Form∈K2, the probability Pjm-≤Πk(m)<jm+∣ξ>kΔ,Π→k-1 can be derived using the similar method in (32) and is omitted. So the first item in (29) can be written as follows:(33)Pj1-≤Πk1<j1+,…,jk1+k2-1-≤Πkk1+k2-1<jk1+k2-1+∣ξ>kΔ,Π→k-1=∑m=1k1+k2Pmax⁡Tj1-,…,Tjk1+k2-1-≤Vk<min⁡Tj1+,…,Tjk1+k2-1+∣ΘkΔ=mcΠk-1m∑j=1k1+k2cΠk-1j.SinceYkΔ-B∣ΘkΔ∈K1~Nμ0-B,Σ0, YkΔ-B∣ΘkΔ∈K2~Nμ1-B,Σ1, the probability in (33) can be calculated using Theorem  3.1 proposed by Provost and Rudiuk [24]. Next, the other transition probabilities in SMDP are derived.IfΠk-11≥Π∗, after full inspection, the probability of system in state 1, that is, true alarm occurs, is given by(34)Pi1,i2,…,ik1+k2-1,PM=Πk-11=Πk3+Πk4.Also, after full inspection, the gear system may be found in healthy state; the probability is given by(35)Pi1,i2,…,ik1+k2-1,1,0,…,0=1-Πk-11=1-Πk3-Πk4.IfΠk-11<Π∗, the probability of failure for the gear system is given by(36)Pi1,i2,…,ik1+k2-1,F=1-RΔ∣Π→k-1=1-Πk-11e-λΔ1+λΔ+p01λΔ22!+p01λΔ33!-Πk-12e-λΔ1+p01λΔ+p01λΔ22!-Πk-13e-λΔ1+λΔ-Πk-14e-λΔ.The remainder transition probabilities are given by(37)PF,1,0,…,0=1,PPM,1,0,…,0=1.IfΠk-11≥Π∗, the corresponding expected cost and expected sojourn time are given by(38)ci1,i2,…,ik1+k2-1=CI+CLP·TI,τi1,i2,…,ik1+k2-1=TI.IfΠk-11<Π∗, the expected cost and expected sojourn time can be computed as follows:(39)ci1,i2,…,ik1+k2-1=Cs∑i1′,i2′,…,ik1+k2-1′∈S-Pi1,i2,…,ik1+k2-1,i1′,i2′,…,ik1+k2-1′+CF+CLP·TF·Pi1,i2,…,ik1+k2-1,F+CAO+CAM∫0Δ∑j∈K2∑i=1jΠk-1iPijtdt,τi1,i2,…,ik1+k2-1=Δ∑i1′,i2′,…,ik1+k2-1′∈S-Pi1,i2,…,ik1+k2-1i1′,i2′,…,ik1+k2-1′+∫0ΔTF+tft∣Π→k-1dt.If the gear system is inPM or F state, the corresponding expected costs and expected sojourn times are given by(40)cPM=CP+CLP·TPM,cF=CF+CLP·TF,τPM=TPM,τF=TF.By setting the maintenance time parametersTI=3 h, TPM=15 h, and TF=30 h and maintenance cost parameters CS=20, CI=100, CPM=300, CF=1000, CLP=300, CAO=100, and CAM=200, we coded the computational algorithm in MATLAB R2015b to calculate the different expected average costs under different control limits Π-. The optimal control limit was Π∗=0.71 and the corresponding average cost was equal to 113.59. The Bayesian control chart of the helicopter gear system is shown in Figure 17. It is observed that based on HSMM the posterior probability Πk1 firstly exceeds the optimal control limit at file 89, which indicates that the system need to be stopped to take full inspection at this moment.Figure 17
Optimal multivariate Bayesian control chart of the gearbox.


## 4.2. Comparison with Other Maintenance Policies

Using the same maintenance time and cost parameters, the optimal control limitΠ∗=0.28 was obtained in the Bayesian control chart based on HMM, in which the sojourn time in each hidden states follows exponential distribution. The corresponding average cost was equal to 121.72. It should be noted that the average cost concluded by HSMM was less than that by HMM. As shown in Figure 17, while based on HMM, the posterior probability Πk1 firstly exceeds the optimal control limit at file 91. The multivariate Bayesian control chart indicates that the proposed Bayesian control scheme can detect the fault of the gear of interest much earlier than the method based on HMM.Next, we considered the traditional age-based replacement policy, which is well known that it does not take into account condition monitoring information. For the helicopter gearbox, the expected cycle cost (CC) under age-based replacement policy is the sum of failure cost, additional operating and maintenance costs while the system operates in unhealthy state and the preventive maintenance cost:(41)EτCC=CF+CLPTF1-Rτ∣Π⃑0+CAO+CAM∫0τP13t+P14tdt+CP+CLPTPRτ∣Π⃑0,where τ is the stopping time.The expected cycle length (CL) is the sum of the expected time under normal operation, replacement time caused by failure, and preventive maintenance time:(42)EτCL=∫0τRt∣Π⃑0dt+TF1-Rτ∣Π⃑0+TPRτ∣Π⃑0.Then the long-run expected average cost per unit time under age-based replacement policy is equivalent to finding the optimal stopping timeτ∗, such that(43)τ∗=arg⁡inf⁡EτCCEτCL.Table5 shows the optimal values obtained under different policies, which are, respectively, the multivariate Bayesian control chart based on HSMM and HMM and the age-based replacement policy. It is interesting to find that the expected average cost using the Bayesian control chart based on the HSMM is considerably lower than the other policies, while the cost obtained by the age-based replacement policy is the highest due to the irrespective of condition monitoring information. This is reasonable because, for the helicopter gear system, the cost incurred by failure is always higher than the cost incurred by inspection and preventive maintenance measures.Table 5
The optimal maintenance cost under different maintenance policies.
Policy
Optimal parameters
Expected average cost
Π
∗
τ
∗
Bayesian (HSMM)
0.71
—
113.59
Bayesian (HMM)
0.28
—
121.72
Aged-based replacement
—
17.24
140.23


## 5. Conclusions

In this paper, we have proposed a novel optimal Bayesian control scheme for the early fault detection of the partially observable helicopter gear system. The GMR signal was selected for data preprocessing to eliminate the influence of load variation from the original signal of the helicopter gearbox. The healthy portion of GMR signal was used to fit a VAR model; thus the residuals of whole data history were obtained. Without the restriction that the sojourn time in each hidden state was exponentially distributed, the general Erlang distribution was considered for modeling the gear system’s sojourn time in each of the hidden states in a 3-state HSMM. The EM algorithm was employed to estimate the unknown parameters of the model. Using HSMM to describe the deterioration process of the gear system, several important quantities for the gearbox residual life prediction, such as conditional reliability distribution function, PDF and MRL were derived in terms of the posterior probability that the system was in each hidden state. The comparison with the life prediction method based on HMM indicated that the relative errors were smaller using HSMM. An optimal multivariate Bayesian control chart with cost-effective criterion was developed. Moreover, the optimal control limit was solved in the SMDP framework. The comparison with the Bayesian control chart based on HMM as well as the age-based replacement policy indicated that the proposed multivariate Bayesian control scheme provided accurate early fault detection to the helicopter gear system with the minimum cost.In this research, we have considered the multivariate Bayesian control chart under fixed sampling interval; further improvement can be obtained by considering two or more sampling intervals or the availability maximization criterion to obtain the optimal control limit, which can be a suitable topic for the future work.



---

*Source: 4682409-2017-03-23.xml*