Parameter Estimation of a Multistate Model for an Aging Piece of Equipment under Condition-Based Maintenance

Abstract

We study a multistate model for an aging piece of equipment under condition-based maintenance and apply an expectation maximization algorithm to obtain maximum likelihood estimates of the model parameters. Because of the monitoring discontinuity, we cannot observe any state's duration. The observation consists of the equipment's state at an inspection or right after a repair. Based on a proper construction of stochastic processes involved in the model, calculation of some probabilities and expectations becomes tractable. Using these probabilities and expectations, we can apply an expectation maximization algorithm to estimate the parameters in the model. We carry out simulation studies to test the accuracy and the efficiency of the algorithm.

1. Introduction

A multistate model may provide more flexibility than a traditional binary state model for modeling equipment conditions. In a multistate model, a piece of equipment is allowed to experience more than two possible states, for example, completely working, partially working, partially failed, and completely failed. Even if every state has an exponential duration distribution, the equipment has a nonexponentially lifetime distribution. Therefore, many authors utilize multistate models for equipment in many disciplines. Examples arise from electric power systems [1], the nuclear science [2], electron devices [3], and so on. Furthermore, although the deterioration process of a system is usually continuous, it is convenient to model it by several discrete states. For example, in [4], a deterioration process is classified into four states: initial, minor deterioration, major deterioration, and failure. Besides, [5, 6] define several discrete states for a system based on different ranges of the hazard rate of the system; that is, a deterioration process is often modeled by a multistate model.

On the other hand, maintenance is an indispensable action to keep the reliability level of industrial equipment. To reduce the total maintenance costs, condition-based maintenance is usually implemented. References [7, 8] are good overviews of condition-based maintenance. In this paper, condition-based maintenance involves three key ideas: a routine inspection, a condition indicator, and a repair. Routine inspections are executed with a fixed benchmark interval. In an inspection, we can acquire the current condition indicator for the equipment. An appropriate repair is carried out when the indicator is not less than a warning value. Moreover, we also record the condition indicator right after the equipment is repaired. This kind of condition-based maintenance is recommended by an IEEE standard [9], and there are many types of equipment operating in such mode [10, 11].

We assume that every state of the equipment has an exponential duration distribution specified by its parameter. In the paper, we investigate the parameter estimation problem, based on the data obtained from routine inspections and repairs.

Such a parameter estimation problem for a multistate model plays an important role in many applications [12, 13]. The problem attracts much attention, and many approaches are proposed. Reference [5] reviews several quite sophisticated empirical approaches. A method of curve fitting is presented in [12] for the probability density function (pdf) of the equipment lifetime. In [14–16], based on proportional hazards models, different approximate methods are proposed to obtain a maximum likelihood estimate (MLE). In these methods, it is assumed that the condition indicators are time independent or the state of the equipment is fixed within two consecutive inspections. Those assumptions are reasonable only for a small benchmark interval or frequent inspections. However, frequent inspections may increase the total costs and may induce a high chance of accidents.

There is an alternative method to obtain an MLE. To better illustrate the idea, we start from a complete data set of every state’s durations through continuous monitoring. As it has an exponential duration distribution, the reciprocal of the average of durations of a state is the MLE of its parameter. This simple example gives a clue to the case of condition-based maintenance. Then, given the incomplete data from routine inspections and repairs, we can calculate conditional expectations of every state’s duration, and we may deduce an MLE. The theoretical foundation of the above discussion is the expectation maximization (EM) algorithms. The EM algorithms are first introduced in [17] to calculate an MLE for incomplete data. Due to its theoretical property, (see, [18, 19] e.g.,), the EM algorithms have good performance in many applications.

Recently, a few researchers have tried to use EM algorithms for problems similar to ours. For example, based on a complicated discussion, [20] applied an EM algorithm to a specific 3-state model under devised condition-based maintenance. However, there is no hint in [20] to extend the result to a general multistate model. In general, although an EM algorithm may have good performance, it is difficult to apply an EM algorithm to a specific model. Besides the choice of hidden random variables or stochastic processes for a model, the calculation of involved conditional expectations is often difficult.

In this paper, we first propose a proper mathematical framework and a special technique of distributions to simplify the computation of conditional expectations. Then we apply an EM algorithm to obtain MLEs of the model parameters, with the observation of routine inspections and repairs. The paper is organized into seven sections. Section 2 establishes a model for our problem. In Section 3, we apply an EM algorithm, which is based on several conditional expectations, to the model. Section 4 establishes formulas of these conditional expectations. Section 5 presents simulation results, and Section 6 presents an application to a real-world dataset. The paper is concluded in Section 7.

2. Models

In this section, we propose the mathematical framework for our problem. In the following Model 1, we present basic concepts of a new multistate model for an aging piece of equipment under condition-based maintenance. Besides, specific hidden random variables are proposed in Model 1.

Model 1. In the multistate model, a piece of equipment with 𝑚 discrete states under condition-based maintenance can be described by a predefined warning value 𝑊≤𝑚 and two families of random variables 𝑇_𝑖=(𝑇¹_𝑖,…,𝑇_𝑖^𝑚−1) and 𝑅_𝑖=(𝑅¹_𝑖,…,𝑅^𝑚_𝑖), 𝑖=0,1,…. The random variable 𝑇^𝑗_𝑖 has an exponential distribution specified by a parameter 𝜃_𝑗, 𝑗=1,…,𝑚−1.The state space of 𝑅^𝑗_𝑖 is {1,…,𝑚}. Finally, random variables 𝑇^𝑗_𝑖 and 𝑅^𝑗_𝑖 are independent.

In Model 1, for both 𝑇^𝑗_𝑖 and 𝑅^𝑗_𝑖, the subscript 𝑖 indicates a time index, and the superscript 𝑗 presents the state of the equipment at the moment. Concretely, we assume that the state of the equipment is 𝑗 at a time instance 𝑖 and a repair takes place, and 𝑅^𝑗_𝑖 is the state of the equipment right after the repair. And 𝑇^𝑗_𝑖 denotes the residual duration of state 𝑗 after a time instance 𝑖, given that there is no maintenance again. As introduced in Section 1, the duration of state 𝑗 has an exponential distribution, and 𝜃_𝑗 is its parameter. Due to the memoryless property of an exponential distribution, the distribution of the residual duration is also exponential. Here and in the following, we assume that the benchmark interval is 1 for simplicity of argument.

