Predicting the Remaining Useful Life of Solid Oxide Fuel Cell Systems Using Adaptive Trend Models of Health Indicators

1. Introduction

SOFC systems represent an attractive technology for the conversion of electrochemical energy into electrical power due to lower emissions, high efficiency, useful waste heat and fuel flexibility. They are being widely applied in auxiliary power unit (APU) and stationary power generators [1]. Yet, the most significant barriers to their broader commercialization and market deployment are insufficient durability, reliability and cost. Reliability and durability, in particular, are significantly impacted by the degradation phenomena that emerge during the operation at high temperatures.

Understanding the mechanisms that cause the deterioration is of paramount importance for the design, operation and management of the SOFC systems. Much of the research endeavour has been focused on modelling the multiphysics and multiscale processes from micro-structural to system level to understand the conditions that cause the degradation and its rate as well as the relationship with the reliability and the life span [2,3]. Due to the complexity of the underlying models and the required computational load, their use in online diagnosis and prognosis is not feasible.

Thanks to the increasing volume of historical data acquired from the rising number of SOFC systems under test in the laboratory and in-situ operation, room for the use of different machine learning approaches to either predict the degradation trend or estimate the RUL has become widely open. Particularly strong research and development of methods of prognosis have been conducted in the related areas of battery research and, to a lesser extent, in the domain of fuel cells. A comprehensive survey of prognostic methods for batteries can be found in [4]. For example in [5] over 20,000 electrochemical impedance spectroscopy (EIS) spectra of commercial Li-ion batteries, collected at different states of health, states of charge and temperatures were applied to a machine learning algorithm to automatically determine which spectral features predict degradation and the remaining life.

Although less intensive compared to the batteries, the amount of research in data-driven prognostics in PEM is still rather solid and a variety of machine learning techniques have been proposed. A good overview can be found in [6].

So far considerably less attention has been paid to the development of algorithms for online prediction of RUL in SOFC systems so the list of references is rather constrained. A brief address can be found in the review paper [7]. In an early attempt in [8] the authors adopt a neural network model in which time is used as an explicit independent variable to describe the decay of voltage as a function of time and operational parameters. The authors in [9] study the RUL prediction under two degradation modes, i.e. anode poisoning and cathode humidification. From a series of simulated life-long tests they learn the parameters of the transition model by the data-driven method (least squares – support vector machine). However, a change in the health metric, voltage in this case, can be directly affected by operating conditions imposed on the system. Therefore, it is difficult to distinguish whether it is caused by the degradation of the fuel cells, a change in load condition or if it is a response of the system due to a controller action.

Deriving the data-driven model of the stack voltage as a function of the process variables has been pursued in [10] by using a dynamic neural network model. The idea is to train the model in the nominal condition and then take the discrepancy between measured and reconstructed stack voltage as a health indicator and anticipate its trend.

In the paper [11] the health indicator is associated with the area specific resistance (ASR), which is largely affected by the stack degradation and can be made mostly invariant to the process variables. The idea is based on a simple premise that the degradation can be viewed as a hidden random walk process. To deconvolute the degradation dynamics from the available process data the authors developed a nonlinear lumped model of the process and the black-box model of the degradation dynamics. By using filtering techniques, the temporal evolution of the ASR is obtained. Based on that the RUL predictions are evaluated using Monte Carlo (MC) simulation. An attempt to include more prior knowledge in the degradation model, in [12] the authors showed that under proper conditions the quality of prognostics can be slightly improved compared to the pure black box model.

The hidden random walk models of the degradation might face two issues. First, it should be noted that a range of phenomenological models have been proposed in [13] based on a thorough investigation of the degradation phenomena under accelerated tests. The decision of which model to apply depends on the correct identification of the degradation mode. Incorrect diagnosis implicates the use of the incorrect degradation model, which might, in turn, lead to misleading results. Second, it can happen in practice that one degradation mode eventually unfolds other degradation modes so that such cascades are immensely difficult to identify.

