# Remaining Useful Lifetime Prognosis of Controlled Systems: A Case of Stochastically Deteriorating Actuator

Authors: Danh Ngoc Nguyen; Laurence Dieulle; Antoine Grall
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper addresses the case of automatic controlled system which deteriorates during its operation because of components’ wear or deterioration. Depending on its specific closed-loop structure, the controlled system has the ability to compensate for disturbances affecting the actuators which can remain partially hidden. The deterioration modeling and the Remaining Useful Lifetime (RUL) estimation for such closed-loop dynamic system have not been addressed extensively. In this paper, we consider a controlled system with Proportional-Integral-Derivative controller. It is assumed that the actuator is subject to shocks that occur randomly in time. An integrated model is proposed to jointly describe the state of the controlled process and the actuator deterioration. Only the output of the controlled system is available to assess its health condition. By considering a Piecewise Deterministic Markov Process, the RUL of the system can be estimated by a two-step approach. In the first step referred as the “Diagnosis” step, the system state is estimated online from the available monitoring observations by using a particle filtering method. In the second step referred as the “Prognosis” step, the RUL is estimated as a conditional reliability by Monte Carlo simulation. To illustrate the approach, a simulated tank level control system is used.


---

## Body

## 1. Introduction

Due to increasing requirements on durability, reliability, and dependability of industrial systems, intensive research activity on maintenance modeling has been developed during the last decades. Based on the available information about the current system state provided by health monitoring process, different condition-based or predictive maintenance decision rules can be proposed so as to optimize the decision-making process, that is, to prevent or correct failures or faults [1, 2]. In condition-based maintenance framework, a deterioration indicator that correctly describes the dynamic of the failure process is required. Usually this efficient indicator can be constructed from collected information on various deterioration-related monitoring parameters, such as vibration, temperature, lubricating oil, and noise levels. Many research efforts have been devoted to deterioration modeling with increasingly sophisticated approaches which consider different deterioration processes and also dynamic environments [3–5]. However, the need of continuous monitoring in cases of dynamic operating condition may increase the systems costs when expensive monitoring devices are required [6, 7]. In this way, a predictive maintenance policy that schedules maintenance actions according to a prognosis activity without specific additional sensors seems to be an appropriate approach [8, 9].Over the last two decades, numerous prognosis approaches have been developed. According to [6] they can be classified into three main categories: statistical approaches, artificial intelligence approaches, and model-based approaches. The assessment of the Remaining Useful Lifetime (RUL) is one of important tasks in Prognosis. Many studies concentrate on the RUL estimation of systems, subsystems, or components, for example, for lithium-ion batteries [10], rotating machinery [11, 12], or car suspension system [13] (see reviews by [14, 15]).In the field of dependability of automated systems and processes, another research aspect concerns fault-tolerant control (FTC) strategies which give the feedback control system the ability to overcome faults [16]. The key objective of FTC system design may not offer optimal performance in a strict sense for normal operation but can generally mitigate effects of system components failures. For instance, in [17] a scheme for integrated fault detection, diagnosis, and reconfigurable control systems against actuator faults is proposed. The FTC strategies work with the fact that system components can fail by fortuitous causes. In effect, when a fault is detected, the reconfigurable controller and the system command input will be designed automatically to achieve desired performance [18]. However, the dynamic evolution of the deterioration which is the origin of these faults has not been worthy of attention. In fact, because of components’ wear or deterioration the control system performance can be gradually decreasing during operation. In such a case, the information provided by the prognosis process about the health of components and/or the system RUL could be useful to allow modifications of the control actions to continue to achieve the control objectives [19, 20]. In [21] the behavior of the considered system is described by a multiple time scale model which was made up of two parts: a fast dynamic behavior part and a slow dynamic behavior one. The slow dynamic behavior part whose structure is known a priori describes the evolution of the damage state. RUL estimation implies identifying parameters of the structure of the slow dynamic behavior part. Nevertheless, according to the best of our knowledge the stochastic deterioration modeling and RUL estimation process for closed-loop dynamic systems such as feedback control systems have not been addressed extensively.The main aim of this paper is to propose a probabilistic framework to assess the RUL of feedback control systems with stochastically deteriorating actuator and random environment. In the conventional closed-loop control system, the measure of the process output is fed back to the controller in order to generate appropriate control actions on the process. The objective of such systems is to maintain the process output within a desired range defined by a desired set-point. In this paper, the system output measurement is considered as the only available data for health assessment reflecting the deterioration phenomenon. The focus is put on the loss of effectiveness of the actuator. Indeed, the actuators are ones of the most important parts of such systems because they represent the physical link between the control law and the controlled process. In particular, the deterioration of actuators in a closed-loop control system can lead to poor performance and, in extreme cases, loss of controllability. In order to describe the interaction between the deterministic behavior of the feedback control system and the stochastic deterioration process, the whole deteriorating closed-loop system is described as a Piecewise Deterministic Markov Process. In this framework, the distribution of the RUL of the system is computed by using a two-step stochastic model-based technique; see [22].The remainder of this paper is organized as follows. Section2 is devoted to the description of the system characteristics and the assumptions about the condition monitoring process which depends on the stochastic evolution of the set-point. Section 3 describes the approach for computing the Remaining Useful Lifetime which is relevant to system state estimation. To illustrate the methodology, a specific case study is introduced in Section 4. Some numerical results are also discussed here. In this section, the performance of the proposed methodology is compared with a standard cumulative damage model in which the deterioration process is perfectly monitored. Finally, conclusion drawn from this work and possible ways for further studies are given.


## 2. General Modeling Framework

This section describes the characteristics of a deteriorating feedback control system whose actuator stochastically degrades through time and the assumptions about the condition monitoring process which relates to the set-point evolution.


### 2.1. Feedback Control System Structure

In practical control applications, the objective of maintaining the output of a specific process within a desired range is usually achieved by using closed-loop control (see Figure1 for a general scheme of a feedback control system). Sensor gives measurements of process output that are used by the controller in order to calculate the appropriate applied command on the actuator in such a way that reduces the difference between the measured value and the desired set-point to zero or to a small deviation.Figure 1
General block diagram of a feedback control system with notations.We consider a process whose evolution of the states(xt)t∈R+ can be described by the differential equation: (1)x˙t=ft,xt,ut,θ,where x(t)∈Rnx is the vector of the process states, u(t)∈R denotes control variable, θ is the vector of process parameters which is typically considered to be constant in practice, and f(·) is the process dynamic function which can be nonlinear.Despite sophisticated filter structures, noise in the measurement process is usually an unavoidable problem. The measurements of outputs(yt)t∈R+ are then related to the state variables by (2)yt=ht,xt,ut,θ+ϵt,where h(·) is the measurement function which defines the outputs y(t)∈Rny. It is assumed that measurement noises (ϵt)t∈R+ are independent random variables with a probability density, not necessarily Gaussian, independent of the process states (xt)t∈R+.This work focuses on the use of the process output in a prognosis purpose. The possible impact of the control law on the remaining useful life of the system is not investigated in this work. Therefore, the generic usual Proportional-Integral-Derivative (PID) controller structure is used. Due to their simplicity and performance, PID controllers are widely used in industrial applications [23]. Defining uc(t) as the controller output at time t, the standard form of the PID controller is given by (3)uct=KPet+1TI∫0t‍eτdτ+TDdetdt,where e(t) is the error signal defined as e(t)=yref(t)-y(t) with yref(t) as the desired set-point (the reference output), KP is the proportional gain, TI is the integral time, and TD is the derivative time of the PID controller. The adjustment of these three parameters for an optimal system response is extensively studied in control system design [23].The actuating signaluc(t) is then used by the actuator to affect the control action on the process. Actuators are physical devices, for example, control valves, pumps, and other control switches. The output of actuator which is the real control variable is defined as a function g depending on the required value uc(t) of the controller and on the actual capacity of actuator C(t). g is a decreasing function with respect to C(t): (4)ut=guct,Ct.At the initial stage of working, the actuator operates perfectly; that is,C(t)=c0, where c0 is the initial nominal capacity of actuator. However, due to the natural ageing or wear of the mechanical and/or electrical parts of the actuator influenced by nondesired effects of the operating condition, the actuator’s effectiveness C(t) decreases in time and subsequently reduces the control system performance. For instance in a piston pump, wear and corrosion during the operating period lead to gradually enlarging the clearance between valve ball and seat, which will result in decreasing flow rate [24].