To introduce some derivational concepts to describe the behavior of the equipment clearly, we introduce several families of auxiliary random variables. These random variables, 𝑌_𝑖, 𝑀_𝑖, and 𝐼_𝑖, are recursively defined. First, 𝑀₀=1, 𝐼₀=0.

For 𝑖≥1,𝑌_𝑖=max𝑠∣^𝑠−1_{𝑗=𝑀𝑖−1}𝑇^𝑗_𝑖−1≤1, (2.1) where 𝑌_𝑖 that indicates the state of the equipment at a time instance 𝑖. In fact, assume that the state of the equipment is 𝑀_𝑖−1 at a time instance 𝑖−1, 𝑌_𝑖 is the maximal number 𝑠 such that the total duration of states 𝑀_𝑖−1,…, and 𝑠−1 is less than 1. In other words, between two consecutive inspections at time indexes 𝑖−1 and 𝑖, the equipment goes through states 𝑀_𝑖−1,…,𝑠−1.

For 𝑘≥1,𝐼_𝑘=min𝑖∣𝑖>𝐼_𝑘−1,𝑌_𝑖≥𝑊, (2.2) where 𝐼_𝑘 indicates the time index when the 𝑘th repair takes place.

Moreover,𝑀_𝐼𝑘=𝑅^𝑌𝐼𝑘_𝐼𝑘, (2.3) and for 𝐼_𝑘−1<𝑖<𝐼_𝑘,𝑀_𝑖=𝑌_𝑖. (2.4) Here 𝑀_𝑖 denotes the state of the equipment right after a time instance 𝑖, whether there is maintenance or not. In fact, if 𝑌_𝑖<𝑊 and no repair takes place, we have 𝑀_𝑖=𝑌_𝑖. If 𝑌_𝑖≥𝑊, there is a proper repair, and 𝑀_𝑖 denotes the state of the equipment right after the repair, as showed in (2.3).

3. EM Algorithm

We now apply the EM algorithm to estimate the parameter 𝜃_𝑗, 𝑗=1,…,𝑚−1, given 𝑛 pairs of observations 𝑦_𝑖 and 𝑚_𝑖−1, 𝑖=1,2,…,𝑛.

The EM algorithm, introduced in [17], is an iterative solution to MLE. Each iteration consists of two steps. In the Expectation step, the following functions 𝑄(𝜃^′∣𝜃) are sought:𝑄𝜃^′∣𝜃=𝐸_𝜃log𝐿𝜃^′∣𝑌_𝑖=𝑦_𝑖,𝑀_𝑖−1=𝑚_𝑖−1,𝑖=1,…,𝑛, (3.1) where parameters 𝜃,𝜃^′∈𝑅₊^𝑚−1. 𝜃^′ stands for parameters to be estimated in the current iteration, and 𝜃 denotes the parameters estimated in the previous iteration. As indicated in [17], the expectation in (3.1) is taken over the hidden random variable 𝑇_𝑖 and 𝑅_𝑖 and 𝐿(𝜃^′) is the joint pdf of 𝑇_𝑖 and 𝑅_𝑖, with the parameter 𝜃^′. Estimates of the parameters 𝜃^′ are obtained in the maximization step by solving the following optimization:max_𝜃′𝑄𝜃^′∣𝜃. (3.2) The estimation formulas are established by setting𝜕𝑄𝜃^′∣𝜃𝜕𝜃^′_𝑗=0,𝑗=1,2,…,𝑚−1. (3.3)

As the pdf of the repairs 𝑅_𝑖 has nothing to do with 𝜃^′ and the joint pdf of 𝑇_𝑖, 𝑖=1,2,…,𝑛 has the following form:^𝑛_𝑖=1^𝑚−1_𝑗=1𝜃_𝑗exp−𝜃_𝑗𝑇^𝑗_𝑖, (3.4)(3.3) leads to an unique solution 𝜃^′=ℛ(𝜃). Here ℛ(𝜃)=(ℛ₁(𝜃),…,ℛ_𝑚−1(𝜃)) and for 𝑗=1,…,𝑚−1,ℛ_𝑗(𝜃)=𝑛^𝑛_𝑖=1𝐸_𝜃𝑇^𝑗_𝑖∣𝑌_𝑖=𝑦_𝑖,𝑀_𝑖−1=𝑚_𝑖−1,𝑖=1,…,𝑛⁻¹. (3.5)

Based on the theory of EM algorithms, (see, [17–19] e.g.), we propose the following algorithm to calculate the MLE of 𝜃_𝑗, 𝑗=1,…,𝑚−1.

Algorithm 3.1. Step 1. Let 𝑡=0, and choose an initial value 𝜃⁰∈𝑅₊^𝑚−1.Step 2. Set 𝜃^𝑡+1=ℛ(𝜃^𝑡).Step 3. Stop and output 𝜃^𝑡 if the stop criterion is satisfied.Step 4. Set 𝑡+1→𝑡 and go to Step 2.

4. Conditional Expectations

To calculate the conditional expectations involved in (3.5), we propose some lemmas and theorems in the following.

Lemma 4.1. Assume that random variables 𝑈₁,…,𝑈_𝑟 are independent, 𝑈_𝑖 has an exponential distribution specified by a parameter 𝜆_𝑖, and, for 𝑖≠𝑗, we have 𝜆_𝑖≠𝜆_𝑗. Then, the summation 𝑆=∑^𝑟_𝑖=1𝑈_𝑖 has a pdf 𝑔(𝑥;𝜆₁,…,𝜆_𝑟). Here, 𝑔𝑥;𝜆₁,…,𝜆_𝑟=^𝑟_𝑖=1𝑐_𝑖𝜆_𝑖exp−𝜆_𝑖𝑥 (4.1) for 𝑥≥0. 𝑐_𝑖 is a function for arguments 𝜆₁,…,𝜆_𝑟, defined by 𝑐_𝑖𝜆₁,…,𝜆_𝑟=∏^𝑟_𝑗=1𝜆_𝑗𝜆_𝑖∏^𝑟_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝜆_𝑖. (4.2) The above denominator is 𝜆_𝑖 when 𝑟=1.