Moreover, from past experimental runs, it seems the degradation progresses mainly as a monotone process [14] or even piecewise linear [15]. Once the health indicator is selected, the prediction of its future evaluation can be anticipated by using time series analysis tools. In quest of finding the optimum trade-off between the simplicity of the prognostic algorithm and the quality of the prognosis, in this paper, we apply a probabilistic linear trend analysis. The idea is to estimate the local trend using ordinary least squares (OLS) and then estimate the RUL as the time when the health indicator is anticipated to hit a limit value that marks the beginning of the failure state. A novelty of this approach is that the distribution of the first-hitting time (FHT) can be expressed in closed form meaning that the computational load is almost nil. Thus the algorithm can be easily applied to industrial controllers with limited computational power.

The paper is structured as follows. In section 2 the method of estimating a linear trend from the data is presented. Then the main idea that leads to the closed form solution of the probability distribution of the FHT from the linear trend model is presented. In section 3, the performance of the method is demonstrated on the different simulated time series scenarios. Finally, in section 4, the method is demonstrated on experimental data obtained from a long-term run of a SOFC system.

2. Methodology

2.1. Local Linear Trend Model

The observed degradation metric is taken as a time series. A rectangular window function of length M is used to select a part of the data - usually the latest M measurements. For each instance of windowed data, we propose the latent state of the health index, denoted as x, to behave as a linear function

$$\begin{array}{c}\hfill x\left(t\right)=n+kt,n,k\in \mathbb{R},\end{array}$$

(1)

while the observed state, denoted as y, is just the latent state with superimposed normally distributed white noise

$$\begin{array}{c}\hfill y\left(t\right)=x\left(t\right)+\eta \left(t\right),\eta \left(t\right)\sim \mathcal{N}(0,{\sigma }_{\eta }^2),{\sigma }_{\eta }\in {\mathbb{R}}^{+}.\end{array}$$

(2)

Since we only get our measurements at discrete times $t_j,j\in (1,\cdots ,N)$, we can rewrite the upper equation to matrix form

$$\begin{array}{c}\hfill \mathit{Y}=\left[\begin{array}{c}y\left(t_{N-M+1}\right)\\ y\left(t_{N-M}\right)\\ ⋮\\ y\left(t_N\right)\end{array}\right]=\mathit{F}\left[\begin{array}{c}n\\ k\end{array}\right]+\left[\begin{array}{c}\eta \left(t_{N-M+1}\right)\\ \eta \left(t_{N-M}\right)\\ ⋮\\ \eta \left(t_N\right)\end{array}\right],\mathrm{where}\mathit{F}=\left[\begin{array}{c}1,t_{N-M+1}\\ 1,t_{N-M}\\ ⋮\\ 1,t_N\end{array}\right].\end{array}$$

(3)

From well-known OLS theory [16] follows that the most likely value of ${[n,k]}^{\top }$, denoted as $\mathit{Z}$, is

$$\begin{array}{c}\hfill \mathrm{E}\left[\begin{array}{c}n\\ k\end{array}\right]=\mathit{Z}={\left({\mathit{F}}^{\top }\mathit{F}\right)}^{-1}{\mathit{F}}^{\top }\mathit{Y},\end{array}$$

(4)

the uncertainty of ${[n,k]}^{\top }$, denoted as $\mathit{R}$, is

$$\begin{array}{c}\hfill \mathrm{Cov}\left[\begin{array}{c}n\\ k\end{array}\right]=\mathit{R}={\left({\mathit{F}}^{\top }\mathit{F}\right)}^{-1}{\sigma }_{\eta }^2,\end{array}$$

(5)

and the noise level is

$$\begin{array}{c}\hfill {\sigma }_{\eta }^2=\frac{1}{M-2}{(\mathit{Y}-\mathit{F}\mathit{Z})}^{\top }(\mathit{Y}-\mathit{F}\mathit{Z}).\end{array}$$

(6)

2.2. Distribution of the Ratio of Two Jointly Normal Variables

Again, the degradation with time, denoted as $x\left(t\right)$, is modelled as a linear function with slope k and intercept n. To account for the uncertainty of the said two parameters, we consider them probabilistic. Their joint distribution is taken to be multivariate-normal

$$\begin{array}{c}\hfill \left[\begin{array}{c}n\\ k\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_n\\ {\mu }_k\end{array}\right],\left[\begin{array}{cc}{\sigma }_n^2,& {\rho }_{n,k}{\sigma }_n{\sigma }_k\\ {\rho }_{n,k}{\sigma }_k{\sigma }_n,& {\sigma }_k^2\end{array}\right]\right),\end{array}$$

