A Transferable Multi-stage Model with
Cycling Discrepancy Learning for
Lithium-ion Battery State of Health Estimation

Performance fading of LiBs includes multiple degradation behaviors over the whole degradation procedure, which is recognized as a multi-stage characteristic [13]. Typically, Jacqueline et al. [14] divides the life cycle into three main stages with the analysis of solid-electrolyte interphase: acceleration, stabilization, and saturation. The concept of stage describes the consecutive discharging cycles following the same degradation pattern. The degradation pattern maintains consistency within the same stage, but is dissimilar between different stages. This description is used to explain the transferable multi-stage idea, and specific stage division and online classification are beyond the scope of this work but can be found in the previous work [13]. Given there are $C$ stages in LiBs, the whole cycling data can be sequentially assigned into each stage accordingly. A three-dimensional data matrix ${\overline{\text{X}}}_c(I_c\times J\times K)$ is arranged for the cycling data collected in the $c^{th}$ stage, where $I_c$ is the number of cycles in the $c^{th}$ stage, $J$ is the number of variables, and $K$ is the length of samplings in each cycle. Correspondingly, the actual capacities of the battery in the $c^{th}$ stage are measured as $Q_c\in R^{1\times I_C}$.

Fig. 1 illustrates the degraded capacity with varying speeds provided by the National Aeronautics and Space Administration (NASA) [15]. Essentially, the variable capacity is the external presentation of the inherent discharging signals. According to a prior knowledge [13], the whole life span is divided into three stages, ranging from the healthy stage and slow-degradation stage to the rapid-degradation stage.

SOH is generally calculated by actual capacity over nominal capacity during cycling. Consequently, timely replacement of a LiB is highly recommended when the actual capacity is below a failure threshold, e.g., $70\%$ of its rated capacity [16]. Typically, a fresh LiB is believed to be 100$\%$ healthy, and its health level will decline over cycles. Considering the stage information, the index $H_{c,i}$ is defined to indicate SOH in the $c^{th}$ stage as below,

where $Q_{c,i}$ is current capacity of the $i^{th}$ cycle in the $c^{th}$ stage, and $Q$ is rated capacity for a battery.

Source domain includes two critical parts: feature space $\chi$ and a marginal probability distribution $P_S\left(X_s\right)$ about samples $X_s$ from this domain. Given this, the purpose of source task is to find a regression $f_{TL}(\cdot )$ and a feature subspace $\chi$ at a certain stage from the source battery, which minimizes the error $E_s$ between the real capacities $Q$ and predictions given below,

With $\chi$ and $f_{TL}(\cdot )$ from the source task, TL attempts to minimize the estimation error for cycle $x_T$ from the target task. However, difficulty lies in that marginal probability distribution $P_T\left(x_T\right)$ of target data is always different from $P_S\left(X_s\right)$. Meanwhile, the amount of data in the target domain is insufficient to well-describe $P_T\left(x_T\right)$. $f_{TL}(\cdot )$ derived in Eq. (2) expects to have a broad generalization ability, which could be used for the target battery at the same stage below,

Although cycling data like discharging voltage look similar, the available time-span is ever-decreasing over cycles because of the unavoidable influence of physicochemical reactions. In light of this, we use cycling consistency among cycles to stand for the similarity and the cycling discrepancy over cycles to indicate the distinctions. Fig. 2 visualizes the basic idea of cycling discrepancy analysis with TL in a three-dimension view, including capacity, stage, and battery.

It is worth pointing out that an assumption lies in the division of the number of stages for both source and target batteries. This assumption is made from the analysis of battery degradation mechanism [17], as many batteries typically follow three states: health status, slight degradation, and rapid degradation. To delivery a coherent framework applied for most batteries, three stages is sufficient in this case.

Phase subspace theory originated from chaotic system analysis, and it is defined as a subspace where all possible states of the system are covered. The most important role of phase subspace allows for analyzing unseen system dynamics through a single observable time series. In this way, a specific state of the system can be represented by a point in the phase space and the time evolution of the system creates a trajectory in the phase space [18]. As for the battery charging and discharging procedure, as shown in Fig. 4, phase subspace aims at describing the mechanics of lithium-ions through observable current, voltage, and temperature, linking the microstates to macrostates. With the flow of lithium-ions in microstates, the current and voltage are generated in macrostates. Although the clear description of the electrochemical reaction is difficult and even impossible, we apply phase space to macroscopic states of the degradation procedure by utilizing discharge current, discharge voltage, and environment temperature. A point in this phase space corresponds to a discharging cycle, with different available capacities as the initial condition.