### 2.2. Actuator Deterioration

As mentioned above, the actuator deterioration process is considered as a source of performance deterioration in physical system. IfD(t) describes the accumulated deterioration of the actuator up to time t (incapacity unit), the capacity of the actuator at time t before its failure can be expressed as (5)Ct=c0-Dt.As can be seen in the literature, the occurrence of partial loss of effectiveness on an actuator is shown as a discrete phenomenon in time [25, 26]. The loss of effectiveness of the actuator is considered to result from the dynamic evolution of the deterioration process. Such deterioration models known as shock deterioration models have been widely used and the process of shocks’ occurrence times is classically modeled by Poisson processes (see [27, 28]). In this work, the actuator is therefore considered to be subject to a discrete-time deterioration process. It means that the occurrence of deterioration is driven by a mechanism in a specified time interval which leads to an increment of damage as described in [29]. More precisely the isolated points in time corresponding to discrete wear amounts which accumulate gradually are supposed to occur according to a homogeneous Poisson process with intensity λ. The amounts of damage per shock are independently and identically distributed (i.i.d.). Let N(t) denote the total number of shocks up to time t≥0. Then the accumulated deterioration of the actuator at time t is (6)D(t)=∑j=0N(t)‍Wj(N(t)=0,1,2,…),where Wj denotes the damage produced at the jth shock and W0=0. Namely, {D(t),t≥0} is a compound Poisson process.An example of such model can be found in [22] where the leak size of pneumatic valve in the BLEED air system is modeled by a random jump process due to historical maintenance records and experts opinion. In [30] the above model is applied to the life of a storage battery whose capacity decreases with time and with each charge and discharge, until it becomes useless. In [29] a compound Poisson process is considered for leakage current modeling of ultra thin gate oxides in nanotechnology. In this work, the system operating mode which defines the evolution of set-point is supposed to be unchanged. That is the reason that a homogeneous Poisson process is well suited to model the discrete shock instants (see [31] for a more general case where the impact of a random change operating mode is taken into account).Under this modeling assumption, the deterioration impacts the actuator only at discrete times. In case where an actuator has a monotone gradual deterioration behavior, other processes should be considered, for example, the homogeneous Gamma process which can be thought of as the accumulation of an infinite number of small shocks [4].


### 2.3. Set-Point Evolution

According to the demand, for example, of the production process, the desired set-point may change through time. For example, in a chemical process control, the controlled variables (e.g., temperature, flow rate, etc.) have to track the time-varying trajectories that depend on the production phases. As another example let us consider a water supply system for a building. In this case the operation of water pump depends on the use of customers with peaks and off-peak periods during a day. The water demand on the morning and the night can be different to other periods in a day and from one day to the other. Hence, the evolution of set-point should be modeled by a stochastic process that allows to take into account the variability of the production plan on a large time horizon. Moreover, the consideration of a stochastic process for set-point evolution gives the ability to take into account the impact of variable environmental conditions (traditionally modeled as random) affecting the predetermined set-point. Take the cement production process as an example. One step of this process consists in mixing clinker with gypsum and other additives with desired proportions. This is realized with the help of the weight belt conveyors which transport the materials from the storages to the cement mill. In practice, clinker is sensitive to atmospheric conditions (e.g., it can easily clot with the humidity). Therefore in order to ensure the desired quality for cement product, the set-points of the velocities of the gypsum and additives conveyors have to be adapted to the real quantity of clinker on the clinker conveyor which is random [32].Hereafter, the random evolution of the set-point is described by a time-homogeneous Markov chain with a finite state spacerset={r1,r2,…,rm} describing, for example, the m production rates. Let Ytref be the set-point at time t. The evolution of the stochastic process {Ytref,t≥0} is characterized by the transition probability matrix P with the (i,j)th element equal to (7)pij(t)=P(Ys+tref=rj∣Ysref=ri).In Figure2 we can see an illustrative example of evolutive set-point with a peak of demand at T2 and another one at T5. The nominal rate activity can be observed, for example, between T1 and T2 whereas two weak activity periods occur, for example, before T1.Figure 2
An example of set-point evolution.In practical situations, the production/process rate is well known in a near future. This is why we suppose that the change dates to change the production rates are directly observed without errors. In this work, only one set of PID controller parameters is chosen for all known values of set-pointrset.


### 2.4. Condition Monitoring Process

The controlled system output is considered here as the only information available for the actuator health assessment. Moreover, this information is available only at a sequence of inspection times following a change instant of set-point. In fact, such a change provides the opportunity to characterize the dynamics of the controlled system. It is known that the response of system with a change of set-point has two periods: transient and steady-state ones. The former takes place in the short period of time immediately after the change, whereas the latter is usually defined to begin when the output entered and remained within a specified error band (2% or 5% of the change size of set-point). A significant part of the dynamic behavior of the system is shown in the transient period. For this reason, we take only observations of system output in this period. The inspection procedure can be illustrated by Figure 3.Figure 3
Inspection policy.After a change of set-point, a finite number of the system output’s observations could be recorded every time intervalΔtinspec. The observation times corresponding to the ith change of set-point are denoted as Ti1,Ti2,…,TiNi. The total number of observations Ni depends on the time duration between two changes of set-point but Ni≤nmax⁡, ∀i. The health information at these instants is then modeled by the random variables Yi1,Yi2,…,YiNi defined from (2) as (8)Yij=hTij,xTij,uTij,θ+ϵTij.Let us introduce the time of predictionTprog>0 which is set equal to the current time or the date at which the last observation has been recorded. It is the time at which the system health can be estimated given all the collected knowledge and a residual lifetime can be derived. In the sequel and for the sake of the paper clarity the reference to the set-point change will be removed as far as possible. As a consequence if n is the total number of observations until Tprog, the observation dates and corresponding system output will be, respectively, denoted as 0<T1<⋯<Tn=Tprog and Y1,Y2,…,Yn. These sequences are used here as the only available information for the prognosis purpose.


## 2.1. Feedback Control System Structure

In practical control applications, the objective of maintaining the output of a specific process within a desired range is usually achieved by using closed-loop control (see Figure1 for a general scheme of a feedback control system). Sensor gives measurements of process output that are used by the controller in order to calculate the appropriate applied command on the actuator in such a way that reduces the difference between the measured value and the desired set-point to zero or to a small deviation.Figure 1
General block diagram of a feedback control system with notations.We consider a process whose evolution of the states(xt)t∈R+ can be described by the differential equation: (1)x˙t=ft,xt,ut,θ,where x(t)∈Rnx is the vector of the process states, u(t)∈R denotes control variable, θ is the vector of process parameters which is typically considered to be constant in practice, and f(·) is the process dynamic function which can be nonlinear.Despite sophisticated filter structures, noise in the measurement process is usually an unavoidable problem. The measurements of outputs(yt)t∈R+ are then related to the state variables by (2)yt=ht,xt,ut,θ+ϵt,where h(·) is the measurement function which defines the outputs y(t)∈Rny. It is assumed that measurement noises (ϵt)t∈R+ are independent random variables with a probability density, not necessarily Gaussian, independent of the process states (xt)t∈R+.This work focuses on the use of the process output in a prognosis purpose. The possible impact of the control law on the remaining useful life of the system is not investigated in this work. Therefore, the generic usual Proportional-Integral-Derivative (PID) controller structure is used. Due to their simplicity and performance, PID controllers are widely used in industrial applications [23]. Defining uc(t) as the controller output at time t, the standard form of the PID controller is given by (3)uct=KPet+1TI∫0t‍eτdτ+TDdetdt,where e(t) is the error signal defined as e(t)=yref(t)-y(t) with yref(t) as the desired set-point (the reference output), KP is the proportional gain, TI is the integral time, and TD is the derivative time of the PID controller. The adjustment of these three parameters for an optimal system response is extensively studied in control system design [23].The actuating signaluc(t) is then used by the actuator to affect the control action on the process. Actuators are physical devices, for example, control valves, pumps, and other control switches. The output of actuator which is the real control variable is defined as a function g depending on the required value uc(t) of the controller and on the actual capacity of actuator C(t). g is a decreasing function with respect to C(t): (4)ut=guct,Ct.At the initial stage of working, the actuator operates perfectly; that is,C(t)=c0, where c0 is the initial nominal capacity of actuator. However, due to the natural ageing or wear of the mechanical and/or electrical parts of the actuator influenced by nondesired effects of the operating condition, the actuator’s effectiveness C(t) decreases in time and subsequently reduces the control system performance. For instance in a piston pump, wear and corrosion during the operating period lead to gradually enlarging the clearance between valve ball and seat, which will result in decreasing flow rate [24].


