A Data-Driven Particle Filter Approach for System-Level Prediction of Remaining Useful Life

PFs are a class of Bayesian filters that use sequential Monte Carlo sampling to approximate the posterior distribution of latent variables with a set of weighted samples or “particles” [7, 14]. While classical filters such as the Kalman Filter are optimal under linear and Gaussian assumptions, PFs offer a more flexible and nonparametric alternative that performs well in nonlinear and non-Gaussian settings. In our work, we adapt the particle filtering framework to sequentially infer the latent degradation parameter $X$ as more observations become available. This allows us to capture both epistemic and stochastic uncertainty in EOL predictions. Below we present the details of this degradation model and the inference procedure.

Bayesian filtering provides a principled framework for sequentially updating our knowledge about a sequence of latent variables as new indirect observations of these variables arrive. At time step $k$, the goal is to estimate the posterior distribution of the latent variable $X_k$ given all observations $y_{0:k}$ up to current time $k$. This posterior is often referred to as the belief distribution and is computed in two recursive steps:

The belief distribution ${\overline{bel}}_k$ thus captures the current probabilistic estimate of $X_k$. Notice that the particle filtering framework depends only on the transition, sensor models:

and the initial distribution of the latent variable $be{l}_0(x)=p(X_0=x)$.

In our setting, the latent process parameters are static, so ideally we have $X_k=X$ for all $k\ge 1$. However, similar to standard particle filtering, we adopt a pseudo-transition model with Gaussian propagation noise to maintain particle diversity and mitigate degeneracy.

Each observation $y_k=(t_k,s_k)$ provides an indirect measurement of the underlying degradation parameters $X$ through the stochastic process $T_s$. This leads to the following sensor and transition models:

where the transition model injects artificial Gaussian noise centered at the previous state. As usual in particle filtering, for simplification, we assume the covariance matrix $\mathrm{\Sigma }$ to be diagonal, which corresponds to injecting independent noise into each state component. Its diagonal entries control the strength of the perturbation, allowing the filter to explore broader regions of the parameter space when uncertainty is high, and to concentrate the posterior distribution as more informative data become available.

The algorithm then follows the standard bootstrap PF framework [7, 8], with a modified correction step tailored for parameter inference. This modification arises from the fact that, unlike in standard particle filtering, not all observations provide the same level of insight into the degradation parameters. Early in the system’s life, observations tend to be less informative, as much of the uncertainty lies in the unobserved future. Consequently, particles that appear unlikely based on current data may still correspond to plausible future scenarios and should not be prematurely discarded. Conversely, near the EOL, observations become more informative and should exert greater influence on particle selection.

To account for this, we introduce a tempering parameter $\kappa >0$, referred to as the observation correction strength, which modulates the influence of the sensor model. We temper the likelihood function by raising it to the power $\kappa$, yielding tempered importance weights:

This approach follows the power posterior formulation in Bayesian inference, where the likelihood is scaled to adjust its impact on the posterior. When , the update is softened to preserve particle diversity under high uncertainty. When $\kappa >1$, the correction becomes sharper and more selective. The standard PF corresponds to the special case $\kappa =1$.

In our framework, both the observation correction strength $\kappa$ and the variances in the diagonal matrix $\mathrm{\Sigma }$ are learned offline from historical data as a function of the observation. Offline data is also used to construct the initial distribution $p(X_0)$. Specifically, we obtain a parameter estimate $x_0^{(i)}$ for each training unit by maximizing the likelihood of its run-to-failure data, for $i=1,2,\mathrm{\dots },n$. These estimates are then augmented with Gaussian noise to generate the desired number of particles $N$, which approximates the initial distribution $p(X_0)$ as:

where ${\delta }_z$ is the Dirac measure centered at $z$. This corresponds to a standard PF initialization, where all particles are equally weighted. The variance of the Gaussian noise used for augmentation is treated as a hyperparameter of the algorithm.

Then the RUL prediction PF we propose can be written as:

1. 1.

   Initialization, $k=0$.

   • $\mathrm{■}$

     For $i=1,\mathrm{\dots },N$, we consider $x_0^{(i)}$ and set $k=1$.

2. 2.

   Prediction step

   • $\mathrm{■}$

     For $i=1,\mathrm{\dots },N$, sample ${\tilde{x}}_k^{(i)}\sim 𝒩(x_{k-1}^{(i)},\mathrm{\Sigma })$