The proof of Lemma 4.1 is based on the properties of polynomials and is presented in Appendix A. Further, we can calculate a useful conditional expectation and some basic probabilities based on Lemma 4.1.

Lemma 4.2. Under the assumption and notations in Lemma 4.1, we have the conditional expectation 𝐸(𝑈_𝑖∣𝑆=𝑥)=𝑏_𝑖(𝑥). Here, for 𝑥≥0, 𝑏_𝑖(𝑥)=1𝑔𝑥;𝜆₁,…,𝜆_𝑟^𝑟_{𝑗=1,𝑗≠𝑖}𝑐_𝑗𝜆_𝑗−𝑥𝑒^{−𝜆𝑖𝑥}+𝑒^{−𝜆𝑗𝑥}−𝑒^{−𝜆𝑖𝑥}𝜆_𝑖−𝜆_𝑗. (4.3)

Proof. Let 𝑆^𝑖=𝑆−𝑈_𝑖, then 𝑆^𝑖 and 𝑈_𝑖 are independent. Based on Lemma 4.1, the pdf 𝑓_𝑆𝑖(𝑥) of 𝑆^𝑖 has a form 𝑓_𝑆𝑖(𝑥)=^𝑟_{𝑗=1,𝑗≠𝑖}𝑐_𝑗,𝑖𝜆_𝑗exp−𝜆_𝑗𝑥 (4.4) for 𝑥≥0, where 𝑐_𝑗,𝑖=∏^𝑟_𝑙=1𝜆_𝑙𝜆_𝑖𝜆_𝑗∏^𝑟_{𝑙=1,𝑙≠𝑖,𝑗}𝜆_𝑙−𝜆_𝑗. (4.5)
Let 𝑓(𝑢,𝑥) be the joint pdf of 𝑈_𝑖 and 𝑆. Since 𝑆^𝑖 and 𝑈_𝑖 are independent and 𝑆=𝑈_𝑖+𝑆^𝑖, we have the following:𝑓(𝑢,𝑥)=𝑓_𝑈𝑖(𝑢)⋅𝑓_𝑆𝑖(𝑥−𝑢). (4.6) Hence, the conditional pdf of 𝑈_𝑖 with respect to 𝑆 is 𝑓_{𝑈𝑖|𝑆}()=𝑓𝑢∣𝑥^𝑈𝑖(𝑢)⋅𝑓_𝑆𝑖(𝑥−𝑢)𝑔𝑥;𝜆₁,…,𝜆_𝑟. (4.7) In the sequel, it follows from the expression of 𝑓_𝑆𝑖(𝑥) that 𝐸𝑈_𝑖=1∣𝑆=𝑥𝑔𝑥;𝜆₁,…,𝜆_𝑟^𝑥₀𝑢𝑓_𝑈𝑖(𝑢)𝑓_𝑆𝑖()=𝜆𝑥−𝑢𝑑𝑢_𝑖𝑔𝑥;𝜆₁,…,𝜆_𝑟^𝑟_{𝑗=1,𝑗≠𝑖}𝑐_𝑗,𝑖𝜆_𝑗𝑥𝑒^{−𝜆𝑖𝑥}𝜆_𝑗−𝜆_𝑖+𝑒^{−𝜆𝑗𝑥}−𝑒^{−𝜆𝑖𝑥}𝜆_𝑗−𝜆_𝑖². (4.8) The result follows from the above equation.

Lemma 4.3. Suppose that a random variable 𝑉 has an exponential distribution specified by a parameter 𝜑, random variables 𝑈₁,…,𝑈_𝑟 and 𝑉 are independent, and 𝜑≠𝜆_𝑖 for 1≤𝑖≤𝑟. Under the assumption and notations of Lemma 4.1, we have 𝑃(𝑆≤1)=𝑝(𝜆₁,…,𝜆_𝑟) and 𝑃(𝑆≤1,𝑆+𝑉>1)=𝑝(𝜆₁,…,𝜆_𝑟;𝜑). Here two functions 𝑝(𝜆₁,…,𝜆_𝑟) and 𝑝(𝜆₁,…,𝜆_𝑟,𝜑) have the following form: 𝑝𝜆₁,…,𝜆_𝑟=^𝑟_𝑖=1𝑐_𝑖1−𝑒^−𝜆𝑖,𝑝𝜆₁,…,𝜆_𝑟;𝜑=^𝑟_𝑖=1𝑐_𝑖𝜆_𝑖𝑒^−𝜆𝑖−𝑒^−𝜑𝜑−𝜆_𝑖. (4.9)

Proof. By computing the definite integral of the pdf of 𝑆 on the interval [0,1], we can obtain the expression of 𝑃(𝑆≤1).
It follows from the independence among 𝑈₁,…,𝑈_𝑟 and 𝑉 that 𝑆=∑^𝑟_𝑖=1𝑈_𝑖 and 𝑉 are independent. Therefore, 𝑃(𝑆≤1,𝑆+𝑉>1)=_𝐷𝑔𝑥;𝜆₁,…,𝜆_𝑟𝑓_𝑉(𝑣)𝑑𝑥𝑑𝑣. (4.10) Here and in the following, 𝐷={(𝑥,𝑣)|𝑥≤1,𝑥+𝑣>1}. Besides, 𝑓_𝑉(𝑥) is the pdf of 𝑉. In the sequel, 𝑃(𝑆≤1,𝑆+𝑉>1)=^𝑟_𝑖=1¹₀𝑐_𝑖𝜆_𝑖𝑒^{−𝜆𝑖𝑥∞}_1−𝑥𝜑𝑒^−𝜑𝑣𝑑𝑣𝑑𝑥 (4.11) follows from the above integral.