## 2.2. Actuator Deterioration

As mentioned above, the actuator deterioration process is considered as a source of performance deterioration in physical system. IfD(t) describes the accumulated deterioration of the actuator up to time t (incapacity unit), the capacity of the actuator at time t before its failure can be expressed as (5)Ct=c0-Dt.As can be seen in the literature, the occurrence of partial loss of effectiveness on an actuator is shown as a discrete phenomenon in time [25, 26]. The loss of effectiveness of the actuator is considered to result from the dynamic evolution of the deterioration process. Such deterioration models known as shock deterioration models have been widely used and the process of shocks’ occurrence times is classically modeled by Poisson processes (see [27, 28]). In this work, the actuator is therefore considered to be subject to a discrete-time deterioration process. It means that the occurrence of deterioration is driven by a mechanism in a specified time interval which leads to an increment of damage as described in [29]. More precisely the isolated points in time corresponding to discrete wear amounts which accumulate gradually are supposed to occur according to a homogeneous Poisson process with intensity λ. The amounts of damage per shock are independently and identically distributed (i.i.d.). Let N(t) denote the total number of shocks up to time t≥0. Then the accumulated deterioration of the actuator at time t is (6)D(t)=∑j=0N(t)‍Wj(N(t)=0,1,2,…),where Wj denotes the damage produced at the jth shock and W0=0. Namely, {D(t),t≥0} is a compound Poisson process.An example of such model can be found in [22] where the leak size of pneumatic valve in the BLEED air system is modeled by a random jump process due to historical maintenance records and experts opinion. In [30] the above model is applied to the life of a storage battery whose capacity decreases with time and with each charge and discharge, until it becomes useless. In [29] a compound Poisson process is considered for leakage current modeling of ultra thin gate oxides in nanotechnology. In this work, the system operating mode which defines the evolution of set-point is supposed to be unchanged. That is the reason that a homogeneous Poisson process is well suited to model the discrete shock instants (see [31] for a more general case where the impact of a random change operating mode is taken into account).Under this modeling assumption, the deterioration impacts the actuator only at discrete times. In case where an actuator has a monotone gradual deterioration behavior, other processes should be considered, for example, the homogeneous Gamma process which can be thought of as the accumulation of an infinite number of small shocks [4].


## 2.3. Set-Point Evolution

According to the demand, for example, of the production process, the desired set-point may change through time. For example, in a chemical process control, the controlled variables (e.g., temperature, flow rate, etc.) have to track the time-varying trajectories that depend on the production phases. As another example let us consider a water supply system for a building. In this case the operation of water pump depends on the use of customers with peaks and off-peak periods during a day. The water demand on the morning and the night can be different to other periods in a day and from one day to the other. Hence, the evolution of set-point should be modeled by a stochastic process that allows to take into account the variability of the production plan on a large time horizon. Moreover, the consideration of a stochastic process for set-point evolution gives the ability to take into account the impact of variable environmental conditions (traditionally modeled as random) affecting the predetermined set-point. Take the cement production process as an example. One step of this process consists in mixing clinker with gypsum and other additives with desired proportions. This is realized with the help of the weight belt conveyors which transport the materials from the storages to the cement mill. In practice, clinker is sensitive to atmospheric conditions (e.g., it can easily clot with the humidity). Therefore in order to ensure the desired quality for cement product, the set-points of the velocities of the gypsum and additives conveyors have to be adapted to the real quantity of clinker on the clinker conveyor which is random [32].Hereafter, the random evolution of the set-point is described by a time-homogeneous Markov chain with a finite state spacerset={r1,r2,…,rm} describing, for example, the m production rates. Let Ytref be the set-point at time t. The evolution of the stochastic process {Ytref,t≥0} is characterized by the transition probability matrix P with the (i,j)th element equal to (7)pij(t)=P(Ys+tref=rj∣Ysref=ri).In Figure2 we can see an illustrative example of evolutive set-point with a peak of demand at T2 and another one at T5. The nominal rate activity can be observed, for example, between T1 and T2 whereas two weak activity periods occur, for example, before T1.Figure 2
An example of set-point evolution.In practical situations, the production/process rate is well known in a near future. This is why we suppose that the change dates to change the production rates are directly observed without errors. In this work, only one set of PID controller parameters is chosen for all known values of set-pointrset.


## 2.4. Condition Monitoring Process

The controlled system output is considered here as the only information available for the actuator health assessment. Moreover, this information is available only at a sequence of inspection times following a change instant of set-point. In fact, such a change provides the opportunity to characterize the dynamics of the controlled system. It is known that the response of system with a change of set-point has two periods: transient and steady-state ones. The former takes place in the short period of time immediately after the change, whereas the latter is usually defined to begin when the output entered and remained within a specified error band (2% or 5% of the change size of set-point). A significant part of the dynamic behavior of the system is shown in the transient period. For this reason, we take only observations of system output in this period. The inspection procedure can be illustrated by Figure 3.Figure 3
Inspection policy.After a change of set-point, a finite number of the system output’s observations could be recorded every time intervalΔtinspec. The observation times corresponding to the ith change of set-point are denoted as Ti1,Ti2,…,TiNi. The total number of observations Ni depends on the time duration between two changes of set-point but Ni≤nmax⁡, ∀i. The health information at these instants is then modeled by the random variables Yi1,Yi2,…,YiNi defined from (2) as (8)Yij=hTij,xTij,uTij,θ+ϵTij.Let us introduce the time of predictionTprog>0 which is set equal to the current time or the date at which the last observation has been recorded. It is the time at which the system health can be estimated given all the collected knowledge and a residual lifetime can be derived. In the sequel and for the sake of the paper clarity the reference to the set-point change will be removed as far as possible. As a consequence if n is the total number of observations until Tprog, the observation dates and corresponding system output will be, respectively, denoted as 0<T1<⋯<Tn=Tprog and Y1,Y2,…,Yn. These sequences are used here as the only available information for the prognosis purpose.


## 3. Remaining Useful Lifetime Assessment Methodology

Due to interaction between stochastic deterioration and deterministic dynamics of the system, the whole integrated closed-loop deteriorating system can be modeled by a Piecewise Deterministic Markov Process (PDMP). Such a class of models has been first introduced by [33]. PDMP is used to model fatigue growth in [34] and corrosion in [35]. Between two successive shocks reducing the actuator capacity, the response of the closed-loop system is described by differential equations which combine the process dynamic characteristics and PID controller behaviors. The randomness impacts the system and intersects its trajectories only at the random discrete times of shock according to the actuator deterioration model description. Interested readers can refer to [22, 33, 36] for a rigorous mathematical definition of a PDMP.Roughly speaking, the characterization of a PDMP requires three basic elements: a probability lawdFz which represents the law of the time Ti+1-Ti before the next jump given the position Zi=z, a Markovian kernel Q(z;du) which represents the probability law driving process position after a jump from position z, and a flow ψ describing the deterministic trajectories between the jumps.


### 3.1. Two-Step Technique for RUL Assessment