(7)

where

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mu }_n\\ {\mu }_k\end{array}\right]=\mathit{Z},{\mu }_n,{\mu }_k\in \mathbb{R},\mathrm{and}\end{array}$$

(8)

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\sigma }_n^2,& {\rho }_{n,k}{\sigma }_n{\sigma }_k\\ {\rho }_{n,k}{\sigma }_k{\sigma }_n,& {\sigma }_k^2\end{array}\right]=\mathit{R},{\sigma }_n,{\sigma }_k\in {\mathbb{R}}^{+},{\rho }_{n,k}\in [-1,1].\end{array}$$

(9)

The threshold, denoted as $\alpha$, is also considered uncertain, meaning

$$\begin{array}{c}\hfill \alpha \sim \mathcal{N}\left({\mu }_{\alpha },{\sigma }_{\alpha }^2\right),{\mu }_{\alpha }\in \mathbb{R},{\sigma }_{\alpha }\in {\mathbb{R}}^{+}.\end{array}$$

(10)

So if $\alpha =kt+n$, then t, interpreted as the FHT, is

$$\begin{array}{c}\hfill t=\frac{\alpha -n}{k}=\frac{\zeta }{k},\end{array}$$

(11)

where $\zeta =\alpha -n.$ The joint distribution of $\zeta$ and k is then

$$\begin{array}{c}\hfill \left[\begin{array}{c}\zeta \\ k\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_{\zeta }\\ {\mu }_k\end{array}\right],\left[\begin{array}{cc}{\sigma }_{\zeta }^2,& {\rho }_{\zeta ,k}{\sigma }_{\zeta }{\sigma }_k\\ {\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta },& {\sigma }_k^2\end{array}\right]\right),\end{array}$$

(12)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\mu }_{\zeta }={\mu }_{\alpha }-{\mu }_n,\\ {\sigma }_{\zeta }^2={\sigma }_{\alpha }^2+{\sigma }_n^2,\\ {\rho }_{\zeta ,k}=\frac{-{\rho }_{n,k}{\sigma }_n}{{\sigma }_{\zeta }}.\end{array}\end{array}$$

(13)

It is clear that $t=\zeta /k$ is a ratio distribution of two (possibly correlated) normal variables. Fortunately, the probability density function for such distributions is available in closed form [17]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_{\mathrm{FHT}}\left(t\right)& =\frac{{\sigma }_k{\sigma }_{\zeta }{\left(1-{\rho }_{\zeta ,k}^2\right)}^{1/2}}{\pi \left({\sigma }_k^2{t}^2-2{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }t+{\sigma }_{\zeta }^2\right)}exp\left[-\frac{1}{2\left(1-{\rho }_{\zeta ,k}^2\right)}\left(\frac{{\mu }_k^2}{{\sigma }_k^2}-2{\rho }_{\zeta ,k}\frac{{\mu }_k}{{\sigma }_k}\frac{{\mu }_{\zeta }}{{\sigma }_{\zeta }}+\frac{{\mu }_{\zeta }^2}{{\sigma }_{\zeta }^2}\right)\right]\hfill \\ \hfill & +\frac{{\mu }_k{\sigma }_{\zeta }^2-{\mu }_{\zeta }{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }+\left({\mu }_{\zeta }{\sigma }_k^2-{\mu }_k{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }\right)t}{\sqrt{2\pi }{\left({\sigma }_k^2{t}^2-2{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }t+{\sigma }_{\zeta }^2\right)}^{3/2}}\hfill \\ \hfill & \times exp\left(-\frac{{\left({\mu }_{\zeta }-{\mu }_{k}t\right)}^2}{2\left({\sigma }_k^2{t}^2-2{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }t+{\sigma }_{\zeta }^2\right)}\right)\hfill \\ \hfill & \times \left[1-2Q\left(\frac{{\mu }_k{\sigma }_{\zeta }^2-{\mu }_{\zeta }{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }+\left({\mu }_{\zeta }{\sigma }_k^2-{\mu }_k{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }\right)t}{{\sigma }_k{\sigma }_{\zeta }{\left(1-{\rho }_{\zeta ,k}^2\right)}^{1/2}{\left({\sigma }_k^2{t}^2-2{\rho }_{\zeta ,k}{\sigma }_k{\sigma }_{\zeta }t+{\sigma }_{\zeta }^2\right)}^{1/2}}\right)\right],\hfill \end{array}\end{array}$$