Then we can calculate several auxiliary important conditional expectations.

Theorem 4.4. Under the assumption and notations of Lemma 4.3, we have 𝐸𝑈_𝑖=𝑔∣𝑆<1^𝑖𝜆₁,…,𝜆_𝑟𝑝𝜆₁,…,𝜆_𝑟,𝐸𝑈_𝑖=𝑔∣𝑆≤1,𝑆+𝑉>1^𝑖𝜆₁,…,𝜆_𝑟;𝜑𝑝𝜆₁,…,𝜆_𝑟,𝐸[]𝜆;𝜑𝑉∣𝑆≤1,𝑆+𝑉>1=ℎ₁,…,𝜆_𝑟;𝜑. (4.12) Here, 𝑔_𝑖(𝜆₁)=(1/𝜆₁)(1−𝑒^−𝜆1)−𝑒^−𝜆1 and, for 𝑟>1, 𝑔_𝑖(𝜆₁,…,𝜆_𝑟) has the following form: ^𝑟_{𝑗=1,𝑗≠𝑖}𝑐_𝑗𝜆_𝑗𝑒^−𝜆𝑖𝜆_𝑖11+𝜆_𝑖−1𝜆_𝑗−𝜆_𝑖+𝑒^−𝜆𝑗𝜆_𝑗𝜆_𝑗−𝜆_𝑖−1𝜆²_𝑖+1𝜆_𝑗𝜆_𝑖. (4.13) Besides, 𝑔_𝑖𝜆₁=𝜆;𝜑¹𝜑−𝜆₁11−𝜑−𝜆₁𝑒^−𝜆1+𝑒^−𝜑𝜑−𝜆₁. (4.14) And, for 𝑟>1, 𝑔_𝑖(𝜆₁,…,𝜆_𝑟;𝜑) has the following form: ^𝑟_{𝑗=1,𝑗≠𝑖}𝑐_𝑗𝜆_𝑗−𝑒^−𝜆𝑖𝜑−𝜆_𝑖−𝑒^−𝜑−𝑒^−𝜆𝑖𝜑−𝜆_𝑖²+1𝜆_𝑖−𝜆_𝑗𝑒^−𝜆𝑗−𝑒^−𝜑𝜑−𝜆_𝑗−𝑒^−𝜆𝑖−𝑒^−𝜑𝜑−𝜆_𝑖. (4.15) Moreover, ℎ(𝜆₁,…,𝜆_𝑟;𝜑) has the following form: 11+𝜑+1𝑝𝜆₁,…,𝜆_𝑟,𝜑^𝑟_𝑗=1𝑐_𝑗𝜆_𝑗𝜑−𝜆_𝑗𝑒^−𝜆𝑗−𝑒^−𝜑𝜑−𝜆_𝑗−𝑒^−𝜆𝑗. (4.16)

Proof. In the following, we calculate three expectations respectively.Case 1. 𝐸[𝑉∣𝑆≤1,𝑆+𝑉>1]. We have the following: 𝐸[]=∬𝑉∣𝑆≤1,𝑆+𝑉>1^𝐷𝑣𝑓_𝑆(𝑥)𝑓_𝑉(𝑣)𝑑𝑥𝑑𝑣𝑃().𝑆≤1,𝑆+𝑉>1 (4.17) Here 𝑓_𝑆(𝑥) and 𝑓_𝑉(𝑣) are pdfs of random variables 𝑆 and 𝑉, respectively. We have 𝑓_𝑆(𝑥)=𝑔(𝑥;𝜆₁,…,𝜆_𝑟). The denominator is presented in Lemma 4.3. Similarly to the proof of Lemma 4.3, the above numerator is _𝐷𝑣𝑓_𝑆(𝑥)𝑓_𝑉(𝑣)𝑑𝑥𝑑𝑣=¹₀^𝑟_𝑖=1𝑐_𝑖𝜆_𝑖𝑒^{−𝜆𝑖𝑥∞}_1−𝑥𝜑𝑣𝑒^−𝜑𝑣=𝑑𝑣𝑑𝑥_𝑟𝑖=1𝑐_𝑖𝜆_𝑖11+𝜑−𝜆_𝑖+1𝜑𝑒^−𝜆𝑖−𝑒^−𝜑𝜑−𝜆_𝑖−^𝑟_𝑖=1𝑐_𝑖𝜆_𝑖𝑒^−𝜆𝑖𝜑−𝜆_𝑖. (4.18) Simplifying the above expression, we can obtain that 𝐸[𝑉∣𝑆]=ℎ(𝜆₁,…,𝜆_𝑟;𝜑).Case 2. We have 𝐸[𝑈_𝑖∣𝑆≤1]. When 𝑟=1, it follows from 𝑆=𝑈₁ that we have 𝐸𝑈₁∣𝑈₁=1≤1𝑃𝑈₁≤1¹₀𝜆₁𝑥𝑒^−𝜆1𝑥1𝑑𝑥=𝜆₁−𝑒^−𝜆11−𝑒^−𝜆1. (4.19) When 𝑟>1, according to the total expectation formula, we have 𝐸𝑈_𝑖∣𝑆≤1=𝐸𝐸𝑈_𝑖∣𝑆∣𝑆≤1. (4.20) Then it follows from Lemma 4.2 and (4.20) that 𝐸𝑈_𝑖∣𝑆≤1⋅𝑝𝜆₁,…,𝜆_𝑟=¹₀𝑏_𝑖(𝑥)𝑓_𝑆(𝑥)𝑑𝑥, (4.21)𝑔_𝑖(𝜆₁,…,𝜆_𝑟) is the above definite integral.Case 3. We have 𝐸[𝑈_𝑖∣𝑆≤1,𝑆+𝑉>1]. When 𝑟=1, it is obvious that 𝐸𝑈₁∣𝑈₁≤1,𝑈₁=∬+𝑉>1^𝐷𝑥𝜆₁𝑒^−𝜆1𝑥𝜑𝑒^−𝜑𝑣𝑑𝑣𝑑𝑥𝑃𝑈₁≤1,𝑈₁+𝑉>1, (4.22) and the result follows.
When 𝑟>1, similarly to the discussion of Case 2, we have𝐸𝑈_𝑖∣𝑆≤1,𝑆+𝑉>1=𝐸𝐸𝑈_𝑖∣𝑆,𝑉∣𝑆≤1,𝑆+𝑉>1. (4.23) As 𝑈₁,…,𝑈_𝑟 and 𝑉 are independent, and 𝑆 is a function of 𝑈₁,…,𝑈_𝑟, we have 𝐸(𝑈_𝑖∣𝑆,𝑉)=𝐸(𝑈_𝑖∣𝑆). Hence, it follows from Lemma 4.2 that 𝐸𝑈_𝑖∣𝑆≤1,𝑆+𝑉>1.𝑝𝜆₁,…,𝜆_𝑟=,𝜑_{10∞1−𝑥}𝜑𝑏_𝑖(𝑥)𝑓_𝑥(𝑥)𝑒^−𝜑𝑣=𝑑𝑣𝑑𝑥_{10𝑟𝑗=1,𝑗≠𝑖}𝑐_𝑗,𝑖𝜆_𝑗𝑒^−𝜑−𝑥𝑒^{(𝜑−𝜆𝑖)𝑥}+𝑒^{(𝜑−𝜆𝑗)𝑥}−𝑒^{(𝜑−𝜆𝑖)𝑥}𝜆_𝑗−𝜆_𝑖𝑑𝑥. (4.24) The result follows from the above equation.