According to the previous description of system and its deterioration process the whole behavior of the deteriorating closed-loop system at timet can be resumed by the extended random vector Zt such that(9)Zt=xtCt,where xt is the state variables of the controlled process and Ct is the actual capacity variable related to the deterioration of the actuator. Note that in the paper modeling framework, the vector Zt defines a time-homogeneous Markov process given that the set-point is a Markovian process. More generally, for example, when the set-point depends deterministically on time, the time t can be included explicitly as a component of Zt for the process to be homogeneous in time [22]. Indeed, if a process Z is Markovian but non-time-homogeneous, then Z~=(Zt,t) is a time-homogeneous Markov process [37]. The considered process (Zt)t∈R+ is a PDMP which is perfectly defined by (i)
the probability lawdFz which is related to the intensity λ of the Poisson process; that is, dFz(v)=λexp⁡(-λv)dv,(ii)
the Markovian kernelQ(z;du) which is the density of the probability law that has been chosen to describe the amounts of damage Wi and which does not depend on z in this paper without any loss of generality,(iii)
the functionψ(z,t) which characterizes the solutions of the ordinary differential equation (1).The kernel of the Markov renewal process is given byN(z,du,dv)=dFz(v)Q(z;du).In the context of the feedback control system, the system failure zoneF is gathering all the unacceptable deterioration states of the actuator. While in a “failed” state the system can still work, but it is unable to fulfill its requirements anymore. The objectives of the control system are not achieved. Practically, the actual capacity of the actuator has to be greater than a minimal capacity level related to the objectives of control system design. The Remaining Useful Lifetime at time tRULt is thus defined as the hitting time of the failure zone F in the process state space; that is (10)RULt=inf⁡s≥t,Zs∈F-t.As stated in [22] the modeling framework of PDMP ensures that the distribution of the system RUL at time Tprog given the online monitoring information up to time Tprog can be written as (11)P(RULTprog>s∣Y1=y1,…,Yn=yn)=∫‍Rzsμy1,…,yndz,where (i)
μy1,…,yn(dz) is the probability law of the system state at time Tprog regarding the available observations y1,…,yn: (12)μy1,…,yn=L(ZTprog∣Y1=y1,…,Yn=yn);(ii)
Rz(s) is the reliability of the system at time s knowing that the initial state value is z: (13)Rz(s)=P(Zu∉F∀u≤s∣Z0=z).The probability density function (pdf) or the mean value ofRULTprog can be derived from (11). The two-step technique detailed in the next paragraphs consists in firstly estimating μy1,…,yn and secondly the conditional reliability knowing ZTprog. Figure 4 illustrates the proposed methodology.Figure 4
Illustration of the prognosis concept.


### 3.2. Step  1: Particle Filtering State Estimation

As mentioned previously the system evolution is modeled using a PDMPZ=(Zt)t∈R+. The sequence, ZT0:k={ZTi,i=0,…,k}, k≤n, where ZT0 is the initial state of the system, is not observed directly but solely and partially through observations Y1:k={Yi,i=1,…,k} as described in Section 2.4. Hence the first objective is to make inferences on the states ZT0:k from the measured values y1:k=y1,…,yk of the observation process Y1:k. More specifically, the main task is to estimate the conditional density, p(zTk∣y1:k), which represents the probability law of the state at time Tk given the observations available up to inspection time Tk. A particle filtering method is considered which allows recursive updates of the density as new observations arrive. First the recursive Bayesian filter is defined by (14)p(zTk∣y1:k-1)=∫‍p(zTk∣zTk-1)p(zTk-1∣y1:k-1)dzTk-1,p(zTk∣y1:k)=p(yk∣zTk)p(zTk∣y1:k-1)p(yk∣y1:k-1),where the quantity p(yk∣y1:k-1) is given by (15)p(yk∣y1:k-1)=∫‍p(yk∣zTk)p(zTk∣y1:k-1)dzTk.The difficulty to implement the recursive Bayesian filter is that the integrals calculations are intractable. Thus, particle filtering is used here to allow for numerical computation of the filtering density. It is a sequential Monte Carlo method particularly useful for optimal estimation and prediction problems in nonlinear non-Gaussian processes [38]. The key idea is to approximate the targeted filtering distribution p(zTk∣y1:k) by a cloud of Ns i.i.d. random samples called particles {zTk(i),i=1,…,Ns} with associated weights {wTk(i),i=1,…,Ns}, which satisfy ∑i=1Ns‍wTk(i)=1. The target distribution at time Tk can be approximated by (16)p(zTk∣y1:k)≈p^(zTk∣y1:k)=∑i=1Ns‍wTk(i)δzTk(i)(dzTk),where δzTk(i)(·) is the Dirac delta mass located in zTk(i).The used particle filter is similar to the Generic Particle Filter in [38] with deterministic resampling method because it seems to be a computationally cheaper algorithm [39]. Indeed, resampling is used to avoid the problem of degeneracy of the algorithm, that is, avoiding the situation that all but one of the importance weights are close to zero [40]. The algorithm uses the prior distribution p(zTk∣zTk-1(i)) derived from (1) to (7) as the importance function.Therefore, the real-time state estimation procedure, given the sequence of measurementy1:k, can be resumed by the algorithm in Algorithm 1.Algorithm 1:Generic particle filter for system state estimation.
Initialization: ∀i=1,…,Ns.
Draw particlezT0(i) according to the initial condition of system
Assign corresponding weightwT0(i)=1/Ns
At step
k
(corresponding to time  Tk): Given zTk-1(i),wTk-1(i)i=1Ns, do
(a) Importance sampling
Based on the system description (derived from (1) to (7)), draw particles z~Tk(i)~p(zTk∣zTk-1(i))
(b) Weight update
Based on the likelihoods of the observationsyk collected (Eq. (8)), assign weights  wTk(i)=wTk-1(i)p(yk∣z~Tk(i))
(c) Weight normalisation
w
T
k
(
i
)
=
w
T
k
(
i
)
∑
i
=
1
N
s
‍
w
T
k
(
i
)
(d) Re-sampling decision
IfN^eff=1/∑i=1Ns‍(wTk(i))2<Nthresh then perform deterministic re-sampling: z~Tk(i),wTk(i)i=1Ns⇒zTki,1/Nsi=1Ns
(e) Distribution
p
(
z
T
k
∣
y
1
:
k
)
≈
∑
i
=
1
N
s
‍
w
T
k
(
i
)
δ
z
T
k
(
i
)
(
d
z
T
k
)
Repeat till the prognosis instant  Tprog  is reached


### 3.3. Step  2: RUL Estimation

The second step of the methodology considered in this paper for the RUL computation requires the estimation of the system reliability starting from the prognosis instantTprog and knowing the approximated pdf of the system state at Tprog as given by (16).Actually, the reliability is computed with the classical Monte Carlo method. It means that the simulation of trajectories of system until its failure is required. The departure point of each trajectory is then randomly selected from the obtained particle set at timeTprog. Each particle is propagated forward to the failure zone by using the future evolution of set-point which is described in Section 2.3. The histogram of the RUL is obtained straightforwardly. The mean value or quantiles of the RUL can also be derived. The procedure is illustrated by Algorithm 2.Algorithm 2:RUL estimation.
GivenzTn(i),wTn(i)i=1Ns, Ndepart number of departure points, Ntraj number of simulation trajectories for each point
For  j=1,…,Ndepart  do
(i) Generate uniform sample:uj~U(0,1)
(ii) Select depart point:
zjselected=zTn(k)with∑l=1k-1‍wTn(l)≤uj<∑l=1k‍wTn(l)
(iii)For  k=1,…,Ntraj  do
Simulate the trajectories according to the system description (derived from (1) to (7))
End
End
Obtain the empirical distribution of RUL


## 3.1. Two-Step Technique for RUL Assessment