(14)

where

$$Q\left(x\right)=\frac{1}{2}\left(error\left(\frac{x}{\sqrt{2}}\right)\right).$$

The error is the complementary error function, defined as $erfc\left(z\right)=1-erf\left(z\right)$, where $erf\left(z\right)$ stands for the error function

$$error\left(z\right)=\frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-t^2}dt.$$

Equation (14) can then be used to calculate the probability density function of the FHT distribution to obtain the estimates for RUL in a closed-form manner. A full derivation is available as supplementary material on this link.

The general idea of the method is showcased in Figure 1, highlighting the use of a data window to correctly predict the distribution of FHT.

It is important to stress that in this paper it is tacitly assumed that the indicator increases together with the degradation level, i.e. it should have a positive trend. In cases of a high noise level combined with a short time window, the slope ${\mu }_k$ of the estimated local trend model at time $t_N$ might turn negative. In that case the intersection (11) occurs at $t<t_N$, which is not a feasible solution and is therefore discarded.

3. Performance Analysis on Simulated Time Series

To demonstrate the performance of the method above, four different simulation scenarios are analysed.

First, a noisy time series with a fixed linear trend in the latent variable $x\left(t\right)$ is presented.

In the second case, the latent variable undergoes an abrupt change in the trend. This case serves to point out the impact of the window length on the adaptivity of the local linear trend model.

In the third case the latent variable suffers from oscillations, modelled by an autoregressive moving-average (ARMA) [18] process. This example aims to demonstrate the ability of the proposed method to account for the additional uncertainty in the FHT caused by the oscillations in the latent variable.

The fourth case is similar to the third one with one difference, that is, the parameters of the latent model are subjected to an abrupt change, once again illuminating the importance of window length.

3.1. Case 1: Time series with a fixed linear trend

In this simplest case the latent variable $x\left(t\right)$ exhibits a fixed linear trend, while the measured health indicator $y\left(t\right)$ is corrupted with noise

$$\begin{array}{c}\hfill x\left(t\right)=n+kt,n,k\in \mathbb{R}\left(\mathrm{Latent}\mathrm{health}\mathrm{metric}\right),\end{array}$$

(15)

$$\begin{array}{c}\hfill y\left(t\right)=x\left(t\right)+\xi \left(t\right),\xi \sim \mathcal{N}(0,{\sigma }_{\xi }^2),{\sigma }_{\xi }\in {\mathbb{R}}^{+}\left(\mathrm{Observed}\mathrm{health}\mathrm{metric}\right).\end{array}$$

(16)

The parameters used in the simulation are given in Table 1.

The simulated time series with a fixed linear trend is plotted in Figure 2. The time between consecutive samples is $t_{i+1}-t_i=1$h. At time $t_N=400$h several local linear trend models are estimated over data windows $M\in \{400,200,60\}$ and then used to evaluate the distribution of FHT. It is not surprising that utilisation of all the available data provides the most accurate prediction of FHT.

3.2. Case 2: Time Series with Linear Trend Subjected to Abrupt Changes

Like case 1, the latent variable follows a linear function, but its parameters undergo an abrupt change

$$\begin{array}{c}\hfill x\left(t\right)=\left(1-\varphi \left(t\right)\right)(k_{1}t+n_1)+\varphi \left(t\right)(k_{2}t+n_2),k_1,k_2,n_1,n_2\in \mathbb{R}(\mathrm{Latent}\mathrm{health}\mathrm{metric}),\end{array}$$

(17)

$$\begin{array}{c}\hfill \varphi \left(t\right)=\frac{1}{1+exp\left(-\kappa (t-\tau )\right)},\kappa \in {\mathbb{R}}^{+},\tau \in \mathbb{R},\end{array}$$

(18)

$$\begin{array}{c}\hfill n_2=(k_1-k_2)\tau +n_1,\end{array}$$

(19)