Now we turn to the following Lemma 4.5 based on the theories of measure and stochastic processes. Please see Appendix B for a proof of Lemma 4.5.

Lemma 4.5. For random variables involved in Model 1, we have the following: 𝐸𝑇^𝑗_𝑘∣𝑀_𝑖−1=𝑚_𝑖−1,𝑌_𝑖=𝑦_𝑖,𝑖=1,…,𝑛=𝐸𝑇^𝑗_𝑘∣𝑀_𝑘=𝑚_𝑘,𝑌_𝑘+1=𝑦_𝑘+1. (4.25) Here, 0≤𝑘≤𝑛−1 and 1≤𝑗≤𝑚−1.

Finally, we arrive at conditional expectations in (3.5), which is necessary to Algorithm 3.1.

Theorem 4.6. We assume that there is a guess 𝜃=(𝜃₁,…,𝜃_𝑚−1) for the parameter vector of Model 1 for which 𝜃_𝑖≠𝜃_𝑗 if 𝑖≠𝑗. Let 𝑒𝑥^𝑗_𝑘=𝐸_𝜃[𝑇^𝑗_𝑘∣𝑌_𝑖=𝑦_𝑖,𝑀_𝑖−1=𝑚_𝑘−1,𝑖=1,…,𝑛] for 1≤𝑗≤𝑚−1 and 0≤𝑘≤𝑛−1. When 𝑦_𝑘+1=𝑚_𝑘, we have the following: 𝑒𝑥^𝑗_𝑘=1𝜃_𝑗,𝑗≠𝑦_𝑘+1,11+𝜃_𝑦𝑗,𝑗=𝑦_𝑘+1. (4.26) When 𝑦_𝑘+1=𝑚, we have 𝑒𝑥^𝑗_𝑘=1𝜃_𝑗,𝑗<𝑚_𝑘,𝑔_{𝑗−𝑚𝑘+1}𝜃_𝑚𝑘,…,𝜃_𝑚−1𝑝𝜃_𝑚𝑘,…,𝜃_𝑚−1,𝑚_𝑘≤𝑗≤𝑚−1. (4.27) Finally, when 𝑦_𝑘+1≠𝑚 and 𝑦_𝑘+1≠𝑚_𝑘, we have 𝑒𝑥^𝑗_𝑘=1𝜃_𝑗,𝑗<𝑚_𝑘𝑜𝑟𝑗>𝑦_𝑘+1,𝑔_{𝑗−𝑚𝑘+1}𝜃_𝑚𝑘,…,𝜃_{𝑦𝑘+1−1};𝜃_𝑦𝑘+1𝑝𝜃_𝑚𝑘,…,𝜃_{𝑦𝑘+1−1};𝜃_𝑦𝑘+1,𝑚_𝑘≤𝑗<𝑦_𝑘+1,ℎ𝜃_𝑚𝑘,…,𝜃_{𝑦𝑘+1−1};𝜃_𝑦𝑘+1,𝑗=𝑦_𝑘+1. (4.28)

Proof. Except some trivial cases, this result follows from Theorem 4.4 and Lemma 4.5.

5. Simulation Results

In this section, we test the efficiency and accuracy of the EM algorithm, based on simulation data.

At first, we investigate the necessity of the condition “𝜃_𝑖≠𝜃_𝑗 for 𝑖≠𝑗” in Theorem 4.6. We can find similar conditions for many theorems and lemmas in Section 4, and there is the root in Lemmas 4.1 and 4.2. Although these conditions make our analysis clear and convenient, they appear to be obstacles for an application.

Now we turn to the assumption of Lemma 4.2. Assume that random variables 𝑈₁, 𝑈₂, and 𝑈₃ are independent, and 𝑈_𝑖 has an exponential distribution specified by a parameter 𝜃_𝑖. Based on Lemma 4.2, we can compute the conditional expectation of 𝑈₁ conditional on the summation 𝑆=∑³_𝑖=1𝑈_𝑖 to be any value, for example, 1.