According to the previous description of system and its deterioration process the whole behavior of the deteriorating closed-loop system at timet can be resumed by the extended random vector Zt such that(9)Zt=xtCt,where xt is the state variables of the controlled process and Ct is the actual capacity variable related to the deterioration of the actuator. Note that in the paper modeling framework, the vector Zt defines a time-homogeneous Markov process given that the set-point is a Markovian process. More generally, for example, when the set-point depends deterministically on time, the time t can be included explicitly as a component of Zt for the process to be homogeneous in time [22]. Indeed, if a process Z is Markovian but non-time-homogeneous, then Z~=(Zt,t) is a time-homogeneous Markov process [37]. The considered process (Zt)t∈R+ is a PDMP which is perfectly defined by (i)
the probability lawdFz which is related to the intensity λ of the Poisson process; that is, dFz(v)=λexp⁡(-λv)dv,(ii)
the Markovian kernelQ(z;du) which is the density of the probability law that has been chosen to describe the amounts of damage Wi and which does not depend on z in this paper without any loss of generality,(iii)
the functionψ(z,t) which characterizes the solutions of the ordinary differential equation (1).The kernel of the Markov renewal process is given byN(z,du,dv)=dFz(v)Q(z;du).In the context of the feedback control system, the system failure zoneF is gathering all the unacceptable deterioration states of the actuator. While in a “failed” state the system can still work, but it is unable to fulfill its requirements anymore. The objectives of the control system are not achieved. Practically, the actual capacity of the actuator has to be greater than a minimal capacity level related to the objectives of control system design. The Remaining Useful Lifetime at time tRULt is thus defined as the hitting time of the failure zone F in the process state space; that is (10)RULt=inf⁡s≥t,Zs∈F-t.As stated in [22] the modeling framework of PDMP ensures that the distribution of the system RUL at time Tprog given the online monitoring information up to time Tprog can be written as (11)P(RULTprog>s∣Y1=y1,…,Yn=yn)=∫‍Rzsμy1,…,yndz,where (i)
μy1,…,yn(dz) is the probability law of the system state at time Tprog regarding the available observations y1,…,yn: (12)μy1,…,yn=L(ZTprog∣Y1=y1,…,Yn=yn);(ii)
Rz(s) is the reliability of the system at time s knowing that the initial state value is z: (13)Rz(s)=P(Zu∉F∀u≤s∣Z0=z).The probability density function (pdf) or the mean value ofRULTprog can be derived from (11). The two-step technique detailed in the next paragraphs consists in firstly estimating μy1,…,yn and secondly the conditional reliability knowing ZTprog. Figure 4 illustrates the proposed methodology.Figure 4
Illustration of the prognosis concept.


## 3.2. Step  1: Particle Filtering State Estimation

As mentioned previously the system evolution is modeled using a PDMPZ=(Zt)t∈R+. The sequence, ZT0:k={ZTi,i=0,…,k}, k≤n, where ZT0 is the initial state of the system, is not observed directly but solely and partially through observations Y1:k={Yi,i=1,…,k} as described in Section 2.4. Hence the first objective is to make inferences on the states ZT0:k from the measured values y1:k=y1,…,yk of the observation process Y1:k. More specifically, the main task is to estimate the conditional density, p(zTk∣y1:k), which represents the probability law of the state at time Tk given the observations available up to inspection time Tk. A particle filtering method is considered which allows recursive updates of the density as new observations arrive. First the recursive Bayesian filter is defined by (14)p(zTk∣y1:k-1)=∫‍p(zTk∣zTk-1)p(zTk-1∣y1:k-1)dzTk-1,p(zTk∣y1:k)=p(yk∣zTk)p(zTk∣y1:k-1)p(yk∣y1:k-1),where the quantity p(yk∣y1:k-1) is given by (15)p(yk∣y1:k-1)=∫‍p(yk∣zTk)p(zTk∣y1:k-1)dzTk.The difficulty to implement the recursive Bayesian filter is that the integrals calculations are intractable. Thus, particle filtering is used here to allow for numerical computation of the filtering density. It is a sequential Monte Carlo method particularly useful for optimal estimation and prediction problems in nonlinear non-Gaussian processes [38]. The key idea is to approximate the targeted filtering distribution p(zTk∣y1:k) by a cloud of Ns i.i.d. random samples called particles {zTk(i),i=1,…,Ns} with associated weights {wTk(i),i=1,…,Ns}, which satisfy ∑i=1Ns‍wTk(i)=1. The target distribution at time Tk can be approximated by (16)p(zTk∣y1:k)≈p^(zTk∣y1:k)=∑i=1Ns‍wTk(i)δzTk(i)(dzTk),where δzTk(i)(·) is the Dirac delta mass located in zTk(i).The used particle filter is similar to the Generic Particle Filter in [38] with deterministic resampling method because it seems to be a computationally cheaper algorithm [39]. Indeed, resampling is used to avoid the problem of degeneracy of the algorithm, that is, avoiding the situation that all but one of the importance weights are close to zero [40]. The algorithm uses the prior distribution p(zTk∣zTk-1(i)) derived from (1) to (7) as the importance function.Therefore, the real-time state estimation procedure, given the sequence of measurementy1:k, can be resumed by the algorithm in Algorithm 1.Algorithm 1:Generic particle filter for system state estimation.
Initialization: ∀i=1,…,Ns.
Draw particlezT0(i) according to the initial condition of system
Assign corresponding weightwT0(i)=1/Ns
At step
k
(corresponding to time  Tk): Given zTk-1(i),wTk-1(i)i=1Ns, do
(a) Importance sampling
Based on the system description (derived from (1) to (7)), draw particles z~Tk(i)~p(zTk∣zTk-1(i))
(b) Weight update
Based on the likelihoods of the observationsyk collected (Eq. (8)), assign weights  wTk(i)=wTk-1(i)p(yk∣z~Tk(i))
(c) Weight normalisation
w
T
k
(
i
)
=
w
T
k
(
i
)
∑
i
=
1
N
s
‍
w
T
k
(
i
)
(d) Re-sampling decision
IfN^eff=1/∑i=1Ns‍(wTk(i))2<Nthresh then perform deterministic re-sampling: z~Tk(i),wTk(i)i=1Ns⇒zTki,1/Nsi=1Ns
(e) Distribution
p
(
z
T
k
∣
y
1
:
k
)
≈
∑
i
=
1
N
s
‍
w
T
k
(
i
)
δ
z
T
k
(
i
)
(
d
z
T
k
)
Repeat till the prognosis instant  Tprog  is reached


## 3.3. Step  2: RUL Estimation

The second step of the methodology considered in this paper for the RUL computation requires the estimation of the system reliability starting from the prognosis instantTprog and knowing the approximated pdf of the system state at Tprog as given by (16).Actually, the reliability is computed with the classical Monte Carlo method. It means that the simulation of trajectories of system until its failure is required. The departure point of each trajectory is then randomly selected from the obtained particle set at timeTprog. Each particle is propagated forward to the failure zone by using the future evolution of set-point which is described in Section 2.3. The histogram of the RUL is obtained straightforwardly. The mean value or quantiles of the RUL can also be derived. The procedure is illustrated by Algorithm 2.Algorithm 2:RUL estimation.
GivenzTn(i),wTn(i)i=1Ns, Ndepart number of departure points, Ntraj number of simulation trajectories for each point
For  j=1,…,Ndepart  do
(i) Generate uniform sample:uj~U(0,1)
(ii) Select depart point:
zjselected=zTn(k)with∑l=1k-1‍wTn(l)≤uj<∑l=1k‍wTn(l)
(iii)For  k=1,…,Ntraj  do
Simulate the trajectories according to the system description (derived from (1) to (7))
End
End
Obtain the empirical distribution of RUL


## 4. Case Study: A Double-Tank Level Control System

In the previous section, a methodology to compute the conditional pdf of the RUL of a closed-loop dynamic system was described. Here, it is illustrated on a well-known feedback control system: a double-tank level control system.


### 4.1. Description of the Case Study