$$\begin{array}{c}\hfill y\left(t\right)=x\left(t\right)+\xi \left(t\right),\xi \sim \mathcal{N}(0,{\sigma }_{\xi }^2),{\sigma }_{\xi }\in {\mathbb{R}}^{+}(\mathrm{Observed}\mathrm{health}\mathrm{metric}),\end{array}$$

(20)

where $\varphi \left(t\right)$ is the sigmoid function which changes the system from $x\left(t\right)=k_{1}t+n_2$ to $x\left(t\right)=k_{2}t+n_2$. With the parameters $\tau$ and $\kappa$ one defines the time of transition and the rate of transition respectively. Note that $n_2$ is not a free parameter to ensure a more natural-looking transition. The parameters used to simulate the second example are provided in Table 1.

In Figure 3 the simulated time series is presented along with the predicted distributions of the FHT obtained from local models estimated on three different windows. Shorter windows better adapt to the changes in the trend. In turn, the predicted FHT are more accurate compared to the models resulting from the longer time window, up to a point.

3.3. Case 3: Time Series Resulting from ARMA Process with Fixed Parameters

A more general case of a time series containing a linear trend and an oscillatory component is analyzed $x\left(t\right)=x_{\mathrm{trend}}\left(t_j\right)+x_{\mathrm{osc}.}\left(t_j\right)$. For that purpose ARMA$(p,q)$ is utilized

$$\begin{array}{c}\hfill x\left(t_j\right)=x_{\mathrm{trend}}\left(t_j\right)+x_{\mathrm{osc}.}\left(t_j\right)(\mathrm{Latent}\mathrm{health}\mathrm{metric}),\end{array}$$

(21)

$$\begin{array}{c}\hfill x_{\mathrm{trend}}\left(t_j\right)=l{t}_j+m,l,m\in \mathbb{R},\end{array}$$

(22)

$$\begin{array}{c}\hfill x_{\mathrm{osc}.}\left(t_j\right)=\sum _{i=1}^p{\varphi }_i{x}_{\mathrm{osc}.}\left(t_{j-i}\right)+\sum _{i=1}^q{\theta }_iϵ\left(t_{j-i}\right)+ϵ\left(t_j\right),ϵ\sim \mathcal{N}(0,{\sigma }_{ϵ}^2),{\varphi }_i,{\theta }_i\in \mathbb{R},{\sigma }_{ϵ}\in {\mathbb{R}}^{+},\end{array}$$

(23)

$$\begin{array}{c}\hfill y\left(t_j\right)=x\left(t_j\right)+\xi \left(t_j\right),\xi \left(t_j\right)\sim \mathcal{N}(0,{\sigma }_{\xi }^2),{\sigma }_{\xi }\in {\mathbb{R}}^{+}(\mathrm{Observed}\mathrm{health}\mathrm{metric}).\end{array}$$

(24)

Following the same steps as in Section 3.1 yields FHT distributions that are too narrow, since the oscillations of the latent variable are falsely regarded as white measurement noise. To circumvent this problem, we note that the scale of oscillations $\mathrm{E}\left[x_{\mathrm{osc}.}^2\right]$ is approximately the same as the scale of the perceived white noise ${\sigma }_{\eta }$. We can now transfer this perceived uncertainty, caused by oscillations, to uncertainty of the threshold, so ${\sigma }_{\alpha }={\sigma }_{\eta }$ (see Eq. 10). Such adaptation of the problem statement rescales the distribution of FHT to a much better width, as shown in Figure 4. Note that this approximation is meaningful only under conditions of relatively high signal-to-noise ratios. Otherwise, the uncertainty of the resulting FHT distributions may be overestimated. Alternatively, if the actual noise level ${\sigma }_{\xi }$ is known, one can subtract the actual noise level from the perceived white noise level, so ${\sigma }_{\alpha }=\sqrt{{\sigma }_{\eta }^2-{\sigma }_{\xi }^2}.$

To adequately test the algorithm, three different ARMA models are employed to simulate oscillating latent variable. Their parameters are presented in Table 3.

Their typical behavior is showcased in Figure 5.

In Figure 6, the algorithm was tested using three data windows of different sizes. For a meaningful comparison to the resulting distributions, the ground-truth FHT distribution was approximated by MC sampling of the known ARMA process and compared to the results.