Let 𝜃₃=4 and𝐵𝜃₁,𝜃₂=𝐸𝑈₁∣𝑆=1. (5.1) We draw the surface of 𝐵(𝜃₁,𝜃₂) in Figure 1. Figure 1 shows that although 𝐵(𝜃₁,𝜃₂) does not exist if 𝜃₁=𝜃₂, 𝜃₁=𝜃₃, or 𝜃₂=𝜃₃, but by filling proper values, 𝐵(𝜃₁,𝜃₂) is continuous. That is, 𝐵(𝜃₁,𝜃₂) has some removable discontinuities. The result is obvious from the viewpoint of engineering; why expectations in Section 4 have “jumps" for any special cases? Similar conclusion can be obtained for the assumption of Lemma 4.1 from the viewpoint of mathematics. In fact, the pdf of 𝑈_𝑖 is an analytical function on 𝜃_𝑖>0. As a convolution, the pdf of 𝑆=𝑈₁+𝑈₂+𝑈₃ is also an analytical function on 𝜃_𝑖>0, 𝑖=1,2,3. And any discontinuity is removable for an expression of such a function.

Figure 1 

𝐵(𝜃1,𝜃2) has some removable discontinuities.

Similar discussions can be applied to other theorems and lemmas in Section 4. In a word, these conditions may not trouble us for any application. For a numerical implement of the EM algorithm, if parameters become the same by accident, a small perturbation can fix the problem.

Secondly, we present a method to simulate the behavior of a 4-state model for an aging piece of equipment under condition-based maintenance, based on some preset parameters. The simulation can provide us a test data set.

Without maintenance, the equipment would run through states 1, 2, 3 orderly to a failure state 4. With maintenance, the straight path to a failure is regularly deflected by inspections and maintenance. In an inspection, if we find that the current state is 1 or 2, no repair is applied. If the state is 3, an appropriate repair is carried out. In this case, the state of the equipment right after the repair is a random variable, which is 1 with probability 0.1, is 2 with probability 0.3, and is 3 with probability 0.6. If the state of the equipment is 4, we must replace it with a new piece of equipment, and the state right after this replacement is 1. The duration of states 1, 2, and 3 have exponential distributions specified by parameters 0.3, 0.29, and 0.5. These parameters are suggested by [21].

Now we turn to five experiments in which we test the efficiency and accuracy of the EM algorithm, based on simulation data provided by the above method.

In the first experiment, we study the relationship between convergence error and the number of iterations. As EM algorithms are recursive, this experiment will provide us with an intuitive stop criterion for Algorithm 3.1. For a data set from 100 inspections, Figure 2 shows that all curves of the estimated parameters 𝜃₁, 𝜃₂, and 𝜃₃ converge while the number of iterations increases. We can conclude from Figure 2 that for 100 or more iterations, we can obtain reasonable results.

Figure 2 

Relationship between convergence error and the number of iterations.

In the second experiment, we run the EM algorithm for different initial values. Figure 3 draws the estimated parameter pairs (𝜃₁,𝜃₂) for 100 different initial values which are chosen randomly from 0.001 to 15. We can see that the EM algorithm converges to the same result for a great range of initial values.

Figure 3 

Estimated parameters for different initial values.

In the third experiment, we run the EM algorithm for different sample sizes. Given the number of inspections 𝐼 = 20, 40, 100, 200, 400, and 1000, we run the simulating equipment and obtain 6 data sets. Then based on these data sets, we estimate parameters 𝜃₁, 𝜃₂, and 𝜃₃, respectively. For the EM algorithm, the initial values of estimated parameters are 𝜃₁=𝜃₂=𝜃₃=0.1, and the number of iterations is 100. Table 1 shows the estimated parameters for different data sets. Of course, the estimated parameters are nonmonotonic for different numbers of inspections, which is inherent to such a random problem.

Table 1 

Estimated parameters with different numbers of inspections.

In the fourth experiment, we repeat the third experiment 100 times to obtain the standard errors of each estimated parameter. Table 2 shows the standard errors of the estimated parameters. Again, all cases demonstrate a decrease in the standard errors, when the number of inspections increases. We can conclude from Tables 1 and 2 that, for 100 or more inspections, we can obtain reasonable estimated parameters with an acceptable standard error.

Table 2 

Standard errors of estimated parameters.

Finally, in the fifth experiment, we compare our method with another method.

As indicated in Section 1, [14–16] used an assumption that the state of the equipment is fixed during two consecutive inspections for their problems. Under this assumption, we can establish an MLE for every parameter. For 𝑛 pairs of observations 𝑌_𝑖 and 𝑀_𝑖−1, 𝑖=1,2,…,𝑛, we write 𝑉₁=0, 𝑆₁=𝑀₀, and for 𝑘≥1,𝑉_𝑘=min𝑖≥𝑉_𝑘−1∣𝑀_𝑖≠𝑆_𝑘or𝑀_𝑖≠𝑌_𝑖,𝑆_𝑘=𝑀_𝑉𝑘. (5.2) Here 𝑉_𝑘 records the moment when the state of the equipment is changed and 𝑆_𝑘 records the corresponding new state. Let 𝒩^𝑗_𝑘 be the number of the set {𝑙≤𝑘∣𝑆_𝑙=𝑗}. Then, under the assumption of [14–16],𝑋^𝑗_𝒩𝑗𝑘=𝑉_𝑘+1−𝑉_𝑘 (5.3) records the 𝒩^𝑗_𝑘th duration of the state 𝑗 of the equipment. Finally, let the total number of changing 𝐾 be the maximal subscript of 𝑉_𝑘. Then, we have the following MLE:𝜃_𝑗=𝒩^𝑗_𝐾^𝒩𝑗𝐾_𝑙=1𝑋^𝑗_𝑙⁻¹. (5.4)

On the other hand, although we assume that, in our method, the benchmark interval is 1 for simplicity of argument, it is not too difficult to implement our method for other benchmark intervals. We run the simulating equipment for 400 time units and obtain inspection data sets for different benchmark intervals, respectively.

Figure 4 shows the estimated parameter 𝜃₁ of two methods. We can conclude that for frequent inspections, (the corresponding benchmark interval is less than 0.2) two methods have similar behaviors. However, we can find that two methods show a striking contrast when the benchmark interval is not less than 0.2. As indicated in [21], the latter is more important for practical situations. As our method is stable for all cases, it is more practical.

Figure 4 