Consider a double-tank level system with cross-sectional area of the first tankS1 and the second one S2. Water or other incompressible fluid (i.e., the mass density of fluid ρ is constant) is pumped into the first tank at the top by pump motor drives. Then, the outflow from the first tank feeds the second tank. The water level of tank 2 is measured by a level measurement sensor and controlled by adjusting the pump motor control input which is calculated by a PID controller. The overall tank level control system is shown in Figure 5.Figure 5
A double-tank level control system.In order to consider the real response of pump motor, the relation between the inlet flow rateqin and the pump motor control input uc is represented as a first order system [41]:(17)dqindt=-1τaqin+Kaτauc,where τa is the time constant of pump motor and Ka is the servo amplifier gain (with the initial gain Kainit). The pump saturates at a maximum input umax⁡ and it cannot draw water from the tank, so uc∈[0,umax⁡].The fluid leaves out at the bottom of each tank through valves with the flow rates according to the Torricelli rule:(18)qj,out=Kvj2ghj,j=1,2,where hj is level of tank j, g is the acceleration of gravity, and Kvj is the specified parameter of the valve j.Using the mass balance equation, the process can be described by the following equations:(19)dh1(t)dt=1S1qin-Kv1S12gh1t,dh2(t)dt=Kv1S12gh1(t)-Kv2S22gh2t.The control objective of the system is to adjust the level of tank 2 according to the set-point evolution. To have simple and comprehensible case study, we suppose that the set-point admits only two valuesr1 and r2 with r1<r2. The sojourn times in the different values of system set-point are characterized by a continuous-time Markov chain whose transition rate matrix is(20)P=-α1α1α2-α2,where α1 and α2 describe transition rate of set-point and the mean sojourn time value ri is equal to 1/αi.Due to deterioration of the pump, its capacityKa stochastically decreases. Each time ξi a shock occurs according to a Poisson process with intensity λ, the capacity of pump C(t)=Ka(t)=Kainit-D(t) decreases by a quantity Wi which follows a uniform distribution on [0;Δ].Under all these considerations, the behavior of water tank level control system can be summed up using the processZ=(Zt)t∈R+, where Zt is given by(21)Zt=Ka(t)h1(t)h2(t).The current state of the system at timet is then a three-component vector Zt, which includes the water levels of both tanks and the current capacity of the pump. Note that only the water level h2(t) of tank 2 is observed.According to (17) and (19), the steady states are obtained at instant tss if (22)uctss=S1S2Kv2Katss2gh2tss.Sinceuc(tss)≤umax⁡, then (23)Ka(tss)≥S1S2Kv2umax⁡2gh2(tss).This condition shows that the required actuator capacity depends on the evolution of set-point process. Indeed, if the set-point takes a small value (i.e., r1), not much controlled actions are needed. However, in order to keep a desired level of safety for the system, especially in the case of a random time-varying set-point, and also to simplify the definition of the zone of failure, the failure threshold is defined as the minimal capacity of actuator which gives the system’s ability to handle all possible values of set-point. The minimal capacity can be defined in the control system design phase. In this case of study, this accepted value is defined as (24)Kamin⁡=S1S2Kv2umax⁡2gmax⁡iri=S1S2Kv2umax⁡2gr2.Thus, the RUL of the system is the remaining time before the processZ enters in the failure zone which is defined as (25)Kat≤Kamin⁡.As mentioned above, the system response (the water level of the tank 2) is considered as the only available health information of the system. Indeed, the water level measurement is recorded for prognosis purposes if a change of the set-point is detected. Only the information of the transient response is used because it is more informative than the steady-state response.


### 4.2. Numerical Illustrations

In order to numerically implement the double-tank level control system, the continuous process model (19) and the actuator model (17) are discretized through the forward Euler scheme (with time step Δt). For PID implementation, the velocity algorithm [42] is used. Therefore, the behavior of control system can be resumed by (26)tk=tk-1+Δt,yk=h2k+ϵk,ykref~P(rk∣rk-1),ek=ykref-yk,uck=uck-1+KP1+ΔtTI+TDΔtek+-1-2TDΔtek-1+TDΔtek-2,uck=umax⁡:uck>umax⁡0:uck<0uck:0≤uck≤umax⁡,Dk=∑j=0N(tk)‍Wj,Kak=Kainit-Dk,qk+1=qk1-Δtτak+KakΔtτakuck,h1k+1=h1k+Δt1S1qk-Kv1S12gh1k,h2k+1=h2k+ΔtKv1S12gh1k-Kv2S22gh2k,where yk is the measurement of the water level at time tk given by a level measurement sensor; the measurement noise ϵk is supposed to be an independent Gaussian random variable with standard deviation σ and mean equal to zero: ϵk~N(0,σ2), N(tk) is the total number of shocks up to time tk≥0. Please note that the white noise is a classical frame for noise modeling and this choice does not affect the performance of the proposed RUL estimation methodology which can handle non-Gaussian noise (see Section 3.2).Numerical values for double-tank level control system are summed up in Table1.Table 1
Double-tank model.
Physical parameters
S
1
=
25
K
v
1
=
8
τ
a
=
1
S
2
=
20
K
v
2
=
6
g
=
9.82
u
max
⁡
=
100
σ
=
0.05
PID controller parameters
K
P
=
4.2519
T
I
=
18.9817
T
D
=
1.6182
Initial condition:t=0
h
1
(
0
)
=
0
h
2
(
0
)
=
0
K
a
init
=
5.0
Natural deterioration
λ
=
1
0
-
3
Δ
=
0.5
Varying set-point
r
1
=
10
α
1
=
0.003
r
2
=
25
α
2
=
0.004Figure6 represents one simulated trajectory of the process Z until the failure of system. The evolution of set-point with successive change of set-point values is illustrated in Figure 6(a). The water levels of tank 1 and tank 2, h1(t) and h2(t), are reflected in Figures 6(b) and 6(c). Figure 6(d) shows a simulated trajectory of the actuator capacity. It is depicted for illustration purpose as if a sensor had been added. Practically the actuator capacity is unavailable (i.e., not measured) on the considered system for diagnosis and prognosis. This is a health indicator that we aim to estimate at some discrete times. The inlet flow rate and the control value applied on the actuator are illustrated by Figure 7.Figure 6
A trajectory of the water tank level control system until failure of actuator: (a) set-point, (b) water level of tank 1, (c) water level of tank 2, and (d) actuator capacity.
(a)
(b)
(c)
(d)Figure 7
Control value applied on the actuator (a) and corresponding inlet flow rate (b).
(a)
(b)As depicted in Figure6, the actuator fails completely (i.e., Ka=0) at 26411.6 time units, but the failure of system here is 18951.8 time units. In effect, one can find that after the system failure instant the water level of tank 2 (the controlled variable) cannot track the evolution of desired set-point.Let us now consider an inspection procedure described in Section2.4. Here, we suppose that the change of set-point is immediately detected and nmax⁡=6 the maximum possible observations are recorded with each time duration Δtinspec=4. The methodology previously described is applied to deduce prognosis about the RUL of system. The first step of the method is to compute the pdf of the system state regarding the available observations until the prognosis time Tprog, for example, at Tprog=15046.8 time units, that is, 238th inspection date. The available health information of the system (the noisy observations of the water level of tank 2) at the inspection times before Tprog is shown in Figure 8.Figure 8
Noisy observations of water level of tank 2.The described particle filtering method is applied in order to estimate the conditional state of the system knowing the noisy measurement ofh2. Approximations of the pdfs are represented in Figure 9(a) for the water level of tank 1, Figure 9(b) for the water level of tank 2, and Figure 9(c) for the actuator capacity with Ns=1000 particles.Figure 9
Conditional distribution of the system state at timeTprog=15046.8 time units given the noisy measurements of h2 for Ns=1000 particles.
(a)
Conditional water level of tank 1
(b)
Conditional water level of tank 2
(c)
Conditional actuator capacity parameterThe last step of the method is to compute the distribution of the RUL of the system starting atTprog knowing the approximated pdf of the system state at Tprog. The RUL distribution has been obtained by Monte Carlo simulation with 1500 trajectories describing the system evolution from its state at the prognosis time until its failure. The resulting RUL is depicted in Figure 10. The point estimate of the RUL can be calculated using the RUL distribution. One can find here that the mean value of RUL is very close to the real value.Figure 10
Remaining Useful Lifetime of the water tank level control system at timeTprog=15046.8 time units.The estimation of the system conditional reliability required in the second step of methodology can also be obtained using two other Monte Carlo estimators based on the decomposition of the PDMP process [22]. Another point of view is to consider the distribution of the RUL as a total probability: (27)Pr⁡(T-t>s∣Y1=y1,…,Yn=yn)=∫‍Pr⁡(T-t>s∣Ka(t)=x)fKat∣y1,…,yn(x)(dx).In (27), Pr⁡(T-t>s∣Ka(t)=x) is the conditional survival function of the sytem considering that the actuator capacity random variable is equal to x at time t. Depending on the deterioration model the conditional survival function may be explicitly obtained. fKat∣y1,…,yn(x) is the estimated pdf of the actuator capacity at time t considering Y1=y1,…,Yn=yn. It can be obtained by the first step of the proposed methodology. This integral in (27) can be numerically computed with classical schemes when an analytic expression of Pr⁡(T-t>s∣Ka(t)) is available.As explained previously, the calculation of the RUL using the numerical integration is not always available because of the complexity ofPr⁡(T-t>s∣Ka(t)=x). In our case study, it can be explicitly deduced with the assumption of a compound Poisson shock process for the actuator’s deterioration process and uniform distribution for damage quantity per shock (cf. (28)). Figure 11 compares the survival function of the RUL at time Tprog=15046.8 time units which is estimated by two methods: the proposed 2-step technique and the numerical integral. The result shows the performance of the methodology.Figure 11
Comparison of the Remaining Useful Lifetime estimation of the water tank level control system at timeTprog=15046.8 time units by two methods: Monte Carlo simulation and numerical computation using the estimated pdf of system state.In cases of complex deterioration models of actuator [31], the presented methodology of RUL estimation in the paper which is based on an intelligent Monte Carlo simulation shows its performance.