This part reconstructs raw cycling data grouped into each stage in a high-dimension space. In practice, discharge current, discharge voltage, and environment temperature are easily measured and widely used to infer battery capacity from a macroscopic perspective. LiBs involve physicochemical reactions that have complex microscopic properties, and it is impossible to reveal full dynamics only using current, voltage, and temperature. A solution is to describe the deep dynamics of a low-dimension system by exploiting a high-dimension embedding. This embedding is expected to uncover unknown and valuable information via the partially observable sensors. Phase subspace reconstruction (PSR) meets this purpose well, which is designed to reveal the underlying dynamics of a system with partial low-dimension measurements [19]. In addition, data reconstruction can graphically detect hidden patterns and dissimilarities across the cycling data.

With dimensions of reconstructed space $r$ and time delay $\tau$, $X_i$ is extended to a high-dimension matrix below,

where,

and $X_i$ is the $i^{th}$ cycle of ${\overline{\text{X}}}_c$; $V_i$ is the reconstructed data for voltage; the reconstructed current $I_i$ and temperature $T_i$ can be obtained in a similar way.

$r$ influences the dimension of reconstructed space, which can be computed using the false neighbor method [19]. A proper $\tau$ is determined when mutual information between $x$ and its lagged vector is minimal. Performing Eq. (4) to each degradation cycle in the same stage, we will obtain the number of $I_C$ reconstructed cycling data using PSR. For clarity, the reconstructed data are still denoted as $X_1$, $X_2$, to $X_C$ corresponding to Cycles 1 to $I_C$, and the dimension of the reconstructed space is still denoted as $J$.

According to the above analysis, $X_i$ can be the superposition of cycling consistency $S_{i,s}$ and cycling discrepancy $S_{i,d}$ when a linear mapping relationship is imposed. Through separating $S_{i,d}$ from $S_{i,s}$, invertible matrices $P_s$ and $P_d$ are corresponding subspaces, to derive the cycling discrepancy and consistency components. Then we have,

Essentially, this problem aims to find the cycling discrepancy features $S_{i,d}$. However, $S_{i,d}$ is always non-stationary along cycle-wise direction. As such, $P_d$ is subject to the nonunique solutions when Eq. (6) is directly solved. An alternative solution is to locate $S_{i,s}$, which almost keeps consistent among cycling data. Consequently, the cycling discrepancy components can be derived by projecting the cycling data onto the orthogonal subspace. The cycling consistency components in Cycle $i$ are decomposed according to Eq. (6) below,

Projecting $X_1$ through $X_{I_C}$ onto the space $P$, a series of features will be computed as $S_1,S_2,\cdots ,S_{I_C}$. Since it is a considerable burden to evaluate the consistency between these features one by one, an efficient way is to employ Gaussian distribution with zero mean and unit variance as an auxiliary indicator, a typical stationary signal. Then KLD is employed to measure the distribution divergence between $S_{i,s}$ and Gaussian distribution. As such, the optimal projection matrix $P_s$ corresponds to the minimal KLD.

With Z-score normalization, the mean and variance of $X_i$ are calculated as $u_i$ and ${\Sigma }_i$, respectively. Introducing a whitening matrix $W_c$ as a part of $P$, the cycling consistency components are obtained below,

where $\Pi$ is the updated projection matrix to be solved, $I_s$ is the identity matrix $I(J\times J)$ truncated to the first $S$ rows, and $W_c$ is calculated by $\sqrt{{\tilde{X}}_c^T{\tilde{X}}_c}$, where ${\tilde{X}}_c$ is the average of all cycling data in the $c^{th}$ stage.

By performing linear transformation $I_s{\Pi W}_c{X_i}^T$ on each cycling data, solving Eq. (8) is transformed to an optimization problem by minimizing the objective function as follows,

where ${\tilde{u}}_{i,s}=I_s{\Pi W}_c{u}_i$ and ${\tilde{\Sigma }}_{i,s}=I_s{\Pi W}_c{\Sigma }_i(I_s{\Pi W}_c{)}^T$ are the mean and variance after linear transformation, respectively; $D_{KL}(\cdot )$ is the KLD function; $N(0,I)$ is Gaussian distribution with zero mean and unit variance; $P_{i,s}\left(m\right)$ is the probability density of $N({\tilde{u}}_{i,s},{\tilde{\Sigma }}_{i,s})$ for the $m^{th}$ sample; $Q\left(m\right)$ is the probability density of $N(0,I)$ for the $m^{th}$ sample.

