Sensitivity of Fractional-Order Recurrent Neural Network with Encoded Physics-Informed Battery Knowledge

Abstract

In the field of state estimation for the lithium-ion battery (LIB), model-based methods (white box) have been developed to explain battery mechanism and data-driven methods (black box) have been designed to learn battery statistics. Both white box methods and black box methods have drawn much attention recently. As the combination of white box and black box, physics-informed machine learning has been investigated by embedding physic laws. For LIB state estimation, this work proposes a fractional-order recurrent neural network (FORNN) encoded with physics-informed battery knowledge. Three aspects of FORNN can be improved by learning certain physics-informed knowledge. Firstly, the fractional-order state feedback is achieved by introducing a fractional-order derivative in a forward propagation process. Secondly, the fractional-order constraint is constructed by a voltage partial derivative equation (PDE) deduced from the battery fractional-order model (FOM). Thirdly, both the fractional-order gradient descent (FOGD) and fractional-order gradient descent with momentum (FOGDm) methods are proposed by introducing a fractional-order gradient in the backpropagation process. For the proposed FORNN, the sensitivity of the added fractional-order parameters are analyzed by experiments under the federal urban driving schedule (FUDS) operation conditions. The experiment results demonstrate that a certain range of every fractional-order parameter can achieve better convergence speed and higher estimation accuracy. On the basis of the sensitivity analysis, the fractional-order parameter tuning rules have been concluded and listed in the discussion part to provide useful references to the parameter tuning of the proposed algorithm.

1. Introduction

As the main energy source and energy storage device, the lithium-ion battery (LIB) and its management have drawn much attention recently. For battery management, safety [1], durability [2], and heterogeneity [3] are the three aspects to be investigated, and the state estimation is the basic function for the three aspects [4,5]. Among the state-of-art literature, the estimation of the state of charge (SOC) [6], state of health (SOH) [7], and remaining useful life (RUL) [8] have been realized by various model-based methods and data-driven methods [9,10]. Model-based methods, especially the electrochemical model [11], can describe inner reactions and lithium concentration changes, such as a loss of lithium inventory (LLI) and loss of active material (LAM) [12,13]. The fractional-order element, also called a Warburg element, has been introduced into the equivalent circuit model (ECM) to construct a fractional-order model (FOM) [14,15]. The FOM can reflect diffusion dynamics in the electrochemical impedance spectrum (EIS), such as the charge transfer reaction and double layer effects in the mid-frequency, and the solid diffusion dynamics in the low-frequency [16,17]. However, the computation burden and complexity of partial derivative equations (PDEs) limit the application of model-based methods in realistic situations. The data-driven method provides a new estimation pattern and replacement for the traditional model-based methods [18]. For realistic battery data with unstable quality [19], data-driven methods are usually designed with a machine learning (ML) algorithm. Among the ML algorithm, neural networks (NNs) are the mainly investigated direction [20], such as the deep neural network (DNN) [21], long short term memory (LSTM) [22], and gated recurrent unit (GRU) [23]. However, the inner battery knowledge cannot be mapped and learned by the state-of-art algorithms due to the lack of interpretability. Hence, current algorithms need to be enhanced to learn the battery’s inner information.

As presented above, model-based methods (white box) and data-driven methods (black box) have been developed independently in two different directions. Hence, physics-informed machine learning has been proposed to combine the advantages of the model-based method and data-driven method, by encoding a physical law of the predicted object. The embedded physics information in ML algorithms were mainly concluded into an observational bias, inductive bias, and learning bias [24]. In this work, we take the advantages of the fractional-order recurrent neural network (FORNN) and fractional-order gradients, which have been investigated theoretically to process time-series data [25,26]. A fractional-order recurrent neural network encoded with physics-informed knowledge of LIB is proposed, called FORNN with PIBatKnow. With encoding battery physical laws, the neural network can be informed by battery knowledge, which acts as a part of the states feedback, the network loss function, and the gradient in backpropagation. The proposed algorithm applies the widely used PDEs to reflect battery physics meanings, which are included in the modeling of LIBs, such as an electrochemical model and ECMs [27,28]. Then, PDEs with physics meanings of the battery are embedded into the network framework; thus, it serves as the representation of battery dynamic characteristic and convergence constraints for network update. Enhanced by physics-informed battery knowledge, nine key fractional-order parameters have been added, and this paper focuses on the sensitivity of the fractional-order parameters to the performance of the proposed FORNN with PIBatKnow. Experiments under dynamic operation conditions are conducted to analyze the sensitivity, and several conclusions have been made in the discussion part to provide useful instructions to the parameter tuning of the proposed algorithm.

2. Preliminaries

2.1. Fractional-Order Derivative

2.1.1. Fractional-Order Derivative Definitions

Fractional calculus is generalized beyond integer calculus by introducing the fractional order $\alpha$, and many definitions of the fractional operator has been developed in the literature. Three commonly used definitions for fractional derivatives are Riemann-Liouville (R-L), Grünwald-Letnikov (G-L), and Caputo [29]. The R-L definition has too complicated initial conditions to calculated in application, thus this work considers the G-L definition and Caputo definition.

Definition 1

(Grünwald-Letnikov (G-L) Derivative [30]).

$${}^{GL}{D}_t^{\alpha }f\left(t\right)=\underset{h\to 0}{lim}{h}^{-\alpha }\sum _{j=0}^{\left[\frac{t-t_0}{h}\right]}{(-1)}^j\left(\begin{array}{c}\alpha \\ j\end{array}\right)f(t-jh),$$

(1)

where ${}^{GL}{D}_t^{\alpha }$ represents the fractional operator in G-L definition, $f\left(t\right)$ is a certain integrable function in $[a,b]$, $\left[\frac{t-t_0}{h}\right]$ represents the approximate recursive terms of integer order parts, $\left(\begin{array}{c}\alpha \\ j\end{array}\right)=\alpha !/j!(\alpha -j)!$ represents the coefficients of the approximate recursive terms.

Lemma 1 (Finite form of Grünwald-Letnikov Derivative).

If the infinite form in (1) is limited as finite terms L, ${}^{GL}{D}_t^{\alpha }$ can be discretized as [31]

$${\Delta }^{\alpha }f\left(x\right)=\sum _{j=0}^L{(-1)}^j\left(\begin{array}{c}\alpha \hfill \\ j\hfill \end{array}\right)f(k-j),$$

(2)

where the sum formula in the right side only selects L finite terms of $f(k-j),j=0,1,\dots ,L$, which are the discrete values in the L previous sampling moments.

Definition 2

(Caputo Derivative [30]).

$${}^C{D}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma (n-\alpha )}\underset{t_0}{\overset{t}{\int }}\frac{f^{\left(n\right)}\left(\tau \right)}{{(t-\tau )}^{\alpha -n+1}}\mathrm{d}\tau ,t>t_0,$$

(3)

where ${}^{\mathrm{C}}{D}_x^{\alpha }$ represents the fractional operator in Caputo definition, $\alpha \in (n-1,n)$, $\Gamma (·)$ is the Gamma function.

Lemma 2.

If $f\left(x\right)$ is smooth with finite fractional-order derivative, Caputo derivative in (3) can be discretized as [32]

$${}_{x_0}^{\mathrm{C}}{D}_x^{\alpha }f\left(x\right)=\sum _{i=n}^{\infty }\frac{f^{\left(i\right)}\left(x_0\right)}{\Gamma (i+1-\alpha )}{(x-x_0)}^{i-\alpha }.$$

(4)

Lemma 3.

In (4), if $f\left(x\right)={(x-x_0)}^{\lambda }$, (4) can be deduced as [33]

$${}_{x_0}^{\mathrm{C}}{D}_x^{\alpha }{(x-x_0)}^{\lambda }=\frac{\Gamma (\lambda +1)}{\Gamma (\lambda -\alpha +1)}{(x-x_0)}^{\lambda -\alpha },$$

(5)

where $\lambda -\alpha +1>0$.

The G-L definition in (2) has a discrete implement form for application, and the Caputo definition in (5) holds a simple form for quadratic function, which can apply the fractional-order chain rule. As more general forms of the derivative than the integer-order one, both the G-L definition and Caputo definition can transfer into the integer-order derivative when $\alpha =1$. Hence, we employ G-L derivative for the fractional-order state feedback and fractional-order constraints, then the Caputo derivative for backpropagation with the fractional-order gradient in this work.

2.1.2. Fractional-Order Gradient