3. 3.

   Correction step

   • $\mathrm{■}$

     For $i=1,\mathrm{\dots },N$, evaluate the importance weights:

     $${\tilde{w}}_k^{(i)}=p{(y_k\mid {\tilde{x}}_k^{(i)})}^{\kappa }.$$

   • $\mathrm{■}$

     Normalise the importance weights.

     $$w_k^{(i)}=\frac{{\tilde{w}}_k^{(i)}}{{\sum }_{j=1}^N{\tilde{w}}_k^{(j)}}.$$

4. 4.

   Selection step

   • $\mathrm{■}$

     Resample with replacement $N$ particles ${\{x_k^{(i)}\}}_{i=1}^N$ from the set ${\{{\tilde{x}}_k^{(i)}\}}_{i=1}^N$ according to the importance weights ${\{w_k^{(i)}\}}_{i=1}^N$.

5. 5.

   Set $k\leftarrow k+1$ and go to step 2.

At each time step $k$, particle filtering maintains an empirical approximation ${\overline{bel}}_k(x)$ of the posterior distribution over the latent degradation parameters. This distribution represents the current belief about the true degradation parameter $X$. In practice, PFs approximate this belief using a discrete posterior distribution composed of Dirac delta functions centered at the particles:

where $w_k^{(i)}$ is the importance weights associated with the particle $x_k^{(i)}$ .

We use this empirical belief distribution to compute the predictive probability distribution of $T_s$, marginalizing over the posterior distribution of $X_k$:

The result is a weighted mixture distribution, where each component corresponds to the performance metric stochastic process where the degradation parameter is the particle, that is, $T_s(x_k^{(i)})$ for each particle $x_k^{(i)}$. This can be interpreted as:

• $\mathrm{■}$

  Each particle $x_k^{(i)}$ represents a distinct hypothesis about the unit’s degradation behavior.

• $\mathrm{■}$

  Its weight $w_k^{(i)}$ reflects the current confidence in that hypothesis.

In particular, the EOL predictive distribution at time step $k$ can be approximated by

From the distribution of the stochastic process ${\{T_s\}}_{s\in [0,1]}$, we compute the EOL prediction at time $k$ as the mean of the EOL predictive distribution:

In general, the distribution $P(T_0=t\mid y_{0:k})$ does not have a closed-form solution, but we can use any numerical method to approximate confidence intervals. To compute a confidence interval at level $\mathrm{\ell }$, we find the quantiles $a$ and $b$ such that

so that .

Finally, the EOL prediction and confidence interval at time step $k$ of the system are computed by applying any aggregation of the corresponding performance predictions. In this work, we use the earliest among the estimated quantities (i.e., mean and quantiles) across all performance metrics as the final prediction and uncertainty bounds.

Figure 1 provides a visual overview of the proposed RUL prediction method at three different time points. Each row corresponds to a distinct stage in the system’s lifecycle.

The left panel in each row displays the predicted RUL of the system, along with a 95% confidence interval computed from the empirical particle distribution. As the system approaches the end of its operational life, the uncertainty decreases and the predictions become more confident. The green dashed line indicates the ground-truth RUL, while the solid blue line and shaded salmon region represent the predicted RUL and its confidence interval, respectively.

The center and right panels depict the predictive probability densities of the stochastic process ${\{T_s\}}_{s\in [0,1]}$, each corresponding to a different performance metric, and are based on their respective current PF estimate of the latent degradation parameters $X$. Yellow, orange, and white points indicate past, current, and future observations, respectively. The horizontal orange bars at the bottom mark the 95% confidence interval of the predicted end-of-life distributions at the current time, and the vertical blue line denotes the predicted RUL. Both are aggregated to produce the RUL prediction and confidence interval shown in the left panel. A brighter point $(t,s)$ in the plot indicates a higher likelihood of reaching performance level $s$ at time $t$. At the second time point, for example, a bimodal predictive distribution emerges, suggesting multiple plausible degradation trajectories consistent with the current particle estimates.

To generate realistic degradation data for evaluating our proposed method, we used a MATLAB implementation of the octo-rotor simulation model introduced in [1]. This high-fidelity framework emulates a fault-tolerant unmanned aerial vehicle (UAV) system by modeling detailed dynamics, control mechanisms, and component degradation under realistic flight conditions. The simulator incorporates empirically validated degradation behaviors for key components such as motor internal resistance, battery capacity, and internal resistance, capturing both component aging and environmental interactions.