The matrix $\Pi$ is solved by conjugate gradient descend [13]. Cycling consistency components with dimensions $S\times K$ are derived below,

Consequently, the cycling discrepancy components with dimensions $(J-S)\times K$ are spanned by the remaining subspace $(\Pi -{\Pi }_s)$. Then the components $S_{i,d}$ in Cycle $i$ can be estimated as follows,

LSTM has been widely applied in various industrial fields for sequential learning, such as risk prediction [20] and performance prediction [21]. Its temporal learning ability is ensured through elaborately designed LSTM cells, which possess memory capability and overcome the vanishing gradient problem.

Taking $S_{i,d}$ as the input, LSTM layers integrate and compress cycling data into the final hidden state for regression. An essential of the LSTM cell is the introduction of the memory and forget mechanism. By storing historical information into the hidden state $h_{k-1}$, the candidate memory mechanism removes useless information for regression and retains the essential part. Next, $h_k$ at time $k$ is iteratively updated with the information from input gate, forget gate, and output gate at the next time, and then essential information contained in $S_{i,d}$ will be compressed into the final vector. Let $H_{i,d}$ be the collection of output vectors of the LSTM layer. This flattened information passes through a fully connected layer below,

where $\hat{Q}$ represents the predictions; $W_f$ and $B_f$ represent the transformation weights and bias of fully connected layer, respectively; and ReLu function is usually selected for $F_f(\cdot )$.

A TemCap model is proposed to enable deep feature representation with a three-layer hierarchical structure for stages with sufficient training data. The first two layers focus on feature representation within each cycle. Specifically, to better present the overall data in a cycle, the part-whole relationship is introduced by breaking the entire cycle into a series of local representations. In the third layer, temporal correlations among cycling data are linked with each other. The information density is from low to high, corresponding from the first layer until the third layer.

The designed hierarchical structure fully utilizes information of discharging data within a cycle and over cycles. To achieve this purpose, we use basic capsules to describe the local information and advanced capsules to indicate the global information. A capsule is defined originally as the fundamental element for feature representation, which has been used to describe properties such as pose, texture, or deformation of an entity or a part of an entity [22]. Mathematically, it is a multidimensional vector achieved through a group of neurons, each of which stands for one dimension. The details of each layer and its formulations are described below:

$\bullet$ Basic Capsule Layer (BCL): Basic capsules are designed to capture the local features of cycling data, which could be denoted as a vector $u_m$ for the $m^{th}$ capsule. From the viewpoint of feature representation, capsules are the basic units in the feature space, and their combinations perform high-level features. Projecting $S_{i,d}$ calculated in Eq. (10) into a space $\Xi$, a series of feature maps are calculated below,

where $\ast$ is convolutional operation, $f(\cdot )$ is the ReLu activation function, $b$ is bias, and $U$ is collection of basic capsules.

$\Xi$ determines the quality of generated local features in an iterative learning way. Initially, $\Xi$ can be generated through convolutional neural layers with randomly generated convolution kernels, but a pooling operation is emitted to avoid information loss. A basic capsule is constructed by grouping values at the same position in feature maps with a length $M$. As such, the continuous $M$ kernels correspond to a certain space $\Xi$. A series of basic capsules with dimensions $1\times M$ will be constructed. This layer transforms the scalar-based outputs into vector-based outputs with the capsule concept.

$\bullet$ Advanced Capsule Layer (ACL): Advanced capsules are high-dimensional features conveying stronger representations compared to basic capsules. For instance, for classification problems, advanced capsules allow distinguishing a dog from a cat by simultaneously viewing the orientation of eyes, nose, ears, mouse, etc. The input of an advanced capsule $j$ is calculated by summarizing all basic capsules below,

where $W_{m,j}$ is the projection matrix to extend low-dimension space to high-dimension space, and $c_{m,j}$ is the coefficient.

As the similar expression of basic capsules, we use a neuron to stand for one dimension of an advanced capsule, and a complete advanced capsule consists of several neurons corresponding to the given dimension. Sequential neurons in fully connected layers are grouped to construct advanced capsules. Essentially, we need to find a high-dimensional space spanned by the advanced capsules. Besides, as the integration of basic capsules distributed at a specific time in a discharging cycle, the advanced capsules contain temporal information.