Estimated parameters of two methods for different benchmark intervals.

6. Application

In this section, we apply the EM algorithm to a real-world dataset. The dataset is a 25-year record of inspections and repairs of a power transformer substation under condition-based maintenance. The power transformer substation consists of two power transformers, and it is a cool backup system. Once the primary transformer incurs faults, the secondary transformer replaces it automatically, and it can generally serve the entire load. The benchmark interval of inspections is three months. During an inspection, a transformer is repaired if it incurs faults, and the repaired primary transformer is put into operation. The repair is perfect because we can replace a transformer if necessary.

The behavior of the power transformer substation can be described by a 3-state model with a warning value 𝑊=2 and a three-month benchmark interval. Here state 1 implies that the primary transformer works, 2 implies that the secondary transformer works, and 3 implies that two transformers incur faults. In Table 3, we present the record of inspections. Here 𝑌 and 𝑀 are defined by (2.1), (2.3), and(2.4)

Table 3 

Inspections of a power transformer substation.

Based on the dataset given in Table 3, we apply Algorithm 3.1 and obtain the following MLE 𝜃₁ and 𝜃₂ for the primary and the secondary transformers, respectively𝜃₁=0.0243,𝜃₂=0.164. (6.1) That is, the mean duration of the primary transformer is 41.2 months, and it is 6.10 months for the secondary transformer.

We use the total time on test (TTT) plot to check the efficiency of the result. The TTT plot is obtained by plotting𝐺𝑟𝑛=∑^𝑟_𝑖=1𝑍_𝑖∶𝑛+(𝑛−𝑟)𝑍_𝑟∶𝑛∑^𝑛_𝑖=1𝑍_𝑖∶𝑛 (6.2)

against 𝑟/𝑛. Here 𝑍_𝑖∶𝑛 is the order statistics of the sample. See [22] for a detailed introduction to TTT plot. From the dataset of Table t010207, the total lifetime of two transformers are 22, 81, 54, 54, 21, 9, and 45 months. We apply Monte Carlo simulations, with the MLE parameters 𝜃₁ and 𝜃₂, to generate 100 independent total lifetime of two transformers. The TTT plots for the real-world, dataset and the simulated dataset is shown in Figure 5. In Figure 5, the circles represent the real data, and the solid points represent the simulated data with the estimated parameters. Figure 5 shows the efficiency of the result.

Figure 5 

The TTT plots for real-world dataset and the simulated dataset.

7. Conclusion

In this paper, we have addressed the parameter estimation problem for the multistate model for aging equipment under condition-based maintenance. Based on the memoryless property of exponential distributions, we proposed a convenient mathematical framework for the problem. In this framework, the calculation of involved conditional expectations became tractable. Then we applied an EM algorithm to obtain MLE of parameters. A sequence of simulation experiments shows that the estimated parameters converge to the preset value of the parameters, even for a moderate number of inspections. Moreover, the EM algorithm is stable for different length of benchmark intervals. Hence, the algorithm can be recommended for practical applications. It is convenient to extend the algorithm to the situation with random benchmark intervals.

In this paper, we assume that the state of the equipment can be directly observed at an inspection. However, it is a difficult task for many types of industrial equipment. In these cases, some other variables may be used to obtain estimates of the states. To establish a full model consisting of equipment, inspections, indirect variables, and investigating parameter estimations, it is challenging and valuable for practical situations.

Appendices

A. Proof of Lemma 4.1

We proceed by an induction on 𝑟 to show the result in Lemma 4.1.

For 𝑟=1, the result clearly holds. And for the convenience of the following discussion, we analyze the case of 𝑟=2. In this case, the pdf of 𝑆 is the convolution of the pdf 𝑓_𝜆1(𝑥) of 𝑈₁ and the pdf 𝑓_𝜆2(𝑥) of 𝑈₂. Thus, for 𝑥≥0,𝑔𝑥;𝜆₁,𝜆₂=𝑓_𝜆1(𝑥)⊗𝑓_𝜆2(𝑥)=𝜆₂𝜆₂−𝜆₁𝜆₁𝑒^−𝜆1𝑥+𝜆₁𝜆₁−𝜆₂𝜆₂𝑒^−𝜆2𝑥. (A.1)

Assume that for 𝑟=𝑘, (4.1) and (4.2) hold. That is, for 𝑥≥0,𝑔𝑥;𝜆₁,…,𝜆_𝑘=^𝑘_𝑖=1𝑐^𝑘_𝑖𝜆_𝑖exp−𝜆_𝑖𝑥, (A.2) where𝑐^𝑘_𝑖=∏^𝑘_𝑗=1𝜆_𝑖𝜆_𝑖∏^𝑘_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝜆_𝑖. (A.3) Let 𝑟=𝑘+1, we have the following:𝑔𝑥;𝜆₁,…,𝜆_𝑘+1=𝑓_𝜆𝑘+1(𝑥)⊗𝑔𝑥;𝜆₁,…,𝜆_𝑘. (A.4) Due to the linearity of convolution, we can derive from (A.1) that, for 𝑥≥0,𝑔𝑥;𝜆₁,…,𝜆_𝑘+1=^𝑘_𝑖=1𝑐^𝑘_𝑖𝑓_𝜆𝑘+1(𝑥)⊗𝑓_𝜆𝑖(𝑥)=^𝑘+1_𝑖=1𝑐_𝑖^𝑘+1𝜆_𝑖𝑒^{−𝜆𝑖𝑥}. (A.5) Here, for 𝑖≤𝑘,𝑐_𝑖^𝑘+1=𝑐^𝑘_𝑖⋅𝜆_𝑘+1𝜆_𝑘+1−𝜆_𝑖=∏^𝑘+1_𝑗=1𝜆_𝑗𝜆_𝑖∏^𝑘+1_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝜆_𝑖. (A.6) Therefore, (A.1) holds for 𝑖≤𝑘. As to 𝑖=𝑘+1, we have𝑐^𝑘+1_𝑘+1=^𝑘_𝑖=1𝑐^𝑘_𝑖⋅𝜆_𝑖𝜆_𝑖−𝜆_𝑘+1=^𝑘_𝑖=1∏^𝑘_𝑗=1𝜆_𝑗𝜆_𝑖∏^𝑘_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝜆_𝑖⋅𝜆_𝑖𝜆_𝑖−𝜆_𝑘+1=∏^𝑘_𝑗=1𝜆_𝑗∏^𝑘_𝑗=1𝜆_𝑗−𝜆_𝑘+1𝑓𝜆_𝑘+1. (A.7)