### 4.3. Impact of Monitoring Information

The presented method proposes a technique in order to compute the RUL of the double-tank level control system based on the partial and imperfect monitoring information. Only the noisy observations of system output that partially reflect the deterioration are available. In order to assess the accuracy of the proposed method, let us consider the case of a direct and perfect monitoring of the deterioration process. In this case, the deterioration level of pump can be perfectly observed which means thatKa(t) is observable and perfectly known at inspection time. The deterioration process of the pump is then modeled by a compound Poisson process (CPP) [43] or a cumulative damage model [28] as illustrated in Figure 12.Figure 12
Cumulative damage model of actuator.Under the supposition that the damage per shock follows a uniform distribution on[0;Δ], the survival function of the RUL in this ideal case knowing that Ka(t)=x with x>Kamin⁡ is given by (see Appendix for the details of calculation) (28)Pr⁡⁡(T-t>s∣Kat=x)=∑n=0∞‍λsnn!e-λs1Δnn!·∑k=0⌊(x-Kamin⁡)/Δ⌋‍-1knkx-Kamin⁡-kΔn,where the notion u is the floor function of u.The Conditional Mean Remaining Useful Lifetime (MRUL) of the pump at aget is then described as(29)ET-t∣Kat=x=1λ∑k=0⌊(x-Kamin⁡)/Δ⌋‍-1k·1k!e((x-Kamin⁡)/Δ-k)x-Kamin⁡Δ-kk.Figure13 depicts the mean time-to-failure calculated by (29) (in case of perfect monitoring information) and the mean time-to-failure estimated by the proposed method (in partial and noisy information case) starting at some different prognosis instants. One can notice the performance of our proposed methodology. For a better lecture of Figure 13 the 95%, 75%, 25%, and 5% quantiles are pointed out.Figure 13
Mean time-to-failure (MTTF) in two cases of monitoring information at some different prognosis instants.


## 4.1. Description of the Case Study

Consider a double-tank level system with cross-sectional area of the first tankS1 and the second one S2. Water or other incompressible fluid (i.e., the mass density of fluid ρ is constant) is pumped into the first tank at the top by pump motor drives. Then, the outflow from the first tank feeds the second tank. The water level of tank 2 is measured by a level measurement sensor and controlled by adjusting the pump motor control input which is calculated by a PID controller. The overall tank level control system is shown in Figure 5.Figure 5
A double-tank level control system.In order to consider the real response of pump motor, the relation between the inlet flow rateqin and the pump motor control input uc is represented as a first order system [41]:(17)dqindt=-1τaqin+Kaτauc,where τa is the time constant of pump motor and Ka is the servo amplifier gain (with the initial gain Kainit). The pump saturates at a maximum input umax⁡ and it cannot draw water from the tank, so uc∈[0,umax⁡].The fluid leaves out at the bottom of each tank through valves with the flow rates according to the Torricelli rule:(18)qj,out=Kvj2ghj,j=1,2,where hj is level of tank j, g is the acceleration of gravity, and Kvj is the specified parameter of the valve j.Using the mass balance equation, the process can be described by the following equations:(19)dh1(t)dt=1S1qin-Kv1S12gh1t,dh2(t)dt=Kv1S12gh1(t)-Kv2S22gh2t.The control objective of the system is to adjust the level of tank 2 according to the set-point evolution. To have simple and comprehensible case study, we suppose that the set-point admits only two valuesr1 and r2 with r1<r2. The sojourn times in the different values of system set-point are characterized by a continuous-time Markov chain whose transition rate matrix is(20)P=-α1α1α2-α2,where α1 and α2 describe transition rate of set-point and the mean sojourn time value ri is equal to 1/αi.Due to deterioration of the pump, its capacityKa stochastically decreases. Each time ξi a shock occurs according to a Poisson process with intensity λ, the capacity of pump C(t)=Ka(t)=Kainit-D(t) decreases by a quantity Wi which follows a uniform distribution on [0;Δ].Under all these considerations, the behavior of water tank level control system can be summed up using the processZ=(Zt)t∈R+, where Zt is given by(21)Zt=Ka(t)h1(t)h2(t).The current state of the system at timet is then a three-component vector Zt, which includes the water levels of both tanks and the current capacity of the pump. Note that only the water level h2(t) of tank 2 is observed.According to (17) and (19), the steady states are obtained at instant tss if (22)uctss=S1S2Kv2Katss2gh2tss.Sinceuc(tss)≤umax⁡, then (23)Ka(tss)≥S1S2Kv2umax⁡2gh2(tss).This condition shows that the required actuator capacity depends on the evolution of set-point process. Indeed, if the set-point takes a small value (i.e., r1), not much controlled actions are needed. However, in order to keep a desired level of safety for the system, especially in the case of a random time-varying set-point, and also to simplify the definition of the zone of failure, the failure threshold is defined as the minimal capacity of actuator which gives the system’s ability to handle all possible values of set-point. The minimal capacity can be defined in the control system design phase. In this case of study, this accepted value is defined as (24)Kamin⁡=S1S2Kv2umax⁡2gmax⁡iri=S1S2Kv2umax⁡2gr2.Thus, the RUL of the system is the remaining time before the processZ enters in the failure zone which is defined as (25)Kat≤Kamin⁡.As mentioned above, the system response (the water level of the tank 2) is considered as the only available health information of the system. Indeed, the water level measurement is recorded for prognosis purposes if a change of the set-point is detected. Only the information of the transient response is used because it is more informative than the steady-state response.


## 4.2. Numerical Illustrations