Examining the same situation from a different standpoint, Figure 7 features the same latent processes as Figure 6, but instead the quantiles of the resulting FHT distributions from various sizes of the data-windows are presented.

It turns out that with enough data and far enough away from the threshold the algorithm stands as a viable alternative to a full regression of an ARMA model for FHT prediction. It is also worth noting that in this case, where the time series has a constant trend, retaining as much information as possible leads to more accurate predictions.

3.4. Case 4: Time Series Resulting from ARMA Process with Abruptly Changing Parameters

In dynamical situations, an important feature of the prediction algorithm is the flexibility to quickly adapt to changes in trends. For that purpose we again utilize limited-size time windows. Similarly to Section 3.2, the algorithm is tested on data generated by ARMA process featuring a change in the trend. The model parameters used for simulation are shown in Table 4.

In Figure 8 and Figure 9 the simulated ARMA process with a change in trend is depicted for two examples. The first example has the data from $t=[0,1000]$ h available and for the second one we use $t=[0,1800]$ h. The change in trend is initiated around $t=500$ h for both examples. For the data interval up to about $t=500$ h the slope of the trend is $0.02$, which changes to $0.03$ afterwards. Resulting FHT distributions derived from different data windows are illustrated in Figure 8 and Figure 9. An obvious improvement in the quality of prognosis can be observed after incorporating only the data after the trend change.

In Figure 10 and Figure 11 we show quantiles of FHT distributions from local trend models obtained from differently sized windows on the data. For example, a window $[800,1000]$h includes $20\%$ of the data up to 1000h. The consequence of such flexibility can be seen in frequent changes in the mean of the FHT distribution acquired from the most recent data, where the trend is not yet obvious. For comparison, ground-truth FHT distribution (dark blue curve) is also provided.

We notice that up to the point, adding data improves precision and accuracy. If we also include data before the trend change, the precision of the predictions does not suffer as much as the accuracy. This highlights a potential limitation in assessing the reliability of the prediction. Without the variance of the prediction ballooning, this could lead to a false sense of confidence in the prediction by the operator.

We also notice a slight bias in predictions, when the time of prediction is somewhat close to the EOL, see Figure 9. Due to the limited time available, any large enough oscillation significantly impacts the point of threshold crossing. The latent variable is therefore unable to regress to the mean before the crossing, which significantly alters the ground truth distribution.

4. Application to the Prognosis of a Laboratory SOFC System

4.1. System Description

The study was done on a 10kW SOFC system, which consists of two interconnected modules, the balance of plant (BoP) and the stack. The latter consists of 80 fuel cells and is fuelled by methane. Water-gas shift (WGS), steam methane reforming (SMR) and hydrogen-oxidation reactions take place inside to produce heat, water and electrical power. Air is fed by the blower to the cathode inlet of the stack. Methane is passed through a reformer in order to get hydrogen, which is then delivered to the anode inlet. The remaining fuel in the cathode output is recycled. During the test, various operational parameters were manipulated for a complete characterization of the system.

4.2. Lumped Model of the Stack

The model of the stack is based on its energy balance [3]

$${\mathrm{K}}_{\mathrm{s}}\frac{d{T}_{\mathrm{out}}}{dt}={\dot{E}}_{\mathrm{in}}\left(T_{\mathrm{in}}\right)-{\dot{E}}_{\mathrm{out}}\left(T_{\mathrm{out}}\right)-IU\left(T_{\mathrm{out}}\right),$$

(25)

where $T_{\mathrm{out}}$ and $T_{\mathrm{in}}$ are the outlet and inlet gas temperatures, ${\dot{E}}_{\mathrm{in}}$ and ${\dot{E}}_{\mathrm{out}}$ denote energy flows in and out of the stack, ${\mathrm{K}}_{\mathrm{s}}$ denotes lumped heat capacity of the stack, U is stack voltage and I is load current.

The inlet/outlet energy flows $\dot{E}$ are defined as [3]

$${\dot{E}}_{\xi }=\sum _i{\dot{n}}_i^{\xi }{h}_i\left(T^{\xi }\right),$$

(26)

where ${\dot{n}}_i^{\xi }$ denotes the molar flow of the $i^{\mathrm{th}}$ gas component, $i\in \{$H², H₂O, CH₄, CO, CO₂, O₂, N₂}, at corresponding locations where $\xi \in \{\mathrm{in},\mathrm{out}\}$ i.e. stack inlet/outlet, and corresponding temperatures $T_i^{\xi }$. $h_i$ is the enthalpy of the gas.