The gradient usually works as the iteration direction to search the optimal minimum point of an unconstrained convex optimization problem [34].

$$\underset{x}{min}f\left(x\right),$$

(6)

where $f\left(x\right)$ is a smooth convex function with a unique global extreme point $x^{*}$. Gradient in continuous form can be presented as

$$\dot{x}=-\rho \nabla f\left(x\right),$$

(7)

where $\rho >0$ is the iteration step size, and $\nabla f\left(x\right)$ is the gradient of $f\left(x\right)$ at x. Equation (7) can be discretized as

$$x_{k+1}=x_k-\rho \nabla f\left(x_k\right),$$

(8)

where k is the iteration step, and $x_k$ is the discrete value at step k. Fractional-order derivative introduces a fractional order $\alpha$ with more flexible performance in fractional-order gradient. The fractional-order gradient in the Caputo definition can be presented as

$$x_{k+1}=x_k-\rho {\nabla }_{x_k}^{\mathrm{C}}{D}_{x_{k+1}}^{\alpha }f\left(x\right),$$

(9)

where $0<\alpha <1$. The convergence of the fractional-order gradient in (9) depends on the fractional order $\alpha$ and the initial value $x_0$, and (9) can converge to a global extreme point $x^{*}$ if the Caputo definition is calculated by the discrete equation (4). In this work, fractional-order gradients are embedded into the backpropagation of FORNN to accelerate the training process and improve prediction accuracy with battery physics information.

2.2. Fractional-Order Equivalent Circuit Model

Compared to the electrochemical model, ECM for LIB holds simpler structure but enough battery physics information for the application to NNs. The fractional order constant phase element (CPE), or called the Warburg element, is introduced due to the capacitance of LIBs in low frequency and mid-frequency, which are not ideal ones with the first-order derivative, which should be modelled as a fractional-order element [30]. Hence, in this work, FOM is employed to describe the three main parts (high-frequency inductive tail, mid-frequency reaction, and low-frequency diffusion dynamics) of the battery EIS [16].

Warburg element or CPE can manifest the phenomenon of capacitance dispersion instead of an ideal capacitor in LIBs [30]. The voltage-current relationship in time domain and the impedance in the frequency domain of a CPE is presented as