$\bullet$ Temporal Representation Layer (TRL): This layer focuses on acquiring temporal information over cycles, linking advanced capsules within each cycle. It allows to receive capacity fading trend over cycles and estimate the possible values with the advanced capsules. As such, we stack all advanced capsules of the $i^{th}$ cycle into one dimension vector $V_i=[v_1,v_2,\dots ,v_D]$ using the flatten operation. Then the collection of advanced capsules in each cycle is sequentially fed into LSTM cells-based regression model. For clarity, the details of LSTM are not repeated here.

$\bullet$ Weight Learning for TemCap: The network parameters include projection subspace $\Xi$ in the basic capsule layer, projection subspace in the advanced capsule layer, LSTM cells, and the connected links between each layer. The proposed TemCap is trained iteratively, which integrates both dynamic routing mechanism and gradient descent to update the network’s weights and bias.

We employ the classical dynamic routing [23] to update the links between basic capsules and advanced capsules below,

where $s_j$ is the output of $j^{th}$ advanced capsule after updating, and $v_j$ has the same meaning as that in Eq. (14).

From Eq. (15), $v_j$ with a small length will shrink to slightly above zero. In contrast, $v_j$ with a large length approaches approximately one. As such, length of the output vector of a capsule indicates the probability of a local feature. Dynamic routing is mainly used for optimizing parameters of basic capsules and advanced capsules. The remaining parameters, including projection subspace to derive basic capsules and advanced capsules, and LSTM layers, can be updated by Gradient descent.

When an online reconstructed cycling matrix $X_{new}$ is available, its belonging to a specific stage is identified by evaluating the stage cycling consistency [13]. Afterwards, cycling consistency components $S_{new,s}$ and cycling discrepancy components $S_{new,d}$ of $X_{new}$ are decomposed below,

Ideally, the well-trained source model expects to gain similar accuracy when it aims at the target task. Considering that the derived cycling discrepancy reflects the differences among LiBs caused by natural reaction and the varying operating conditions, the defined similarity analysis contributes to judge whether the developed source model in a specific stage can be directly applied considering manufacturing differences in the consistent phase. Naturally, this kind of manufacturing difference performs constant behaviors during the capacity fading procedure and influences the consistent features over cycles. The stable behaviors may affect either all stages or partial stages, which inspires us to conduct similarity checking from stage to stage. We establish the control limit for cycling consistency components $S_{i,s}$ from the training dataset, which follows the $F$ distribution below,

where ${\Lambda }_s$ is the covariance matrix of ${\overline{S}}_s$, which is the average of all components $S_{i,s}$ in a certain stage, ${\alpha }_s$ is the significant level (here is 0.05) to derive the 95$\%$ confidence limit, and $N_c$ is the total number of samples in Stage $c$.

If $S_{new,s}$ is covered by the control limit $C{L}_s$, the developed source models could be directly adopted by feeding $S_{new,d}$. Otherwise, we will update the compensation information fromn the target data into the source model.

The updating strategy allows to rapidly transfer the well-developed source model from a battery cell to others, handling the estimation error caused by the steady discrepancy between the source and target domains. A handful of discharging cycles at the beginning of the stage are denoted as $X_{new,1}$ to $X_{new,T}$. The correspondingly actual capacities are remarked as $Q_T$, and then their predictions ${\hat{Q}}_T$ of the original source model can be calculated. The vector of estimation error is computed as $E_T={\hat{Q}}_T-Q_T$.

A shallow neural network with a fully connected layer is adopted to regress $E_T$ on the decomposed consistent components $S_{new,d}$. In a parallel structure, final prediction for the target task is the sum of the prediction from the source model and the error compensation from the shallow neural network.

Remark: The training time complexity of the proposed method is discussed, mainly focusing on two kinds of neural networks. First, for the lightweight LSTM network, the overall time complexity per time step is $o\left(W\right)$ [24], where $W$ is the number of all weights in the network. Second, for TemCap, the time complexity of the basic capsule layer is $o(n_b\times s_b^2\times m_b^2)$ [25], where $n_b$ is the number of filters, $s_b$ is the spatial size of the filter, and $m_b$ is the spatial size of the output feature map. The time complexity of the advanced capsule layer is calculated as $o(M\times n\times m)$, where $M$ is the dimension size of a basic capsule, $n$ is the number of the basic capsule, and $m$ is the number of advanced capsules. With an LSTM layer, the total time complexity during training of TemCap is $o(n_b\times s_b^2\times m_b^2+M\times n\times m+W)$.