It is obvious that𝑓(𝑥)=^𝑘_𝑖=1∏^𝑘_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝑥∏^𝑘_{𝑗=1,𝑗≠𝑖}𝜆_𝑗−𝜆_𝑖 (A.8) is a polynomial in a variable 𝑥, and the degree of 𝑓(𝑥) is less than 𝑘. Moreover, we have 𝑓(𝜆_𝑖)−1=0, 𝑖=1,…,𝑘. Since a polynomial with a degree less than 𝑘 must be a zero polynomial if it has 𝑘 different roots. Hence, 𝑓(𝑥)−1≡0, or 𝑓(𝑥)≡1. Hence,𝑐^𝑘+1_𝑘+1=∏^𝑘_𝑗=1𝜆_𝑗∏^𝑘_𝑗=1𝜆_𝑗−𝜆_𝑘+1. (A.9)

So the result is true for 𝑟=𝑘+1. This completes the proof.

B. Proof of Lemma 4.5

For a time instance 𝑘, let 𝜏₁=max{𝑖|𝐼_𝑖≤𝑘−2} and 𝜏₂=min{𝑖|𝐼_𝑖≥𝑘−1}. We will consider the following filtrations:ℱ₁=𝜎𝑇^𝑗_𝑖,𝑀_𝐼𝑙∣𝑗=1,…,𝑚−1,𝑖≤𝑘−2,𝑙≤𝜏₁,ℱ₂=𝜎𝑇^𝑗_𝑘−1,ℱ∣𝑗=1,…,𝑚−1₃=𝜎𝑇^𝑗_𝑖,𝑀_𝐼𝑙∣𝑗=1,…,𝑚−1,𝑖≥𝑘,𝑙≥𝜏₂. (B.1)

Let 𝐴₁={𝑀_𝑘−1=𝑟} and 𝐵₁={𝑌_𝑘=𝑠} for given 1≤𝑟,𝑠≤𝑚. We will prove that for every events 𝐴∈ℱ₁, 𝐵∈ℱ₂, and 𝐶∈ℱ₃,𝑃𝐴𝐵𝐶∣𝐴₁𝐵₁=𝑃𝐴𝐶∣𝐴₁𝐵₁⋅𝑃𝐵∣𝐴₁𝐵₁. (B.2) That is, ℱ₁∨ℱ₃ and ℱ₂ are independent, conditional on {𝑀_𝑘−1=𝑟} and {𝑌_𝑘=𝑠}.

It is obvious that 𝐴₁∈ℱ₁. Moreover, as𝐵₁=𝑇^𝑘_𝑟>1,𝑠=𝑟,_{𝑚−1𝑗=𝑟}𝑇^𝑗_𝑘−1≤1,𝑠=𝑚,_{𝑠−1𝑗=𝑟}𝑇^𝑗_𝑘−1≤1,^𝑠_𝑗=𝑟𝑇^𝑗_𝑘−1>1,otherwise, (B.3) we have 𝐵₁∈ℱ₂. Hence, due to the independence among ℱ₁, ℱ₂, and ℱ₃, we have the following:𝑃𝐴₁𝐵₁=𝑃𝐴₁𝑃𝐵₁,𝑃𝐴𝐶𝐴₁𝐵₁=𝑃𝐴𝐶𝐴₁𝑃𝐵₁. (B.4) Therefore, we have the following:𝑃𝐴𝐶∣𝐴₁𝐵₁=𝑃𝐴𝐶𝐴₁𝐵₁𝑃𝐴₁𝐵₁=𝑃𝐴𝐶∣𝐴₁,𝑃𝐵∣𝐴₁𝐵₁=𝑃𝐵∣𝐵₁. (B.5) Further, as ℱ₁∨ℱ₃ and ℱ₂ are independent, we obtain that𝑃𝐴𝐵𝐶𝐴₁𝐵₁=𝑃𝐴𝐶𝐴₁𝑃𝐵𝐵₁. (B.6) It follows from (B.5), and the above equation that𝑃𝐴𝐵𝐶∣𝐴₁𝐵₁=𝑃𝐴𝐶∣𝐴₁𝐵₁𝑃𝐵∣𝐴₁𝐵₁. (B.7) This implies (B.2).

Then, we have the following:𝐸𝑇^𝑗_𝑘−1∣𝜎𝑀_𝑘−1,𝑌_𝑘∨ℱ₁∨ℱ₃=𝐸𝑇^𝑗_𝑘−1∣𝑀_𝑘−1,𝑌_𝑘. (B.8) As𝑌_𝑘=𝑦_𝑘,𝑀_𝑘−1=𝑚_𝑘−1∈𝜎𝑀_𝑘,𝑌_𝑘+1∨ℱ₁∨ℱ₃ (B.9) for every 𝑖≤𝑛, taking the {𝑌_𝑖=𝑦_𝑖,𝑀_𝑖−1=𝑚_𝑖−1,𝑖=1,…,𝑛} projection of both sides of (B.8), we have the result𝐸𝑇^𝑗_𝑘−1∣𝑌_𝑖=𝑦_𝑖,𝑀_𝑖−1=𝑚_𝑖−1,1≤𝑖≤𝑛=𝐸𝑇^𝑗_𝑘−1∣𝑌_𝑘=𝑦_𝑘,𝑀_𝑘−1=𝑚_𝑘−1. (B.10)

Acknowledgments

The authors acknowledge the funding of the National Natural Science Foundation of China (Grant Number 50977073). Qihong Duan acknowledges the funding of the National Natural Science Foundation of China (Gant Number 70971109).