$$\left\{\begin{array}{c}i\left(t\right)=C_{CPE}\frac{{\mathrm{d}}^{\alpha }u\left(t\right)}{\mathrm{d}{t}^{\alpha }},0<\alpha <1,t\ge 0\\ Z_{CPE}\left(s\right)=\frac{U\left(s\right)}{I\left(s\right)}=\frac{1}{C_{CPE}·s^{\alpha }}=\frac{1}{C_{CPE}{\left(j\omega \right)}^{\alpha }}\end{array}\right.,$$

(10)

where $Z_{CPE}$ is the complex impedance, $C_{CPE}$ is the capacity coefficient, j is the imaginary unit, $\alpha$ is the fractional order related to capacitance dispersion and $\omega$ is the angular frequency.

Figure 1 presents the mainly used fractional-order ECM in four different forms. If $n=1$ in Figure 1a, ECM in Figure 1a becomes a fractional-order Thevenin model (also called 1-RC model), which only considers high-frequency and mid-frequency reactions, and ignores the Warburg element of diffusion effects [35]. Figure 1b and Figure 1c are fractional-order Partnership for a New Generation of Vehicles (PNGV) and fractional-order “Randles” model, respectively, and both of them are systems with two fractional orders. Fractional-order Randles model in Figure 1c is proposed to mainly reflect the double layer effect in mid-frequency, which shows like a semi-ellipse rather than a semi-circle due to capacitance dispersion, while the fractional-order PNGV model in Figure 1b is a widely used FOM due to the full-scale reflection of LIB dynamics in all frequency range. The corresponding state-space model (SSM) of the fractional-order PNGV model is [16]

$$\begin{array}{c}\left[\begin{array}{c}\frac{d^{{\alpha }_1}{U}_{CP{E}_1}}{d{t}^{{\alpha }_1}}\\ \stackrel{}{\frac{d^{{\alpha }_2}{U}_{CP{E}_2}}{d{t}^{{\alpha }_2}}}\end{array}\right]=\left[\begin{array}{cc}-\frac{1}{R_1{C}_{CP{E}_1}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{U}_{CP{E}_1}\\ U_{CP{E}_2}\end{array}\right]+\left[\begin{array}{c}\frac{1}{C_{CP{E}_1}}\\ \stackrel{}{\frac{1}{C_{CP{E}_2}}}\end{array}\right]I,\hfill \\ U_t=\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{U}_{CP{E}_1}\\ U_{CP{E}_2}\end{array}\right]+R_{\mathrm{ohm}}I+E,\hfill \end{array}$$

(11)

where $U_{CP{E}_2}$ and $U_{CP{E}_1}$ are the voltages of Warburg element $CP{E}_2$ and $CP{E}_1$; $C_{CP{E}_2}$ and $C_{CP{E}_1}$ are the capacitance of $CP{E}_2$ and $CP{E}_1$; ${\alpha }_1$ and ${\alpha }_2$ are the fractional orders of $CP{E}_2$ and $CP{E}_1$, respectively; I is the current. For higher modeling accuracy, some extended high-order FOMs for LIB are proposed by adding low-frequency component and replacing ideal capacitor with CPEs [15], which results in a high-order FOM with three fractional orders, as shown in Figure 1d.

3. Fractional-Order Recurrent Neural Network with Physics-Informed Knowledge

To make algorithm informed by physics and mechanism, the fractional-order derivative is applied to the recurrent neural network to construct the fractional-order recurrent neural network with physics-informed knowledge, simplified as FORNN with PIBatKnow in the following. The integral architecture of FORNN with PIBatKnow is presented in Figure 2. Three aspects of FORNN with PIBatKnow can be improved as fractional-order state feedback for forward propagation, fractional-order constraints for loss function, and fractional-order gradients for backpropagation, which are introduced in every subsection of this part, respectively.

The basis of FORNN with PIBatKnow in Figure 2 is an RNN, which contains an input layer with m neurons (first layer in Figure 2), several hidden layers $(h_1,\dots ,h_{p-1},h_p,\dots ,h_l)$ with $n_p(p=1,2,\dots ,l)$ neurons (middle part in Figure 2), and an output layer with q neurons (last layer in Figure 2), respectively. Suppose the training input dataset is $({\mathbf{x}}_i,{\mathbf{y}}_i),(i=1,2,\dots ,N)$, where ${\mathbf{x}}_i={(x_{i1},x_{i2},\dots ,x_{im})}^{\mathrm{T}}$ is the network input, and ${\mathbf{y}}_i={(y_{i1},y_{i2},\dots ,y_{iq})}^{\mathrm{T}}$ is the network target (training label). To simplify expression, vectors ${\mathbf{x}}_i$ and ${\mathbf{y}}_i$ are presented as x and y in the rest of this paper. The specific part of RNN is the chain structure inside the hidden layers, which hold hidden states feedback in time series. Assume $W_p$ and $b_p$ be the weight and bias matrix connecting the $(p-1)\mathrm{th}$ hidden layer to the $p\mathrm{th}$ hidden layer, $W_{hp}$ be the weights for memory updates of the $p\mathrm{th}$ hidden layer in the chain structure of RNN, $g(·)$ be the activation functions, and $L\left(g\right(x),y)$ be the loss function. Within the pre-set epochs threshold of training process, RNN would go through forward propagation and backpropagation with training data, and the forward propagation starting from the input layer can be presented as

$$\left\{\begin{array}{c}{a}_p\left(x\right)=W_p{h}_{p-1}\left(x\right)+b_p,\hfill \\ h_p\left(x\right)=g\left(a_p\left(x\right)\right),p=1,2,\dots ,l\hfill \end{array}\right.$$

(12)

where $a_p\left(x\right)$ and $h_p\left(x\right)$ are the input and the output of the $p\mathrm{th}$ hidden layer, respectively. Equation (9) is the basic iterative equation of the proposed FORNN with PIBatKnow.

3.1. Fractional-Order State Feedback

Fractional-order derivative can be introduced into the state feedback in the chain structure of the proposed FORNN, resulting in fractional-order state feedback, as shown in Figure 3b. Figure 3a is the integer-order state feedback, which can be presented as

$${}_{t_0}{D}_t{h}_p\left(t\right)=-W_{hp}{h}_p\left(t\right)+{W^{\prime }}_{p}g\left(h_{p-1}\left(t\right)\right)+d_p,p=1,2,\dots ,l$$

(13)

where ${W^{\prime }}_p$ is a coefficient related with weights $W_p$, $d_p$ is a constant related with bias $b_p$, and memory weights $W_{hp}$ are shared among all the moments in the time series.

Fractional-order state feedback in Figure 3b is extended from the integer-order one in Figure 3a, thus fractional-order state feedback can be presented as

$${}_{t_0}{D}_t^{\alpha }{h}_p\left(t\right)=-W_{hp}{h}_p\left(t\right)+{W^{\prime }}_{p}g\left(h_{p-1}\left(t\right)\right)+d_p,p=1,2,\dots ,l$$

(14)

where ${}_{t_0}{D}_t^{\alpha }$ is the fractional-order derivative of $h_p\left(t\right)$ in $[t_0,t]$, and (14) includes the integer-order derivative (13) when $\alpha =1$. Equation (14) can be discretized by the G-L definition in (2), which is

$${\Delta }^{\alpha }{h}_p\left(t\right)=\sum _{j=0}^L{(-1)}^j\left(\begin{array}{c}\alpha \hfill \\ j\hfill \end{array}\right)h_p(k-j)$$

(15)

Take (15) into (14), we obtain that

$$\sum _{j=0}^L{(-1)}^j\left(\begin{array}{c}\alpha \hfill \\ j\hfill \end{array}\right)h_p(k-j)=-W_{hp}\left(k\right)h_p(k-1)+{W^{\prime }}_p\left(k\right)g\left(h_{p-1}\left(k\right)\right)+d_p$$

(16)

Note that $h_p\left(k\right)$ is included in the first item ($j=0$) of (15), and the second item of (15) is $-h_p(k-1)$ when $j=1$. Hence, expanding the first item and the second item in (15), then normalizing the coefficients, the discrete form of fractional-order state feedback can be deduced as

$$\begin{array}{cc}\hfill h_p\left(k\right)=& (1-W_{hp}\left(k\right))h_p(k-1)+\sum _{j=2}^L{(-1)}^{j+1}\left(\begin{array}{c}\alpha \\ j\end{array}\right)h_p(k-j)\hfill \\ \hfill & +{W^{\prime }}_p\left(k\right)g\left(h_{p-1}\left(k\right)\right)+d_p\hfill \end{array}$$

(17)

The upper limit L of the finite discrete form in (17) should be selected as a suitable constant when realistic application. From (17), fractional-order state feedback is essentially a fractional-order differential of the hidden states $h_p\left(k\right)$. Before discretization, the fractional-order system presented by (14) still has a Mittag-Leffler stable problem, which would influence the embedding into the network and the training process of the proposed FORNN with PIBatKnow.

3.2. Fractional-Order Constraints

As provided in the introduction Section 1, the embedded physics information in ML algorithms were mainly concluded into observational bias, inductive bias, and learning bias [24]. A typical example of learning bias is that physics are enforced via soft penalty constraints into the loss function of NNs. In this subsection, the soft penalty constraint is realized by a fractional-order PDE reflecting battery knowledge and embedded into the loss function of RNN, as shown in Figure 4. Except for the common loss of the supervised measurement of outputs, a fractional-order PDE driven by physics-informed knowledge is added as an unsupervised fractional-order loss, constructing the final loss for training backpropagation.

3.2.1. Constraint in Fractional-Order PDE Form

As shown in Figure 4, the unsupervised constraint is in a fractional-order PDE form, which can include various types of fractional-order derivatives, such as PDE loss, initial condition (IC) loss, and boundary condition (BC) loss. The fractional-order PDE constraint is firstly presented in this part, then the specific fractional-order PDE for LIB is presented in the following part. The general form of a PDE with initial condition and boundary condition is given in Definition 3.

Definition 3.

PDE with initial and boundary conditions [36].

$$\begin{array}{cc}\hfill & \mathcal{L}\left[y\right]\left(x\right)=f\left(x\right),x\in \Omega \times [0,T],\hfill \\ \hfill & \mathcal{B}\left[y\right]\left(x\right)=y_b\left(x\right),x\in \partial \Omega \times [0,T],\hfill \\ \hfill & y(x,0)=y_0\left(x\right),x\in \Omega ,\hfill \end{array}$$

(18)

where $\mathcal{L}$ and $\mathcal{B}$ represents the differential operators, Ω is the spatio-temporal domain, y is the solution, $y_b$ and $y_0$ are the boundary and initial functions.

The exact solution $y^{*}$ of the PDE in Definition 3 usually lies in an infinite dimensional space, and the FORNN with PDE constraint in this work would parameterize the solution $y^{*}$ of as y (18) to approximate the ground truth in a numerical way. If the differential operator in Definition 3 is chosen as integer-order derivative, it can obtain an integer-order PDE constraint as stated in Lemma 4.

Lemma 4.

If the differential operators in (18) is an integer-order derivative, a simple physics-informed PDE constraint can be presented as

$$f_{\mathrm{PDE}}(x,t)=\frac{\partial y}{\partial t}+a_1\frac{\partial y}{\partial x}+a_2\frac{{\partial }^{2}y}{\partial x^2},$$

(19)

where $a_1$ and $a_2$ are the coefficients of boundary condition loss (BC loss) and initial condition loss (IC loss). The PDE constraint in (19) satisfies a certain physics-informed law of the object, and (19) can be extended to a complex form with extra information by adding other types of derivatives more than IC loss or BC loss.

Remark 1.

The PDE constraint in Lemma 4 could be integer-order PDEs, integro-differential equations, fractional-order PDEs, or stochastic PDEs [24].

As stated in Remark 1, the PDE constraint in Lemma 4 may vary from different physical equations according to the investigated object in different problems, such as the viscous Burgers’ equation $\frac{\partial y}{\partial t}+y\frac{\partial y}{\partial x}=v\frac{{\partial }^{2}y}{\partial x^2}$ with a initial condition and Dirichlet boundary conditions. Considering the embedding battery knowledge represented by FOMs in this work, the integer-order PDE constraint in Lemma 4 would be extended to a fractional-order one, as stated in Lemma 5.

Lemma 5.

Combining with the various forms of fractional-order derivatives in Figure 4, and if the differential operators in (18) is a fractional-order derivative, the fractional-order PDE constraint with physics-informed knowledge can be presented as

$$f_{\mathrm{FOPDE}}(x,t)=\underset{\mathrm{PDE}\mathrm{loss}}{\underbrace{\frac{{\partial }^{\alpha }y}{\partial t^{\alpha }}}}+\underset{\mathrm{condition}\mathrm{loss}}{\underbrace{a_1\frac{{\partial }^{\alpha }y}{\partial x_1^{\alpha }}+a_2\frac{{\partial }^{\alpha }y}{\partial x_2^{\alpha }}}}+\underset{\mathrm{input}\mathrm{constraint}}{\underbrace{b_1\frac{{\partial }^{\alpha }{x}_1}{\partial t^{\alpha }}+b_2\frac{{\partial }^{\alpha }{x}_2}{\partial t^{\alpha }}}}+\underset{\mathrm{coupled}\mathrm{rate}}{\underbrace{c_1\frac{{\partial }^{\alpha }{x}_1}{\partial x_2^{\alpha }}}},$$

(20)

where $a_i,i=1,2$, $b_i,i=1,2$, and $c_1$ are the coefficients of condition loss (such as IC loss and BC loss), input constraint, and coupled rate, respectively. The fractional-order PDE constraint $f_{\mathrm{FOPDE}}(x,t)$ in (20) normalizes the coefficient of PDE loss (${\partial }^{\alpha }y/\partial t^{\alpha }$) as 1, and sets the fractional order α as the same, which can be various in application. The fractional-order PDE constraint $f_{\mathrm{FOPDE}}(x,t)$ in (20) satisfies a certain fractional-order physics-informed law of the object.

The fractional-order PDE constraint in Lemma 5 contains PDE loss, condition loss, input constraint, and coupled constraint. Among the four parts, input constraint may be the input derivative in time series, and coupled constraint may be the coupled rate or coupled relationships among inputs, which usually exists in physical systems. Note that the fractional-order PDE in (20) covers the integer-order one in (19) when the fractional-order is set as 1.

Lemma 6.

According to [24], and referring to the fractional-order PDE constraint in (20), the final loss function can be deduced as

$$\mathcal{L}=w_{\mathrm{data}}{\mathcal{L}}_{\mathrm{data}}+w_{\mathrm{PDE}}{\mathcal{L}}_{\mathrm{PDE}},$$

(21)

where

$${\mathcal{L}}_{\mathrm{data}}=\frac{1}{N_{\mathrm{data}}}{\sum }_{i=1}^{N_{\mathrm{data}}}{(y(x_i,t_i)-y_i)}^2$$

(22)

is the supervised loss of data output error, and

$$\begin{array}{c}{\mathcal{L}}_{\mathrm{PDE}}=\frac{1}{N_{\mathrm{PDE}}}{\sum }_{j=1}^{N_{\mathrm{PDE}}}{\left[f_{\mathrm{FOPDE}}(x,t)\right]}^2{|}_{(x_j,t_j)}\hfill \\ =\frac{1}{N_{\mathrm{PDE}}}{\sum }_{j=1}^{N_{\mathrm{PDE}}}{(\frac{{\partial }^{\alpha }y}{\partial t^{\alpha }}+a_i\frac{{\partial }^{\alpha }y}{\partial x_m^{\alpha }}+b_i\frac{{\partial }^{\alpha }{x}_m}{\partial t^{\alpha }}+c_i\frac{{\partial }^{\alpha }{x}_m}{\partial x_n^{\alpha }})}^2{|}_{(x_j,t_j,m=1,2,m\ne n,n=1,2)}\hfill \end{array}$$

(23)

is the unsupervised loss of fractional-order PDE constraint $f_{\mathrm{FOPDE}}(x,t)$.

In (21)–(23), $w_{\mathrm{data}}$ and $w_{\mathrm{PDE}}$ are the weights to balance the two loss ${\mathcal{L}}_{\mathrm{data}}$ and ${\mathcal{L}}_{\mathrm{PDE}}$; $y_i$ are values of y at $(x_i,t_i)$; $(x_i,t_i)$ and $(x_j,t_j)$ are two sets of points sampled at the initial/boundary locations and in the entire domain, respectively.

With Lemmas 5 and 6, the general form of a fractional-order PDE constraint is proposed to encode physics-informed knowledge into the loss function of FORNN. For a specific object, the physics-informed PDE in (20) should be extracted from certain equations by object modeling.

3.2.2. Constraint of Battery Terminal Voltage Derivative Equation

As introducing in Section 2.2, the physics-informed laws of LIBs can be described by fractional-order elements in battery FOM, which is a suitable choice to be transferred to the fractional-order PDE constraint as stated in Lemma 5. A second-order FOM is firstly considered in this work as shown in Figure 5 [37,38].

In Figure 5, the impedance of the second-order FOM in the left side can be presented as

$$Z=R_{\mathrm{ohm}}+Z_{\mathrm{ARC}}+Z_{\mathrm{Warburg}}=R_{\mathrm{ohm}}+\frac{R_{ct}}{R_{ct}{C}_1{s}^{{\alpha }_1}+1}+\frac{1}{C_2{s}^{{\alpha }_2}}$$

(24)

where E is the open circuit voltage (OCV), $U_t$ is the terminal voltage, and $R_{\mathrm{ohm}}$, RC tank $Z_{\mathrm{ARC}}$ ($R_{ct}$ and $CP{E}_1$), and $Z_{\mathrm{Warburg}}$ ($CP{E}_2$) represent the ohmic resistance in high-frequency, double layer effects and charge transfer reactions in mid-frequency, and solid diffusion dynamics in low-frequency of EIS of LIBs, respectively. Specifically, the time constant of $R_{\mathrm{ohm}}$ and RC tank $Z_{\mathrm{ARC}}$ is less than 0.05 s and $Z_{\mathrm{Warburg}}$ is longer than 50 s in the time scale [37], while the sampling frequency of general BMS is about $1\sim 10$ Hz, which cannot accurately identify the transient response of $R_{\mathrm{ohm}}$ and the dynamic response of $Z_{\mathrm{ARC}}$. Considering $R_{\mathrm{ohm}}$ and RC tank $Z_{\mathrm{ARC}}$ as one process, the second-order FOM is simplified as a RC tank as shown in the right side of Figure 5, which has the impedance as

$$\frac{U_t\left(s\right)-E\left(s\right)}{I\left(s\right)}=R_0+Z_{\mathrm{Warburg}}=R_0+\frac{1}{C{s}^{\alpha }},$$

(25)

where $R_0=R_{\mathrm{ohm}}+Z_{\mathrm{ARC}}$. Equation (25) acts as the basis of the fractional-order constraints for the loss function of PIRNN, which is furtherly deduced in the latter part.

Take the simplified FOM in Figure 5, and the corresponding voltage-current equation (25) can be presented as

$$\begin{array}{c}{U}_t\left(s\right)-E\left(s\right)=(R_0+\frac{1}{C{s}^{\alpha }})I\left(s\right)\\ \Rightarrow C{s}^{\alpha }(U_t\left(s\right)-E\left(s\right))=(R_{0}C{s}^{\alpha }+1)I\left(s\right),\end{array}$$

(26)

where $U_t$ is the terminal voltage, E is the OCV, and I is the current; $R_0$ and C are the resistance part and capacity value of the FOM in Figure 5, respectively.

Applying inverse Laplace transform, transfer (26) into time domain as

$$C·D^{\alpha }(u_t\left(t\right)-E\left(t\right))=(R_{0}C·D^{\alpha }+1)i\left(t\right),$$

(27)

where $u_t\left(t\right)$, $E\left(t\right)$, and $i\left(t\right)$ are the time-series terminal voltage, OCV, and current in time domain. Rewrite (27) as

$$C\frac{{\partial }^{\alpha }{u}_t\left(t\right)}{\partial t^{\alpha }}=C\frac{{\partial }^{\alpha }E\left(t\right)}{\partial t^{\alpha }}+R_{0}C\frac{{\partial }^{\alpha }i\left(t\right)}{\partial t^{\alpha }}+i\left(t\right),$$

(28)

In (28), the terminal voltage $u_t\left(t\right)$ and current $i\left(t\right)$ can be directly obtained from realistic sampled EV data, while OCV $E\left(t\right)$ is an inside variable of the battery that can be hardly measured in real operation conditions. However, SOC have a relationship with OCV, which can be approximately calculated as [39]

$$E\left(t\right)=f\left(SOC\right)={\sum }_{k=0}^{\infty }{d}_{k}SO{C}^k\left(t\right),$$

(29)

where $d_k$ is the coefficients of $SO{C}^k\left(t\right)$. Considering (28) as the fractional-order PDE constraint and a finite terms k in (29), some certain assumptions are proposed to replace OCV in (28) by SOC, as stated in Lemma 6.

Lemma 7.

According to the measured SOC-OCV curves in literature, let k in (29) holds $k=1$ and assume OCV $E\left(t\right)$ is monotonic with $SOC\left(t\right)$, so that the fractional-order derivative of $E\left(t\right)$ holds the relationship with the fractional-order derivative of $SOC\left(t\right)$ as

$$\frac{{\partial }^{\alpha }E\left(t\right)}{\partial t^{\alpha }}=d_{soc}\frac{{\partial }^{\alpha }SOC\left(t\right)}{\partial t^{\alpha }},$$

(30)

where $d_{soc}$ is the ratio of $({\partial }^{\alpha }E\left(t\right)/\partial t^{\alpha })$ to $({\partial }^{\alpha }SOC\left(t\right)/\partial t^{\alpha })$.

Take (30) into (28), we obtain the fractional-order PDE constraint from the FOM of battery as

$$C\frac{{\partial }^{\alpha }{u}_t\left(t\right)}{\partial t^{\alpha }}=d_{soc}C\frac{{\partial }^{\alpha }SOC\left(t\right)}{\partial t^{\alpha }}+R_{0}C\frac{{\partial }^{\alpha }i\left(t\right)}{\partial t^{\alpha }}+i\left(t\right).$$

(31)

For FORNN with PIBatKnow, the inputs and outputs should be carefully selected and combined with the fractional-order constraint (31). Moreover, it should consider the available variables from realistic data, the influence of temperature, and the aim of the proposed FORNN. For example, taking state estimation as target, then the inputs are defined as $\mathbf{x}\left(t\right)={[x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)]}^{\mathrm{T}}={[u_t\left(t\right),i\left(t\right),T_{temp}\left(t\right)]}^{\mathrm{T}}$, and the output as $y\left(t\right)=SOC\left(t\right)$, thus the fractional-order PDE constraint can be deduced as

$$f_{SOC}(x,t)=d_{soc}C\frac{{\partial }^{\alpha }y\left(t\right)}{\partial t^{\alpha }}+R_{0}C\frac{{\partial }^{\alpha }{x}_2\left(t\right)}{\partial t^{\alpha }}+x_2\left(t\right)-C\frac{{\partial }^{\alpha }{x}_1\left(t\right)}{\partial t^{\alpha }},$$

(32)

Combining (32) with (23), the loss of fractional-order PDE constraint can be calculated in the form of a mean sum of square (MSS), and embedded into the loss function as a fractional-order unsupervised loss. G-L definition can be employed to discrete the fractional-order derivatives in (32) for calculation. Select finite terms $L=5$ in (2), and the fractional-order derivatives of inputs $x\left(t\right)$ and output $y\left(t\right)$ in (32) can be deduced as

$${\Delta }^{\alpha }M\left(t\right)=\sum _{j=0}^5{(-1)}^j\left(\begin{array}{c}\alpha \hfill \\ j\hfill \end{array}\right)M(k-j),$$

(33)

where $M\left(t\right)=x_1\left(t\right),x_2\left(t\right),\mathrm{or}y\left(t\right)$, and k is the discrete step. Note that the discrete equation in (33) requires the history data in previous $\{5-j,(j=0,1,\dots ,5\left)\right\}$ moments, which can be found in the time-series data of the realistic sampled voltage, current, and the output SOC. With (32) and (33), a FORNN with a fractional-order PDE constraint is constructed and encoded by battery voltage equation, then the proposed FORNN can be applied to conduct SOC estimation or further life prediction.

3.3. Fractional-Order Descent Methods

Besides fractional-order state feedback and fractional-order constraint, the FORNN can also be enhanced in the backpropagation process by fractional-order gradients. In the opposite direction of forward propagation, backpropagation starts from the output layer to the input layer, mainly updating weights $W_p$, $W_{hp}$, and bias $b_p$, $(p=1,2,\dots ,l)$ in the gradient descent direction, which is called the gradient descent (GD) method, as shown in Figure 6a. It always makes the weights $W_p=W_{hp}$ in a realistic application, and the bias $b_p$ holds the same update with weight $W_p$, thus only weight $W_p$ is discussed in the following. During backpropagation, the layer state gradients are calculated, then weight gradients are based on the layer gradients to update weight $W_p$. The gradients of the $p\mathrm{th}$ hidden layer for the $k\mathrm{th},(k=1,2,\dots ,N)$ input sample are presented as

$$\left\{\begin{array}{cc}\hfill \frac{\partial L\left(g\right(x),y)}{\partial h_p^k\left(x\right)}=& W_{p+1}^k\frac{\partial L\left(g\right(x),y)}{\partial a_{p+1}^k\left(x\right)}\hfill \\ \hfill \stackrel{}{\frac{\partial L\left(g\right(x),y)}{\partial a_p^k\left(x\right)}}=& \frac{\partial L\left(g\right(x),y)}{\partial h_p^k\left(x\right)}·g^{\prime }\left(a_p^k\left(x\right)\right)=W_{p+1}^k\frac{\partial L\left(g\right(x),y)}{\partial a_{p+1}^k\left(x\right)}·g^{\prime }\left(a_p^k\left(x\right)\right),\hfill \end{array}\right.$$

(34)

where $\partial L\left(g\right(x),y)/\partial h_p^k\left(x\right)$ and $\partial L\left(g\right(x),y)/\partial a_p^k\left(x\right)$ are the gradients of the output $h_p^k\left(x\right)$ and the input $a_p^k\left(x\right)$ of the $p\mathrm{th}$ hidden layer, respectively. Since the backpropagation starts from the output layer, it supposes that the gradients $\partial L\left(g\right(x),y)/\partial h_t^k\left(x\right)$ and $\partial L\left(g\right(x),y)/\partial a_t^k\left(x\right)$ of the $t\mathrm{th},(t=p+1,p+2,\dots ,l,l+1)$ layer are already known. With layer gradients in (34), fractional-order derivative is introduced from the output layer to calculate the fractional-order gradients of loss function $L\left(g\right(x),y)$; then, all the weights gradients of hidden layers are transferred to fractional-order ones, which turns out to be a fractional-order gradient descent (FOGD) method, as shown in Figure 6b. For time-series battery dataset as the inputs of the proposed FORNN, if considering the gradient descent method with momentum (GDm) shown in Figure 6c, the FOGD method can be extended to the FOGD method with momentum (FOGDm), as presented in Figure 6d.

3.3.1. Fractional-Order Gradient Descent

FOGD employs the fractional-order gradients of the loss function and weights $W_p$ instead of the integer-order gradients. Based on the fractional-order gradient in (9), the FOGD for the updates of weight $W_p$ is presented as

$$W_p^{k+1}=W_p^k-\eta {·}_{W_p^0}{D}_{W_p^k}^{\alpha }L(f\left(x\right),y)=W_p^k-\eta ·\frac{{\partial }^{\alpha }L(f\left(x\right),y)}{\partial {\left(W_p^k\right)}^{\alpha }}$$

(35)

According to $a_p\left(x\right)=W_p{h}_{p-1}\left(x\right)+b_p$ in (12) and the chain rule, (35) is deduced as

$$W_p^{k+1}=W_p^k-\eta ·\frac{\partial L\left(g\right(x),y)}{\partial a_p^k\left(x\right)}·\frac{{\partial }^{\alpha }{a}_p^k\left(x\right)}{\partial {\left(W_p^k\right)}^{\alpha }}$$

(36)

where $\partial L\left(g\right(x),y)/\partial {\left(W_p^k\right)}^{\alpha }$ is the fractional-order gradients of weight $W_p$ to the loss function $L\left(g\right(x),y)$, and $\eta$ is the learning rate (iteration step size). It should be noted that, according to Lemma 3 in Section 2.1.1, the chain rule in (36) for fractional-order partial derivatives only validates when the loss function $L\left(g\right(x),y)$ holds an exponential form as the function $f\left(x\right)$ in (5). In (36), the gradient of the input $a_p^k\left(x\right)$ to the loss function ($\partial L\left(g\right(x),y)/\partial a_p^k\left(x\right)$) can be obtained by (34), and ${\partial }^{\alpha }{a}_p^k\left(x\right)/\partial {\left(W_p^k\right)}^{\alpha }$ can be calculated by (5). Hence, combining (34) and the discrete Caputo definition in (5), the fractional-order gradients of weight $W_p$ to the loss function $L\left(g\right(x),y)$ can be deduced as

$$\frac{{\partial }^{\alpha }L(g\left(x\right),y)}{\partial {\left(W_p^k\right)}^{\alpha }}=\frac{h_{p-1}^k\left(x\right)}{\Gamma (2-\alpha )}·\frac{\partial L\left(g\right(x),y)}{\partial a_p^k\left(x\right)}{\left(W_p^k-W_p^0\right)}^{1-\alpha }$$

(37)

where $W_p^0$ is the initial values of weight $W_p$. Take (37) into (36), we can obtain the FOGD method for backpropagation as

$$W_p^{k+1}=W_p^k-\eta ·\frac{h_{p-1}^k\left(x\right)}{\Gamma (2-\alpha )}·\frac{\partial L\left(f\right(x),y)}{\partial a_p^k\left(x\right)}{\left(W_p^k-W_p^0\right)}^{1-\alpha }$$

(38)

where $\alpha$ is a fractional order that may be related to battery knowledge and sensitive to the training results. The sensitivity of the fractional order $\alpha$ and the corresponding influence on network output is discussed in the experimental part.

3.3.2. Fractional-Order Gradient Descent with Momentum

FOGDm method is extended from FOGD method by adding momentum, which demonstrates the changing direction of the fractional-order gradients. Given the integer-order gradient descent (GD) method with momentum as

$$\left\{\begin{array}{cc}\hfill V_p^{k+1}& =\mu ·V_p^k-\eta ·\frac{\partial L\left(g\right(x),y)}{\partial W_p^k}\hfill \\ \hfill W_p^{k+1}& =W_p^k+V_p^{k+1}\hfill \end{array}\right.$$

(39)

where $\mu$ is the momentum factor, which is the proportion of history gradient data, and $V_p$ is the weight momentum. Then the weight update by FOGDm method can be presented as

$$\left\{\begin{array}{cc}\hfill V_p^{k+1}& =\mu ·V_p^k-\eta ·\frac{{\partial }^{\alpha }L(g\left(x\right),y)}{\partial {\left(W_p^k\right)}^{\alpha }}\hfill \\ \hfill W_p^{k+1}& =W_p^k+V_p^{k+1}\hfill \end{array}\right.$$

(40)

Take the fractional-order gradients of weight $W_p$ to the loss function $L\left(g\right(x),y)$ in (37) into (40), we obtain the discrete updating equation of weight $W_p$ by the FOGDm method as

$$\left\{\begin{array}{cc}\hfill V_p^{k+1}& =\mu ·V_p^k-\eta ·\frac{h_{p-1}^k\left(x\right)}{\Gamma (2-\alpha )}·\frac{\partial L\left(g\right(x),y)}{\partial a_p^k\left(x\right)}{\left(W_p^k-W_p^0\right)}^{1-\alpha }\hfill \\ \hfill W_p^{k+1}& =W_p^k+V_p^{k+1}\hfill \end{array}\right.$$

(41)

With (41) to encode the momentum of fractional-order gradient into the backpropagation process, the FOGDm method accelerates the convergence of the proposed FORNN with faster updates of weight $W_p$ and bias $b_p$, and adjusts fractional-order gradients more rapidly than the FOGD method.

4. Experiment Setups

4.1. Battery and Operation Conditions

The proposed FORNN with PIBatKnow in this work may contain fractional-order state feedback, fractional-order loss constraint, and fractional-order gradient for backpropagation, which bring more parameters for algorithm tuning. Hence, the sensitivity analysis is necessary to identify suitable values and maximize the effectiveness for the practical estimation of LIB. For sensitivity analysis, experiments under five temperatures (5 ${}^{\circ }\mathrm{C}$, 15 ${}^{\circ }\mathrm{C}$, 25 ${}^{\circ }\mathrm{C}$, 35 ${}^{\circ }\mathrm{C}$, 45 ${}^{\circ }\mathrm{C}$) are conducted for a 18650 cell to collect battery data, and FUDS operation condition is employed to simulate the dynamic working condition. The parameters of the 18650 cell is shown in

Table 1. LIB 18650 cell parameters.

Parameters	Values	Parameters	Values
Rated capacity (0.5${\mathrm{C}}_5$A)	2000 mAh	Rated voltage	3.7 V
Max charge voltage	4.2 V	Discharge cut-off voltage	2.75 V
Working temperature (charge)	0 ${}^{\circ }\mathrm{C}$–45 ${}^{\circ }\mathrm{C}$	Working temperature (discharge)	−20 ${}^{°}\mathrm{C}$–60 ${}^{°}\mathrm{C}$

. The experimental data is sampled by a BTS-4 battery tester, and the cell is put in an incubator to control the testing temperature. The cell is fully charged by constant-current constant-voltage (CC-CV) method before applying FUDS, then fully discharged to the cut-off voltage 2.75 V under cycling FUDS periods.

4.2. Dataset and Initialization

With the collected data, SOC is set as the estimation target of the proposed FORNN with PIBatKnow, so that the relationship between the input and the output can be constructed and analyzed clearly. The real SOC values are calculated by the ampere-hour integral method. Take the influence of temperature into consideration, this work for the sensitivity analysis takes the current, voltage, and temperature as the three main dimensions related to SOC. Hence, the inputs are selected as $\mathbf{x}\left(t\right)={[x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)]}^{\mathrm{T}}={[u_t\left(t\right),i\left(t\right),T_{temp}\left(t\right)]}^{\mathrm{T}}$, and the output as $y\left(t\right)=SOC\left(t\right)$. Figure 7 is the current and voltage data after preprocess. The total amount of data is 27,213 points, which can support the big-data learning during training process for the proposed algorithm. The collected data is divided according to the five temperatures, that is, 45 ${}^{\circ }\mathrm{C}$ is selected as the testing data, and data under 5 ${}^{\circ }\mathrm{C}$, 15 ${}^{\circ }\mathrm{C}$, 25 ${}^{\circ }\mathrm{C}$, and 35 ${}^{\circ }\mathrm{C}$ is divided into training data and validation data. The division ratio is 0.75:0.048:0.202, resulting in training data (20,410 points), validation data (1306 points), and testing data (5497 points). Since the aim of this work is the sensitivity analysis of parameters, only the data of one cell is enough, and any one of five temperatures (not specific 45 ${}^{\circ }\mathrm{C}$) can be selected as testing data. Moreover, some unchanged parameters are initialized and listed in

Table 2. Unchanged parameters of the proposed FORNN with PIBatKnow.

Name	Value/Range	Name	Value/Range
Hidden layers	1	Hidden neurons	12
max epoch $E_{\max }$	300	Performance function	MSE
train:valid:test	0.75:0.048:0.202	Training goal	1.6 $\times {10}^{-4}$

.

5. Sensitivity Analysis and Estimation Results

After the dataset division and initialization of the proposed algorithm, the sensitivity of the main parameters related to the fractional-order knowledge are investigated in this work. According to Section 3, nine fractional-order parameters are investigated and divided into three categories, as shown in

Table 3. Sensitivity categories of the fractional-order constraint and the fractional-order gradient in the proposed algorithm.

Type	Name	Value/Range	Attribution
Fractional-order Gradient sensitivity	Fractional order ${\alpha }_1$	$[0.1,1]$	FOGDm in (41)
	Momentum weight $\mu$	$[0.1,1]$
	Learning rate $\eta$	$[0.08,0.26]$
Impedance sensitivity	Fractional order ${\alpha }_2$	$[0.1,1]$	FO PDE $f_{SOC}(x,t)$ in (32)
	Ratio of OCV-SOC $d_{soc}$	$[5,50]$
	Capacitance $C_{\mathrm{bat}}$ (unit: C)	$[2.5,25]$
	ohm resistance $R_0$ (unit: $\Omega$)	[5 $\times {10}^{-3}$, 6 $\times {10}^{-3}$]
Loss weight sensitivity	Loss weight $w_{\mathrm{data}}$	$[0.1,1]$	final loss $\mathcal{L}$ in (20)
	Loss weight $w_{\mathrm{PDE}}$	$1-w_{\mathrm{data}}$

, that is, parameters’ sensitivity in the FOGD and FOGDm method (fractional-order gradient sensitivity), parameters’ sensitivity in fractional-order PDE constraint (impedance sensitivity), and weights in loss calculation (loss weight sensitivity). A certain possible range of the parameters are also provided in Table 3, in which only capacitance $C_{\mathrm{bat}}$ and ohm resistance $R_0$ have physic units, that is, C and $\Omega$, respectively, and other seven fractional-order parameters hold the unit 1. For example, the parameter range in impedance sensitivity is determined according to the battery FOM in (25) and Figure 5. With the three categories in Table 3 and other initialized network parameters in Table 2, the sensitivity analysis is conducted in the following sections, respectively.

5.1. Estimation with Fractional-Order Gradient Sensitivity

The first sensitivity category is the fractional-order gradient sensitivity, which contains the fractional order ${\alpha }_1$, the momentum weight $\mu$, and the learning rate $\eta$ in the FOGD and FOGDm method. Fractional order reflects and embeds battery fractional-order characteristics. In this section, only fractional order ${\alpha }_1$, momentum weight $\mu$, and learning rate $\eta$ change in the range provided in Table 3, respectively. Other parameters stay unchanged as shown in

Table 4. Default values of the nine main fractional-order parameters.

${\mathit{\alpha }}_1$	$\mathit{\mu }$	$\mathit{\eta }$	${\mathit{\alpha }}_2$	${\mathit{d}}_{\mathit{soc}}$	${\mathit{C}}_{\mathbf{bat}}$	${\mathit{R}}_0$	${\mathit{w}}_{\mathbf{data}}$	${\mathit{w}}_{\mathbf{PDE}}$
0.9	0.75	0.18	0.9	40	20	0.005	0.8	0.2

, and the default values of fractional order ${\alpha }_1$, momentum weight $\mu$, and learning rate $\eta$ are also provided when acting as the unchanged values in the other two sensitivity categories.

To verify the estimation effects, SOC is taken as the output of the proposed FORNN with PIBatKnow. Figure 8 presents the SOC estimation results under the sensitivity of the fractional order ${\alpha }_1$ (${\alpha }_1$ = 0.1:0.1:1) in FOGD and FOGDm methods. Expression variable = a:interval:b denotes that $variable$ changes from a to b in the range $[a,b]$ and takes the values in every interval $interval$, then this expression is applied to all the parameter sensitivity in the following part. Combined with Section 4.2, Figure 8b contains five discharging snippets under FUDS condition in five temperature (5 ${}^{\circ }\mathrm{C}$, 15 ${}^{\circ }\mathrm{C}$, 25 ${}^{\circ }\mathrm{C}$, 35 ${}^{\circ }\mathrm{C}$, 45 ${}^{\circ }\mathrm{C}$). Snippet in 45 ${}^{\circ }\mathrm{C}$ is the testing dataset, whose outputs and errors are enlarged in Figure 8c and Figure 8e, respectively.

The testing loss in Figure 8a shows that loss decreases and converges faster with larger ${\alpha }_1$, but the training process turns unstable when $\alpha >0.5$, especially $\alpha =1$. Moreover, the outputs of all dataset (training, validation, and testing) and errors (errors = target − outputs) are presented in Figure 8b and Figure 8d, respectively. Figure 8b,d also show the unstable accuracy performance. For fractional order ${\alpha }_1$, there exists a trade-off between performance and stability. With the outputs and the corresponding errors, it proves that the training dataset and testing dataset have similar performance, and the network does not overfit or underfit.

Similarly, the sensitivity of momentum $\mu$ and learning rate $\eta$ are provided in Figure 9 and Figure 10, respectively. Since the testing prediction is the main point, and performance of the training or validation datasets are consistent with the testing dataset, only the testing outputs and testing errors are provided to focus on the testing performance, which is different from Figure 8.

From Figure 9, larger momentum weight $\mu$ would cause lower convergence speed and more instability, especially when $\mu =1.0$, in which the network cannot converge, but the accuracy is enhanced by a larger $\mu$. It should be noted that only the FOGDm method has the momentum weight $\mu$, which determines the ratio of the previous momentum $V_p$, as shown in (40). As to the learning rate $\eta$ in Figure 10, it is very sensitive to the convergence speed and algorithm stability but not so sensitive to the accuracy, because the testing loss in Figure 10a varies a lot with different $\eta$, while the outputs and errors in Figure 10b,c do not show too much difference. A value of $\eta$ larger than 0.2 is not suitable for the proposed FORNN with PIBatKnow.

5.2. Estimation with Impedance Sensitivity

The second sensitivity category is the impedance sensitivity. For the fractional-order PDE constraint in loss function of the proposed FORNN with PIBatKnow, the main parameters include ${\alpha }_2$, $d_{soc}$, $C_{\mathrm{bat}}$, and $R_0$ of the battery FOM, and the sensitivity of the four impedance parameters are investigated in this section.

Figure 11, Figure 12, Figure 13 and Figure 14 are the testing loss, testing outputs, and corresponding errors under the four impedance parameters ${\alpha }_2$, $d_{soc}$, $C_{\mathrm{bat}}$, and $R_0$, respectively. To better illustrate the loss value, the testing loss is enlarged within the early stage of epochs, and a subfigure enlarging the stable loss in later epochs is added in every figure of testing loss, like Figure 11a. Part of the training process would be terminated before the max epoch (300) when reaching the training goal (1.6 $\times {10}^{-4}$) or other validation conditions. For example, in the subfigure of Figure 11a, the training process with ${\alpha }_2=0.1$ stopped at 159th epoch. From the accuracy aspect, the outputs and errors in Figure 11, Figure 12, Figure 13 and Figure 14 do not show large fluctuations when the impedance parameters change. Most of the errors with various ${\alpha }_2$, $d_{soc}$, $C_{\mathrm{bat}}$, and $R_0$ maintain the magnitude within $0.15$.

From Figure 11, the testing loss decreases and becomes more stable with larger fractional order ${\alpha }_2$, and the errors also decrease to within $0.05$, as shown in Figure 11c. From Figure 12, the testing loss and outputs accuracy do not show certain sensitive pattern when the ratio $d_{soc}$ between OCV and SOC changes, but the testing loss becomes more unstable with larger $d_{soc}$. From Figure 13, the testing loss converges faster when the capacitance $C_{\mathrm{bat}}$ is in the middle value (such as $15,17.5,20$). The stable loss becomes flatter and the errors become smaller with larger $C_{\mathrm{bat}}$. From Figure 14, a large resistance $R_0$ would cause unstable convergence process and large testing loss, but the estimation accuracy in high SOC range may be improved.

5.3. Estimation with Loss Weight Sensitivity

The third sensitivity category is the loss weight sensitivity. The loss weight $w_{\mathrm{data}}$ and $w_{\mathrm{PDE}}$ can determine the ratio of the supervised measurement loss and the unsupervised fractional-order loss calculated by physics-informed PDE in Figure 4. Since loss weight $w_{\mathrm{data}}$ and $w_{\mathrm{PDE}}$ hold the relationship $w_{\mathrm{PDE}}=1-w_{\mathrm{data}}$, the two parameters may be considered as one parameter sensitivity, as presented in Figure 15. $w_{\mathrm{data}}$ = 0.1:0.1:1 corresponds to $w_{\mathrm{PDE}}$ = 0.9:0.1:0. The testing loss in Figure 15a demonstrates that larger $w_{\mathrm{data}}$ makes the network process converge faster and more stable, which reflects the relative instability of the physics-informed knowledge embedding. However, the algorithm can achieve a trade-off with sensitivity analysis and suitable tuning. As to the estimation accuracy, the outputs and errors in Figure 15b,c do not show much improvement both in the small value and large value of $w_{\mathrm{data}}$.

5.4. Correlation Analysis

Correlation analysis is provided in this section to present the dynamic relationship between parameters and the proposed algorithm performance. Based on the possible range in Table 3, the three categories (nine parameters) are sampled as ten groups of values randomly, as shown in

Table 5. Ten groups of radom values of the nine parameters in the three categories for correlation analysis.

No.	${\mathit{\alpha }}_1$	$\mathit{\mu }$	$\mathit{\eta }$	${\mathit{\alpha }}_2$	${\mathit{d}}_{\mathit{soc}}$	${\mathit{C}}_{\mathbf{bat}}$	${\mathit{R}}_0$	${\mathit{w}}_{\mathbf{data}}$	${\mathit{w}}_{\mathbf{PDE}}$	$\mathit{perf}$	$\mathit{trperf}$	$\mathit{vperf}$	$\mathit{tperf}$
1	1	0.72	0.125	0.5	11.75	22.75	5.24 $\times {10}^{-3}$	0.89	0.11	0.001349	0.001158	0.002434	0.001799
2	0.47	0.6	0.251	0.23	39.2	18.025	5.34 $\times {10}^{-3}$	0.66	0.34	0.001327	0.001177	0.001309	0.001889
3	0.54	0.98	0.152	0.16	23.45	24.325	5.45 $\times {10}^{-3}$	0.38	0.62	0.218343	0.232547	0.005316	0.216220
4	0.08	0.66	0.1484	0.32	33.35	11.725	5.12 $\times {10}^{-3}$	0.3	0.7	0.002358	0.002098	0.003900	0.002954
5	0.97	0.15	0.2186	0.6	8.15	23.425	5.52 $\times {10}^{-3}$	0.29	0.71	0.029044	0.032181	0.003943	0.023360
6	0.62	0.13	0.1088	0.47	36.95	9.925	5.85 $\times {10}^{-3}$	0.89	0.11	0.010673	0.009013	0.011097	0.016737
7	0.45	0.91	0.233	0.78	45.5	3.625	5.67 $\times {10}^{-3}$	0.15	0.85	0.004387	0.004279	0.000754	0.005652
8	0.91	0.92	0.1394	0.82	16.25	16.675	5.93 $\times {10}^{-3}$	0.91	0.09	0.000861	0.000894	0.000692	0.000782
9	0.81	0.14	0.1934	0.28	45.5	6.775	5.45 $\times {10}^{-3}$	0.52	0.48	0.024046	0.025149	0.003720	0.024783
10	0.12	0.04	0.0854	0.66	12.2	2.95	5.19 $\times {10}^{-3}$	0.34	0.66	1.044450	1.235187	0.046351	0.573336

. Mean square error (MSE) is selected as the evaluation function in this work, thus the change of MSE is considered as the indicator of the algorithm performance, marked as $perf$ and presented in Table 5. Despite the total MSE $perf$, the performance of training dataset, validation dataset, and testing dataset can be calculated as $trperf$, $vperf$, and $tperf$, respectively.

To quantize the dynamic trend between the nine fractional-order parameters and algorithm performance, the Pearson correlation coefficients are calculated with Table 3 and listed in

Table 6. Correlation coefficients of performance with nine main fractional-order parameters.

MSE	${\mathit{\alpha }}_1$	$\mathit{\mu }$	$\mathit{\eta }$	${\mathit{\alpha }}_2$	${\mathit{d}}_{\mathit{soc}}$	${\mathit{C}}_{\mathbf{bat}}$	${\mathit{R}}_0$	${\mathit{w}}_{\mathbf{data}}$	${\mathit{w}}_{\mathbf{PDE}}$
$perf$	0.5103	0.3907	0.5173	0.1686	0.3843	0.3870	0.3737	0.2855	0.2855
$trperf$	0.5104	0.3977	0.5161	0.1794	0.3839	0.3957	0.3744	0.2821	0.2821
$vperf$	0.5108	0.5532	0.5948	0.1963	0.3357	0.4748	0.3196	0.1887	0.1887
$tperf$	0.5019	0.3272	0.5174	0.0812	0.3818	0.3115	0.3632	0.3091	0.3091

The numbers are in bold to highlight the strong correlation coefficients larger than 0.5 or near 0.5 (smaller than 0.5 but within 0.05).

. From the correlation coefficients in Table 6, the fractional order ${\alpha }_1$ in gradients and learning rate $\eta$ have a strong influence on the algorithm, while the momentum weight $\mu$ in the gradient method mainly contributes to the validation performance. The impedance (mainly $d_{soc}$, $C_{\mathrm{bat}}$, and $R_0$) of battery FOM (25) have strong correlations with the validation performance and testing performance, while the fractional order ${\alpha }_2$ has weaker influence on the testing dataset. The loss weight $w_{\mathrm{data}}$ and $w_{\mathrm{PDE}}$ have weaker influence on training performance than the other parameters, but the two of them have the same correlation coefficients due to the relationship $w_{\mathrm{PDE}}=1-w_{\mathrm{data}}$. The total performance $perf$ and training performance $trperf$ have very close values because the division ratio is 0.75:0.048:0.202, which demonstrates that training performance constructs the main part of the total performance.

6. Discussion

Combined with the estimation results under the three sensitivity categories and the correlation analysis, the sensitivity and correlation of the nine main fractional-order parameters to the algorithm performance can be concluded as shown in

Table 7. Positive and negative correlation of the nine main fractional-order parameters to the algorithm performance.

index	${\mathit{\alpha }}_1$	$\mathit{\mu }$	$\mathit{\eta }$	${\mathit{\alpha }}_2$	${\mathit{d}}_{\mathit{s}\mathit{o}\mathit{c}}$	${\mathit{C}}_{\mathbf{bat}}$	${\mathit{R}}_0$	${\mathit{w}}_{\mathbf{data}}$	${\mathit{w}}_{\mathbf{PDE}}$
speed 1	↗	↘	↘	↘	↗	middle 4	↘	↗	↘
loss 2	↗	↗	-	↗	-	-	↘	-	-
accuracy 3	↗	↗	-	↗	-	↗	↗	-	-
stability	↘	↘	↘	↗	↘	↗	↘	↗	↘

¹ speed means convergence speed. ² Eloss means testing loss. ³ accuracy means estimation accuracy. ⁴ middle means that the algorithm achieves faster convergence speed in the middle values of C_bat. The direction of arrows means postive (upward arrow) or negative (downward arrow) correlation trends.

.

Take convergence speed, testing loss, estimation accuracy, and stability as the algorithm performance indexes, several notes can be listed as references and instructions for the tuning of the nine main fractional-order parameters in the proposed FORNN with PIBatKnow as follows.

• For the convergence speed, it would be boosted by a larger value of ${\alpha }_1$, $d_{soc}$, and $w_{\mathrm{data}}$, and by a smaller value of $\mu$, $\eta$, ${\alpha }_2$, $R_0$, and $w_{\mathrm{PDE}}$;
• For the testing loss, it would be improved by a larger value of ${\alpha }_1$, $\mu$, and ${\alpha }_2$, and by a smaller value of $R_0$, but does not show sensitivity to $\eta$, $d_{soc}$, $C_{\mathrm{bat}}$, $w_{\mathrm{PDE}}$, and $w_{\mathrm{data}}$;
• For the estimation accuracy, it would be improved by a larger value of ${\alpha }_1$, $\mu$, ${\alpha }_2$, $C_{\mathrm{bat}}$, and $R_0$, but does not show sensitivity to $\eta$, $d_{soc}$, $w_{\mathrm{PDE}}$, and $w_{\mathrm{data}}$;
• For the algorithm stability, it would be enhanced by a larger value of ${\alpha }_2$, $C_{\mathrm{bat}}$, and $w_{\mathrm{PDE}}$, and by a smaller value of ${\alpha }_1$, $\mu$, $\eta$, $d_{soc}$, $R_0$ and $w_{\mathrm{PDE}}$;
• According to the correlation analysis, if the nine fractional-order parameters can be tuned adaptively, the fractional order ${\alpha }_1$, the learning rate $\eta$, and the capacitance $C_{\mathrm{bat}}$ have the most dynamic correlation to the performance;
• For the fractional order ${\alpha }_1$ in FOGD and FOGDm methods, a value in range $(0.5,1)$ is suitable and a trade-off would be made between performance and stability;
• The ratio of the previous momentum $V_p$, larger $\mu$ means larger inertia of the previous convergence, which makes the speed slow but improves the learning ability to achieve higher accuracy;
• The proposed algorithm achieves faster convergence speed in the middle values of the capacitance $C_{\mathrm{bat}}$ in battery FOM;
• The loss weight $w_{\mathrm{PDE}}$ has opposite tuning direction to the loss weight $w_{\mathrm{data}}$.

7. Conclusions

This paper presents the sensitivity analysis of a fractional-order recurrent neural network (FORNN) with physics-informed knowledge of lithium-ion battery. A recurrent neural network is physics-informed by encoding with fractional-order battery knowledge, resulting in the proposed algorithm, called FORNN with PIBatKnow. Three types of physics-informed patterns are proposed as the fractional-order state feedback for forward propagation, the fractional-order PDE constraint for loss function, and the FOGD methods for backpropagation. Then, the sensitivity of nine main fractional-order parameters from the FOGD and FOGDm methods, the battery FOM, and the loss calculation are taken into consideration. Both SOC estimation experiment under parameter sensitivity and correlation analysis are conducted. Several notes and discussions are presented to provide possible instructions for the tuning of the fractional-order parameters in the proposed FORNN with PIBatKnow. The discretization of fractional-order derivatives is the key point to embed fractional-order dynamics into neural network. For example, the fractional-order state feedback can be achieved by changing the discrete state calculating equations of the hidden layer; the fractional-order constraints need to rewrite the loss function code with the GL discrete form; the FOGD methods need to rewrite the training function code to introduce fractional-order gradients in the discrete form. As presented in the theoretical section Section 3, the fractional-order state feedback, the fractional-order constraints, and the FOGD methods can also be applied to other objects more than batteries. It only needs to transfer the fractional-order constraints related to the fractional-order physics knowledge of the other objects, and all methods in this work can be easily employed.

In future research, an enhanced FORNN can be furtherly investigated by improving fractional-order dynamics of the RNN system. The Caputo and GL derivatives used in this work can partly capture the battery mechanism dynamics to achieve faster convergence speed or higher estimation accuracy, and further discussion about the employed fractional-order definitions in the FORNN can be conducted in our future work. Specifically, the fractional-order derivatives can be furtherly improved, such as employing the fractional-order Caputo-Fabrizio or Atangana-Baleanu-Caputo derivative instead of the Caputo and GL derivative, the explanations of fractional-order derivative in the FORNN, and so on. Moreover, the fractional-order neural network is the combination of battery fractional-order modeling with machine learning; thus, more fractional-order information may be added into the network design to develop a physics-informed system. Hence, the fractional-order neural network can be enhanced by designing a specialized architecture of RNN. For example, FORNN with physics embedded in the loss function can be extended into multi-task learning algorithms, by introducing a set of physical constraints from various LIBs.