Using this platform and the original degradation profiles, we created a custom dataset aligned with our experimental goals. Each UAV is modeled with nine degrading components: eight motors (each with degrading internal resistance) and one battery (modeled with degrading charge capacity and internal resistance). We simulated run-to-failure trajectories for 20 UAVs, each performing repeated missions along a square 1km flight path at a fixed altitude of 50meters. The missions include full take-off and landing cycles and are conducted under a constant 10 m/s wind disturbance in the negative $x$-direction. This simulation environment provides a controlled yet complex test bed for validating prognostic methods under realistic operational and degradation conditions.

Two performance metrics were monitored throughout each UAV’s run-to-failure simulation: a component-level metric, the battery state of charge (SOC), and a system-level metric, the cumulative position error over the flight trajectory. The position error was computed as the accumulated distance between the UAV’s measured position and its projection onto the intended straight-line path between consecutive waypoints. After each flight, these metrics were recorded, and the simulation continued until either metric exceeded a predefined failure threshold. Specifically, the failure threshold for SOC was set to $0.3$, while the threshold for cumulative position error was defined as the equivalent of $1$ meter per second of expected flight time.

Figure 2 illustrates the evolution of both SOC and cumulative position error across the 20 simulated UAVs, with dashed lines indicating the failure thresholds for each metric.

As a post-processing step, the final flight in which failure occurred was removed from each trajectory, and both performance metrics were normalized to the interval $[0,1]$. The maximum healthy value observed across the dataset was mapped to $1$, and the respective failure threshold to $0$. Among the 20 UAVs, failure was triggered by the cumulative position error in 9 cases, and by SOC depletion in 11 cases. We randomly selected 15 UAVs for training and 5 for testing.

The first step of our method is to estimate the initial distribution of the EOL particle filter, denoted $P(X_0)$. To do this, we learn an individual parameter vector $x_0^{(i)}$ for each training unit $i=1,2,\mathrm{\dots },n$ by maximizing the likelihood of its degradation data and using the performance degradation model defined in equation (9). In our dataset, the number of training units is $n=15$. So the initial estimation of the latent distribution $P(X_0)$ is given by equation (5) for the learned values $x_0^{(i)}$, $i=1,2,\mathrm{\dots },15$.

Figure 3 shows the predictive distribution of the particle filter deduced in equation (6),

for the initial step $k=0$ and without considering the augmentation step (i.e. $N=n$). We include the augmentation step in the next subsection because the augmentation noise variance is a hyperparameter.

To learn the correction strength and the variance of the PF propagation noise, we learned a function $f$ that maps the current observation to these six parameters: $\kappa$ (as defined in equation (4)) and the diagonal entries of $\mathrm{\Sigma }$ as defined in equation (3)) corresponding to the noise of the degradation model parameters $(\beta ,\lambda ,\rho ,\mu ,\alpha )$. The function is linear, except for a softplus activation applied to the output to ensure positivity. Since the observation is two-dimensional, $y_k=(t_k,s_k)$, and we included bias terms, we obtained a total of $3\times 6=18$ hyperparameters for learning this function. An additional hyperparameter is used for the variance of the augmentation noise employed to build the initial distribution, resulting in a total of 19 hyperparameters.

We used Optuna [2] for hyperparameter optimization, employing its implementation of the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [4], which is well-suited for efficiently exploring complex and non-convex search spaces. We applied leave-one-out cross-validation. In each fold, the initial distribution was estimated using 14 out of the 15 training units, and the hyperparameters were optimized by evaluating performance on the remaining unit. For the experiments, we set the number of particles at $n_{\text{particles}}=3010$ (we added 10 to make the number of particles a multiple of the 14 training unit fold). We applied our EOL particle filter method independently to both performance metrics: SOC and cumulative position error and at each time step, we evaluated the predictive accuracy using the log-likelihood of future observations under the predictive distribution of the EOL particle filter, as defined in equation (6).

To assess performance, we used three standard prognostic metrics: one for RUL prediction quality (root mean square error, RMSE) and two for uncertainty estimation (prediction interval coverage probability, PICP, and prediction interval normalized width, PINAW). RMSE captures overall prediction accuracy; PICP measures how often the ground-truth RUL falls within the prediction interval; and PINAW captures how tight the interval is, normalized over the range of RUL values.

Figure 4 shows the predicted RUL trajectories along with 95% confidence intervals for both predictive and uncertainty variance across the 5 test UAVs. The experiment was repeated 10 times to assess statistical robustness. The model effectively captures degradation trends and represents both aleatoric and epistemic uncertainty. As the system approaches EOL, the uncertainty naturally decreases due to the accumulation of informative observations. Table 1 reports the numerical results over the full lifetime and the last 50 flights, highlighting improved accuracy in late-stage predictions.