For the purpose of constructing such a simple model, SMR, WGS and hydrogen oxidation reactions are considered to be in thermodynamic equilibrium. This assumption is often employed in SOFC systems and can be justified [3,19,20,21,22]. Assuming total methane reforming and using the WGS equilibrium constant from [23], the reforming equilibrium can be solved [24], yielding the concentration of each outlet gas component, which we take to be proportional to its flow rates.

The stack voltage U can be expressed as

$$U={\mathrm{N}}_0\left(U_0-\mathrm{ASR}\times \frac{I}{\mathrm{A}}\right),$$

(27)

where A is the active area of a single cell, ${\mathrm{N}}_0$ is the number of cells and $U_0$ is the Nerst voltage [25]

$$\begin{array}{c}\hfill U_0=1.2586-0.000252T+\frac{\mathrm{R}T}{2\mathrm{F}}ln\frac{p_{{\mathrm{H}}_2}{p}_{{\mathrm{O}}_2}^{0.5}}{p_{{\mathrm{H}}_2\mathrm{O}}},\end{array}$$

(28)

where R is the universal gas constant.

For an approximation of the partial pressure of gases inside the stack, the average of input and output molar flows is used

$$\begin{array}{c}\hfill p_j=\frac{1}{2}\left(\frac{{\dot{n}}_j^{\mathrm{in}}}{{\dot{n}}_{\mathrm{tot}.}^{\mathrm{in}}}+\frac{{\dot{n}}_j^{\mathrm{out}}}{{\dot{n}}_{\mathrm{tot}.}^{\mathrm{out}}}\right),\end{array}$$

(29)

where $j\in \left\{\begin{array}{c}\hfill {\mathrm{H}}_2,{\mathrm{H}}_2\mathrm{O},{\mathrm{O}}_2\end{array}\right\}$ and ${\dot{n}}_{\mathrm{tot}.}^{\mathrm{in}/\mathrm{out}}$ are total anode or cathode inlet/outlet molar flows.

The temperature variation of lumped Ohmic ASR of a single cell can be expressed by the second-order Steinhart–Hart equation [20,25]

$$\mathrm{ASR}={\mathrm{ASR}}_{0}exp\left[\frac{30000}{\mathrm{R}}\left(\frac{1}{T}-\frac{1}{T_0}\right)\right],$$

(30)

where ${\mathrm{ASR}}_0$ is ASR at reference temperature $T_0$. It is ${\mathrm{ASR}}_0$ we are actually after, since we are in a need of temperature-invariant health index.

Putting all the necessary parts together, we arrive at the full lumped dynamic model of the stack

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{K}}_{\mathrm{s}}\frac{dT}{dt}& ={\dot{E}}_{\mathrm{in}}\left(T_{\mathrm{in}}\right)-{\dot{E}}_{\mathrm{out}}\left(T\right)-\hfill \\ \hfill & -I{\mathrm{N}}_0\left(1.2586-0.000252T+\frac{\mathrm{R}T}{2\mathrm{F}}ln\frac{p_{{\mathrm{H}}_2}{p}_{{\mathrm{O}}_2}^{0.5}}{p_{{\mathrm{H}}_2\mathrm{O}}}-{\mathrm{ASR}}_{0}exp\left[\frac{30000}{\mathrm{R}}\left(\frac{1}{T}-\frac{1}{T_0}\right)\right]\times \frac{I}{\mathrm{A}}\right).\hfill \end{array}\end{array}$$

(31)

Since it was only possible to measure ${\mathrm{ASR}}_0$ indirectly and (31) is non-linear, an unscented Kalman filter [26], highly suitable to such problems, was used to filter out as much noise from the measurements as possible.

Lastly, an adequate criterion for the EOL for the stack is needed. While we will be using some rather arbitrary value of ${\mathrm{ASR}}_0$ that signals the EOL, to get a physically meaningful value, one needs to define the lowest appropriate stack voltage at the same typical operating point. This minimum voltage can be easily converted to an ${\mathrm{ASR}}_0$ value using (27).