Since training cycles of B7 in Stages 1 and 3 are less than 50, a lightweight LSTM is preferred as the backbone to develop source models. The first 70$\%$ cycles of each stage are used as training data, i.e., the first 20 cycles in Stage 1 and the first 40 cycles in Stage 3. Estimation results are tested with the remaining cycles with the well-trained model.

With the network configuration listed in Table II, the specific estimation results for each stage of B7 are plotted in Figs. 8(a), 8(d), and 8(g), including training and testing phases. By looking into the results, we have the following understandings:

• The ground truth of battery capacity is non-stationary, decreasing with the increase of cycles. The predictions of the proposed method are capable of capturing the degradation trend very well in the testing phase, signifying the good generalization ability of the proposed method. Moreover, the regeneration phenomenon of battery capacity can be well predicted.

• The proposed method runs ten times to evaluate the average performance and robustness. The generated estimation is close to its ground truth with a narrow variation band (95$\%$ confidence level).

Using the index $RMSE$, traditional LSTM and gated recurrent unit (GRU) [27] are adopted for quantitative comparison with the same stage division information. The network configurations of LSTM and GRU are listed in Table II for clarity, which both are lightweight models.

Contributing to finding the degradation-related information with the help of phase space reconstruction and cycling discrepancy components, the estimation errors of the proposed method are below $0.5\%$ for B7 in both Stages 1 and 2. For a fair comparison, the counterparts, including GRU and LSTM, are repetitively run ten times to evaluate the average performance with a standard variation. Overall, the proposed method outperforms both GRU and LSTM in Stages 1 and 3 regarding accuracy and robustness, as shown in Table III. Specifically, in Stage 3, the estimation accuracy of GRU and LSTM decrease significantly. Furthermore, two non-stage models treating the entire degradation process in a single behavior are employed for comparison, i.e., LSTM-TL [8] and GPR-CC [26]. Both models take the first 70$\%$ cycles of the whole process as training data and the remaining ones as testing data. As the testing cycles mainly belong to the divided Stage 3, it makes sense to compare the estimation results in Stage 3 of the proposed method with that of LSTM-TL [8] and GPR-CC [26]. As such, the proposed method still outperforms both counterparts with the training data only from Stage 3. It is more reasonable to infer future capacity with historical cycling data belonging to the same stage. We have observed the same conclusions for B6 and B5. Specifically, the estimation results in Stage 3 of the proposed method are much less than that of GPR-CC [26]. Detailed estimation results for B5 and B6 are visualized in Figs. 8(b), 8(h), 8(c), and 8(i).

Since the network structure of LSTM and TemCap is not complex, empirical knowledge is used to configure the hyper-parameter [5] [30], and good estimation results are obtained without many tuning burdens. Alternatively, automatic hyper-parameter tuning is regarded as a promising solution, such as neural architecture search [31], Bayesian optimization [32].

The number of training data of Stage 2 is more than 50 cycles. The designed TemCap is adopted for Stage 2 to capture the deep features, following the same data preparation with the other two stages.

According to the network structure of Stage 2 specified in Table II, the specific estimation results of the proposed method for B7, B6, and B5 are plotted in Figs. 8(d), 8(e), and 8(f), respectively. The tracking trend of predictions meets well with the ground truth. Moreover, for B6 and B5, predictions reach the threshold of the end of life at the actual cycle in the testing phase. As such, degraded batteries can be timely replaced without any delay. Following the network configurations given in Table II, the average estimation results of GRU and LSTM for B7 are 1.58 and 2.57, which are larger than the 1.18 yielded by the proposed method. Besides, their estimation variations are also larger than that of the proposed method. As mentioned earlier, the testing data of LSTM-TL [8] and GPR-CC [26] are mainly distributed in Stage 3. Therefore, a large estimation error of the proposed method in Stage 2 doesnot mean a performance reduction. On the contrary, the proposed method yields even better results for B6 and B5 in Stage 2. The detailed comparisons between the proposed method and its counterparts are summarized in Table III, in which the results marked in blue indicate that the proposed method yields better performance.