In order to numerically implement the double-tank level control system, the continuous process model (19) and the actuator model (17) are discretized through the forward Euler scheme (with time step Δt). For PID implementation, the velocity algorithm [42] is used. Therefore, the behavior of control system can be resumed by (26)tk=tk-1+Δt,yk=h2k+ϵk,ykref~P(rk∣rk-1),ek=ykref-yk,uck=uck-1+KP1+ΔtTI+TDΔtek+-1-2TDΔtek-1+TDΔtek-2,uck=umax⁡:uck>umax⁡0:uck<0uck:0≤uck≤umax⁡,Dk=∑j=0N(tk)‍Wj,Kak=Kainit-Dk,qk+1=qk1-Δtτak+KakΔtτakuck,h1k+1=h1k+Δt1S1qk-Kv1S12gh1k,h2k+1=h2k+ΔtKv1S12gh1k-Kv2S22gh2k,where yk is the measurement of the water level at time tk given by a level measurement sensor; the measurement noise ϵk is supposed to be an independent Gaussian random variable with standard deviation σ and mean equal to zero: ϵk~N(0,σ2), N(tk) is the total number of shocks up to time tk≥0. Please note that the white noise is a classical frame for noise modeling and this choice does not affect the performance of the proposed RUL estimation methodology which can handle non-Gaussian noise (see Section 3.2).Numerical values for double-tank level control system are summed up in Table1.Table 1
Double-tank model.
Physical parameters
S
1
=
25
K
v
1
=
8
τ
a
=
1
S
2
=
20
K
v
2
=
6
g
=
9.82
u
max
⁡
=
100
σ
=
0.05
PID controller parameters
K
P
=
4.2519
T
I
=
18.9817
T
D
=
1.6182
Initial condition:t=0
h
1
(
0
)
=
0
h
2
(
0
)
=
0
K
a
init
=
5.0
Natural deterioration
λ
=
1
0
-
3
Δ
=
0.5
Varying set-point
r
1
=
10
α
1
=
0.003
r
2
=
25
α
2
=
0.004Figure6 represents one simulated trajectory of the process Z until the failure of system. The evolution of set-point with successive change of set-point values is illustrated in Figure 6(a). The water levels of tank 1 and tank 2, h1(t) and h2(t), are reflected in Figures 6(b) and 6(c). Figure 6(d) shows a simulated trajectory of the actuator capacity. It is depicted for illustration purpose as if a sensor had been added. Practically the actuator capacity is unavailable (i.e., not measured) on the considered system for diagnosis and prognosis. This is a health indicator that we aim to estimate at some discrete times. The inlet flow rate and the control value applied on the actuator are illustrated by Figure 7.Figure 6
A trajectory of the water tank level control system until failure of actuator: (a) set-point, (b) water level of tank 1, (c) water level of tank 2, and (d) actuator capacity.
(a)
(b)
(c)
(d)Figure 7
Control value applied on the actuator (a) and corresponding inlet flow rate (b).
(a)
(b)As depicted in Figure6, the actuator fails completely (i.e., Ka=0) at 26411.6 time units, but the failure of system here is 18951.8 time units. In effect, one can find that after the system failure instant the water level of tank 2 (the controlled variable) cannot track the evolution of desired set-point.Let us now consider an inspection procedure described in Section2.4. Here, we suppose that the change of set-point is immediately detected and nmax⁡=6 the maximum possible observations are recorded with each time duration Δtinspec=4. The methodology previously described is applied to deduce prognosis about the RUL of system. The first step of the method is to compute the pdf of the system state regarding the available observations until the prognosis time Tprog, for example, at Tprog=15046.8 time units, that is, 238th inspection date. The available health information of the system (the noisy observations of the water level of tank 2) at the inspection times before Tprog is shown in Figure 8.Figure 8
Noisy observations of water level of tank 2.The described particle filtering method is applied in order to estimate the conditional state of the system knowing the noisy measurement ofh2. Approximations of the pdfs are represented in Figure 9(a) for the water level of tank 1, Figure 9(b) for the water level of tank 2, and Figure 9(c) for the actuator capacity with Ns=1000 particles.Figure 9
Conditional distribution of the system state at timeTprog=15046.8 time units given the noisy measurements of h2 for Ns=1000 particles.
(a)
Conditional water level of tank 1
(b)
Conditional water level of tank 2
(c)
Conditional actuator capacity parameterThe last step of the method is to compute the distribution of the RUL of the system starting atTprog knowing the approximated pdf of the system state at Tprog. The RUL distribution has been obtained by Monte Carlo simulation with 1500 trajectories describing the system evolution from its state at the prognosis time until its failure. The resulting RUL is depicted in Figure 10. The point estimate of the RUL can be calculated using the RUL distribution. One can find here that the mean value of RUL is very close to the real value.Figure 10
Remaining Useful Lifetime of the water tank level control system at timeTprog=15046.8 time units.The estimation of the system conditional reliability required in the second step of methodology can also be obtained using two other Monte Carlo estimators based on the decomposition of the PDMP process [22]. Another point of view is to consider the distribution of the RUL as a total probability: (27)Pr⁡(T-t>s∣Y1=y1,…,Yn=yn)=∫‍Pr⁡(T-t>s∣Ka(t)=x)fKat∣y1,…,yn(x)(dx).In (27), Pr⁡(T-t>s∣Ka(t)=x) is the conditional survival function of the sytem considering that the actuator capacity random variable is equal to x at time t. Depending on the deterioration model the conditional survival function may be explicitly obtained. fKat∣y1,…,yn(x) is the estimated pdf of the actuator capacity at time t considering Y1=y1,…,Yn=yn. It can be obtained by the first step of the proposed methodology. This integral in (27) can be numerically computed with classical schemes when an analytic expression of Pr⁡(T-t>s∣Ka(t)) is available.As explained previously, the calculation of the RUL using the numerical integration is not always available because of the complexity ofPr⁡(T-t>s∣Ka(t)=x). In our case study, it can be explicitly deduced with the assumption of a compound Poisson shock process for the actuator’s deterioration process and uniform distribution for damage quantity per shock (cf. (28)). Figure 11 compares the survival function of the RUL at time Tprog=15046.8 time units which is estimated by two methods: the proposed 2-step technique and the numerical integral. The result shows the performance of the methodology.Figure 11
Comparison of the Remaining Useful Lifetime estimation of the water tank level control system at timeTprog=15046.8 time units by two methods: Monte Carlo simulation and numerical computation using the estimated pdf of system state.In cases of complex deterioration models of actuator [31], the presented methodology of RUL estimation in the paper which is based on an intelligent Monte Carlo simulation shows its performance.


## 4.3. Impact of Monitoring Information

The presented method proposes a technique in order to compute the RUL of the double-tank level control system based on the partial and imperfect monitoring information. Only the noisy observations of system output that partially reflect the deterioration are available. In order to assess the accuracy of the proposed method, let us consider the case of a direct and perfect monitoring of the deterioration process. In this case, the deterioration level of pump can be perfectly observed which means thatKa(t) is observable and perfectly known at inspection time. The deterioration process of the pump is then modeled by a compound Poisson process (CPP) [43] or a cumulative damage model [28] as illustrated in Figure 12.Figure 12
Cumulative damage model of actuator.Under the supposition that the damage per shock follows a uniform distribution on[0;Δ], the survival function of the RUL in this ideal case knowing that Ka(t)=x with x>Kamin⁡ is given by (see Appendix for the details of calculation) (28)Pr⁡⁡(T-t>s∣Kat=x)=∑n=0∞‍λsnn!e-λs1Δnn!·∑k=0⌊(x-Kamin⁡)/Δ⌋‍-1knkx-Kamin⁡-kΔn,where the notion u is the floor function of u.The Conditional Mean Remaining Useful Lifetime (MRUL) of the pump at aget is then described as(29)ET-t∣Kat=x=1λ∑k=0⌊(x-Kamin⁡)/Δ⌋‍-1k·1k!e((x-Kamin⁡)/Δ-k)x-Kamin⁡Δ-kk.Figure13 depicts the mean time-to-failure calculated by (29) (in case of perfect monitoring information) and the mean time-to-failure estimated by the proposed method (in partial and noisy information case) starting at some different prognosis instants. One can notice the performance of our proposed methodology. For a better lecture of Figure 13 the 95%, 75%, 25%, and 5% quantiles are pointed out.Figure 13
Mean time-to-failure (MTTF) in two cases of monitoring information at some different prognosis instants.


## 5. Conclusion

This work is a proposal for a positioning of the problem of the RUL evaluation of a dynamic control system with a stochastically deteriorating actuator and aims to combine the dynamic and the stochastic part system modeling using only the output of the system. The present paper proposes a modeling framework using PDMP that shows the ability to combine the deterministic behavior of a feedback control system with the stochastic deterioration process for the actuator. In this framework, the loss of effectiveness of actuator is modeled by the random gaps which intersect the deterministic trajectory of closed-loop system only at random discrete times. Particle filtering technique is used to estimate online the state of considered system regarding only the noisy observations of closed system output. By using a methodology based on the assumption of Markov property, the Remaining Useful Lifetime can be deduced with Monte Carlo simulation. A simulated double-tank level control system was used as a case study to illustrate the efficiency of the proposed approach.Future research will be focused on the use of the estimation of the system state and the RUL in order to optimize the decision-making process. The decision related to inspections and preventive/corrective actions should be considered using the RUL information for the purpose of the cost reduction. Another perspective relates to deterioration modeling of actuator. On one hand, the actuator is less efficient through time because of natural deterioration process. On the other hand, the set-point level impacts also the deterioration process of actuator. For example, in a centrifugal pump, an increased demand of pump flow will cause bearing friction and impeller wear to increase at a faster rate. Hence, the impact of the mission profile (the evolution of set-point) should be addressed.



---

*Source: 356916-2015-06-08.xml*