4.3. Prognostics of the Remaining Useful Life Based on ASR as a Health Indicator

In Figure 12 operational data of almost 1850h is presented. The system was operated in a non-stationary regime during which there were frequent changes in the load. Therefore, the stack voltage is not the best candidate for the health indicator. Instead, ASR was chosen as the health metric since it is invariant to changes in process variables (technically, we have chosen ASR₀ as the metric according to (31)). For validation of the algorithm the threshold of ASR was set to $0.11\mathsf{\Omega }$mm² (a red horizontal line in Figure 12).

Validation of the Proposed Method by Comparison to an Alternative Prognostic Algorithm

To assess the new algorithm on experimental data, where the latent model is unknown, we decided to compare it to an ARMA process whose parameters were obtained by the maximum-likelihood method. Specifically, we have used statsmodels v0.14.1 [27], a Python package for inference for statistical models. After comparing several different variations of the model via the Akaike information criterion, we have landed on ARMA$(2,0)$ with a deterministic linear function as an exogenous variable (same framework as in (21)).

4.4. Results

Results for the experimental data can be seen in Figure 13. The improvement in the quality of resulting FHT distribution is proportional to the number of data points included in the prognosis, which is expected since the ASR time signal includes no obvious trend changes. The alternative FHT distribution, sampled from fitted ARMA model, can be seen in Figure 13 as well. It was smoothed using kernel density estimation (KDE) [28] for easier comparison.

To further compare the quality of the derived algorithm to an ARMA model, Figure 14 shows the evolution of variance and mean of both predictive distributions over the amount of data taken into account. The time of perceived EOL is shown with a horizontal red line.

Again, the presented algorithm shows comparable results to the alternative. The larger variance at the start can be attributed to the lack of sufficient data. Such uncertainty in the prediction can be valuable for the operator, as it suggests a need for caution when interpreting the prognosis. Over time we can observe a slow regression of the mean to the ground truth FHT, together with the variance shrinking to a useful level.

On the other hand, the alternative prediction from an ARMA model suffers from overconfidence issues in the beginning, since only a point-wise fit is considered for the parameters. Other than that, it also slowly converges to the vicinity of the actual FHT.

Similarly to Figure 14, Figure 15 shows quantiles of FHT distributions derived throughout the experiment. Abrupt changes in the mean of the FHT distributions can be observed, which indicate volatile oscillations of ASR during the measurement process. In the magnified inset of the Figure we again compare the actual FHT to the prediction. Upon closer examination, we conclude that the predicted FHT distributions successfully encompass the eventual EOL as early as $t=800$h.

5. Conclusion

In this paper, a novel algorithm for predicting the remaining useful life for SOFC systems based on adaptive trend estimation of the time-evolving health indicator is proposed. A distinctive feature of the algorithm is that the predicted distribution of the remaining useful life is provided in a closed form. Hence the computational time is insignificant, thus rendering the algorithm appropriate for implementation in industrial low-cost platforms with limited computational resources.

The performance of the algorithm is first demonstrated on four different artificial types of time series that are intended to mimic the potential evolutions of the health indicators met in practice.

The study shows that in the simple case of a constant trend time series larger size data windows deliver more accurate predictions of FHT. If on the other hand time series exhibit variable trends, shorter data windows offer better flexibility to account for the variations.

Tests on time series that depart from the strongly monotone features due to internal dynamics modelled by ARMA process show that the proposed algorithm consistently captures the distribution of FHT without the need for inferring the coefficients of the ARMA process.

Finally, the performance of the proposed algorithm is compared with a more elaborated algorithm ARMA$(2,0)$ on experimental data. Specifically, ASR measurements from a laboratory SOFC stack were modelled as a time series. It turns out that the proposed algorithm consistently delivers the same level of quality of prognosis throughout the experiment even though estimation of ARMA model parameters is computationally more demanding than operating with linear local trend.

To sum up, the proposed prognostic algorithm is usable as a plug-and-play tool appropriate for SOFC applications lacking specific prognostic models and without prior knowledge of the degradation. In the follow up it would be of interest to explore adaptive strategies to dynamically adjust the threshold values based on evolving system conditions. This adaptability could enhance the algorithm’s performance in handling varying degradation patterns and unexpected changes in the health indicator.