**1. Introduction**

Due to the high operating voltage and high energy density of lithium-ion batteries, both grid and off-grid applications using lithium-ion batteries have gained a lot of attention [1], especially in the field of electric vehicles [2]. In addition to the advancement of battery manufacturing technology, a sufficiently accurate and agile battery management system (BMS) is important for the further application of lithium-ion batteries [3].

During battery charge/discharge cycles, the external characteristics of batteries behave as a time-varying nonlinear system due to the nonlinear relationship between electric potential versus State of Charge (*SOC*), as well as the voltage drop due to various polarization phenomena (e.g., ohmic polarization, activation polarization, concentration polarization) [4]. It is important to clarify that polarization is a universal phenomenon that occurs inside the lithium-ion battery during operation, with different types of polarization occurring at different locations inside the battery. Furthermore, the relative proportion of each type of polarization varies with the change of external excitation [5]. Polarization hinders lithium intercalation and deintercalation kinetics, leading to a decline in energy efficiency and performance of the battery [6]. Therefore, polarization is of concern during the entire phase from battery material development to the end-of-life phase of batteries. In addition, higher charge/discharge rates, extremely low ambient temperatures, and increased cycles, which are common scenarios encountered during the use of power batteries, can all lead to increased battery polarization [7].

Polarization is concerned in different stages of battery development and application.

**Citation:** Xia, B.; Ye, B.; Cao, J. Polarization Voltage Characterization of Lithium-Ion Batteries Based on a Lumped Diffusion Model and Joint Parameter Estimation Algorithm. *Energies* **2022**, *15*, 1150. https:// doi.org/10.3390/en15031150

Academic Editor: Alon Kuperman

Received: 12 January 2022 Accepted: 31 January 2022 Published: 4 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).


charging process to make a balance between charging time and temperature rise and combined it with a Genetic Algorithm (GA) to find the optimal charging current trajectory. Verification experiments of battery aging proved that this charging method has similar capacity retention to the 0.5C Constant current Constant voltage (CC-CV) charging method, but with improved charging efficiency [22,23].

The main characterization methods of battery polarization can be roughly divided into experiment-based methods and model-based methods. For model-based methods, the electrochemical model and equivalent circuit model are two of the most common models adopted.

The experiment-based method probes the polarization characteristics of a battery cell using invasive or noninvasive methods through physical and chemical probing techniques. Ex-situ X-ray diffraction technology is adopted to observe polarization behavior during the initial lithiation process in the oxide of WO3 and derivatives, and polarization arises from conversion reaction alleviated using the WO3−*<sup>x</sup>* oxide compared to the WO3 oxide, due to higher electrical conductivity [9]. Xu et al. used synchrotron X-ray tomography analysis and microstructure-resolved computational modeling to analyze the morphological defects of electrodes from multiple spatial scales, combined with battery frequency domain analysis, to explore the correlation between polarization and morphological defects of electrode structures [24]. The GITT method uses the subtraction of the quasi-open-circuit voltage (QOCV) from the closed-circuit voltage (CCV) as a characterization of electrode polarization [9,12,25]. Noelle D.J. et al. imposed abuse conditions on a battery and direct current (DC) internal resistance analysis was used to quantitatively characterize the ohmic and polarization resistance during thermal runaway [13]. Furthermore, a destructive intrusion test was conducted to investigate the relationship between electrolyte concentration and the polarization internal resistance of the battery in Noelle's research.

An electrochemical model is a common tool in battery modelling and simulation [26]. Nyman A. et al. adopted the Pseudo-2D model to locate and quantify the various types of polarization occurring inside the cell, but in this study only the polarization occurrence at specific *SOC* levels (40% and 80%) was investigated [5]. Huang et al. developed a coupled electrochemical-thermal model based on a one-dimensional electrochemical model with COMSOL Multiphysics software to study the effect of discharge rate on heat production of the cell, including the heat of polarization [14]. Li et al. improved the traditional Single Particle Model (SPM) and identified the model parameters by excitation response analysis and conducted experiments to confirm the model's performance at large charge/discharge C-rate conditions (up to 4C discharge condition for LiCoO2 batteries, and up to 5C discharge for LiFePO4 batteries) [3]. Yan et al. developed a 3D model that preserves the effect of inhomogeneous geometrical characteristics on global polarization as well as the local polarization in the electrodes, which can provide more information on the polarization characteristics of the real cell compared to the above-mentioned Pseudo-2D model [27]. Qiu et al. investigated the effects of ambient temperature, charge/discharge rate, and the number of cycles on the polarization characteristics of batteries based on an electrochemical-thermal coupling model [7]. In [16], for planar electrodes of pouchtype lithium-ion batteries, an analytical model was established and the concentrationindependent polarization expression was derived. In addition, the potential and current distribution of the electrodes during the constant-current discharge process were studied, but limited only to the constant-current discharge condition. However, for the Pseudo-2D model, SPM, and their derivatives, parameter identification of the models often requires customization of specific cycle data and cannot be based on real battery cycle data. As in [28], the parameters in improved SPM are identified by frequency response analysis. As well as in [29], seven model parameters are to be identified after rederivation of the original SPM, those seven parameters are divided into three groups, and different customized cycle conditions are designed to identify these three groups of model parameters. In addition, the control equations of this type of model are of very high order and require high calculation power, which is not suitable for real-time applications [30].

An equivalent circuit model which consists of passive components such as resistors, capacitors, inductors, and W-resistors is widely used in BMS [4,18,19,31], but fails to reveal the essential information since passive components do not directly correspond to the internal components or reaction processes inside real batteries [32]. Parameter identification methods for ECMs can be divided into two categories: the time-domain approach and the frequency-domain approach. For example, Li X. et al. used a second-order equivalent circuit model and correlated the small time-constant RC component and the activation polarization and the large time constant RC component with the concentration polarization [19]. The problem is that the two RC components couple with each other, and such a simple distinction does not strictly distinguish between these two types of polarization in a physical sense and on a time scale when using time-domain identification methods. As for the frequency domain approach, the identification of equivalent circuit model parameters using electrochemical impedance spectra is a common approach [11,33–35]. In addition, the DRT approach can distinguish more precisely voltage drops caused by each type of polarization under the frequency domain, and then assign a reasonable time constant to each polarization loss. Furthermore, Zhou et al. combined the DRT method and a physics-based impedance model to separate the solid-phase diffusive polarization voltage drop and the liquid-phase diffusive polarization voltage drop [36]. However, the frequency domain analysis method requires a large input excitation frequency span, and the current BMS on board is not able to meet the requirement. In addition, electrochemical impedance spectroscopy involves a specific *SOC* level [11], and continuous identification during battery cycling is not possible.

Characterizing polarization by changes in cell potential under battery operating conditions is a common method [17], in which the directly measurable cell terminal voltage is used as a measure of characterization accuracy. In this paper, we propose a quantitative battery polarization characterization tool based on a lumped diffusion model (LDM) [32,37] with a joint parameter identification algorithm consisting of the Particle swarm optimization (PSO) algorithm and Levenberg-Marquardt (L-M) method and demonstrate the effectiveness of the proposed method through accuracy verification experiments. Furthermore, a hardware platform was built to demonstrate that the proposed method is capable of quantitative real-time characterization of three types of polarization voltage drops, and has a good tracking performance for the terminal voltage. Compared with the equivalent circuit model, this model preserves the internal physicochemical processes of the battery. Compared with other electrochemical models, this model has fewer parameters to be identified and has good prospects for online applications. The remainder of this paper is organized as follows. First, LDM and the joint parameter identification algorithm are introduced, then the accuracy and effectiveness of the proposed method are verified based on two sets of real-world testing data. Finally, the implementation and results analysis of the hardware platform for online applications are introduced.

#### **2. Battery Modeling and Joint Algorithm Scheme**

#### *2.1. Battery Model Description*

A Pseudo-2D model was proposed by Doyle et al. in 1993 to describe the behavior of lithium-ion batteries based on porous electrode theory and concentrated solution theory [38]. The physicochemical equations that constitute the Pseudo-2D model are the electrochemical reaction process at the critical surface between the active particle surface and the electrolyte solution in both electrode regions, the solid-phase diffusion process, the liquid-phase diffusion process, the solid-phase ohmic resistance, and the liquid-phase ohmic resistance [39,40]. A single particle model was proposed by B. Haran in 1998, which was obtained by further simplifying the assumptions based on the Pseudo-2D model, neglecting the differences in the liquid-phase concentration distribution in the thickness dimension of the electrode sheet, so that one spherical particle can be used to represent the whole electrode [37,41].

The model adopted in this study is LDM, which is a further simplified version of the SPM and Pseudo-2D model, and no longer distinguishes the difference in the spatial distribution of the same type of polarization at different locations of the cell while retaining the control equations reflecting the real physicochemical processes inside the cell. The voltage drop under the battery operating condition is attributed to three types of polarization impedance: ohmic polarization impedance, activation polarization impedance, and concentration polarization impedance. The battery terminal voltage *Ecell* can be obtained by:

$$E\_{cell} = E\_{OCP}(SOC) + \eta\_{olm} + \eta\_{act} + \eta\_{con} \tag{1}$$

where *EOCP*(*SOC*) is the open-circuit voltage, which is a function of the average *SOC* of electrodes. *ηohm*, *ηact*, and *ηcon* are ohmic polarization overpotential, activation polarization overpotential, and concentration polarization overpotential, respectively. *SOC* at a certain moment *SOC*(*t*) can be obtained by the ampere-hour integral method:

$$SOC(t) = SOC(t\_0) + \frac{\int\_{t\_0}^{t} i(t)dt}{Q\_{\text{in}}} \tag{2}$$

where *Qn* is the battery capacity, and *i*(*t*) is the instantaneous current. The time-discrete expression of the above equation takes the form:

$$\text{SOC}(k) = \text{SOC}(0) + \frac{\sum\_{t=0}^{t=k} I(t) \cdot T\_s}{Q\_n} \tag{3}$$

where *Ts* is the sampling period, and *I*(*t*) is the applied current at time point *t*. Ohmic polarization overpotential is defined as:

$$
\eta\_{\text{ohm}} = R\_{\text{ohm}} \cdot I \tag{4}
$$

where *Rohm* is the ohmic resistance. Under the lithium deintercalation/intercalation kinetics assumption on the electrode particle surface in the Pseudo-2D model, the relationship between current density, lithium concentration, and intercalation overpotential is given by the Butler-Volmer formula [42] is expressed in the form:

$$\frac{I}{I\_{1C}} = f\_0 \left( \exp\left(\frac{(1-\alpha)F}{RT}\eta\_{act}\right) - \exp\left(\frac{-\alpha F}{RT}\eta\_{act}\right) \right) \tag{5}$$

where *J*<sup>0</sup> is the dimensionless charge exchange current, used to describe the charge transfer reaction rate on the surface of both electrodes, *I*1*<sup>C</sup>* is the value of the applied current taken at 1 C rate, which is related to the cell capacity, *R* is the molar gas constant, *F* is the Faraday constant, *T* is the reference temperature, and *α* is the charge transfer coefficient. In LDM, the difference in the spatial distribution of current density is neglected and the charge transfer coefficients of both electrodes take the value of 0.5. The overpotential of the two electrodes is considered as a whole, and a single equation expresses the reaction overpotential of the whole cell about the input current excitation:

$$\eta\_{act} = \frac{2RT}{F} \text{asinh} \left(\frac{I}{2f\_0 I\_{1\text{C}}}\right) \tag{6}$$

In this model, the electrode is idealized as a spherical particle, and the electrode local State of Charge *iSOC*(*X*, *t*) varies with both time *t* and the dimensionless spatial variable *X*. The partial differential control equation is obtained by reformulation of Fick's law and solved using a spherically symmetric solution, expressed as:

$$\text{tr}\frac{\partial iSOC}{\partial t} = -\frac{\partial}{\partial X}\left(-\frac{\partial}{\partial X}(iSOC)\right) \tag{7}$$

where *τ* is the diffusion time constant in *s*, and *X* ∈ [0, 1] denotes the dimensionless spatial variable in particle size scale. For this partial differential equation, the initial value condition is

$$i\text{SOC}(X,0) = \text{SOC}\_{0\prime} \quad \text{t} = 0 \tag{8}$$

Both the left and right side boundary conditions are Neuman boundary conditions:

$$\frac{\partial}{\partial X}(iSOC(0,t)) = 0, \quad X = 0 \tag{9}$$

$$\frac{\partial}{\partial X}(iSOC(1,t)) = \frac{\tau I}{N\_{shape}Q\_n}, \quad X = 1 \tag{10}$$

where, for spherical particles, the dimension number *Nshape* takes the value of 3. When State of Charge of the electrode particle surface is defined as *SOCsur f* :

$$SOC\_{surf}(t) = iSOC(1, t), \quad X = 1 \tag{11}$$

The electrode particle average State of Charge *SOCave* reflected the molarity of lithium ions inside the particle, which can be obtained by integrating the local State of Charge *iSOC* over the particle volume with:

$$SOC\_{\text{ave}}(t) = \frac{\int\_0^1 iSOC(X, t) \cdot 4\pi X^2 dX}{\int\_0^1 4\pi X^2 dX} \tag{12}$$

Therefore, the concentration polarization overpotential is expressed as:

$$
\eta\_{\rm conv} = E\_{\rm OCP} \left( \rm SOC\_{surf} \right) - E\_{\rm OCP} \left( \rm SOC\_{\rm avc} \right) \tag{13}
$$

After introducing the above equation, the battery terminal voltage *Ecell* can be reformulated as:

$$E\_{cell} = E\_{OCP} \left( \text{SOC}\_{surf} \right) + \eta\_{olm} + \eta\_{act} \tag{14}$$

### *2.2. Joint Parameter Estimation Algorithm Design*

Three parameters can be identified in the above-mentioned LDM, namely, the ohmic resistance *Rohm* associated with ohmic polarization, the dimensionless charge exchange current *J*<sup>0</sup> associated with activation polarization, and the diffusion time constant *τ* associated with concentration polarization. The objective of the parameter identification algorithm is to find the optimal solution of the state parameters by solving for the minimum of the objective error function so that the model output voltage is as close as possible to the real-world terminal voltage.

Many algorithms were adopted for the identification of model parameters, which in general can be divided into gradient-free methods (i.e., PSO algorithm) and gradient methods (i.e., L-M method). The fitting accuracy of the PSO algorithm is often inferior to that of the L-M method, while the initial value of the L-M method is crucial to determine will fall into a local optimum [43]. Therefore, in this case, we first used the PSO algorithm for the prediction of the model parameters and used the identification results as the initial values of the L-M method to establish a joint algorithm for the accurate estimation of the parameters in LDM.

The PSO algorithm was proposed by Eberhart and Kennedy in 1995 [44]. The population contains a certain number of particles, each of which represents a possible solution. The initial positions of the particles are generally determined randomly. The fitness function associated with those model parameters to be optimized represents the distance of each particle from the optimal solution. The particle velocity determines the direction and distance of each particle's motion during each iteration. Based on this set of rules, particles in the population search for the optimal solution in the solution space.

In the PSO algorithm, the fitness function established based on the above LDM is defined as:

$$f = \text{abs}(\mu\_t(k) - \hat{u}\_t(k))\tag{15}$$

where *ut* is the real-world terminal voltage and *u*ˆ*<sup>t</sup>* is the model output voltage.

In the three-dimensional search space *<sup>S</sup>* <sup>⊆</sup> *<sup>R</sup>*3, the population contains *<sup>P</sup>* particles [1,2, ... ,*p*, ... ,*P*] with the maximum iteration number *G*. The position vector *xp* = *xp*,*τ*, *xp*,*invJ*0, *xp*,*ηIR*,1*<sup>C</sup>* ∈ *S*, and the velocity vector is *vp* = *vp*,*τ*, *vp*,*invJ*0, *vp*,*ηIR*,1*<sup>C</sup>* ∈ *S* for the *p*th particle. Note that to avoid a divide-by-zero error in the calculations, the dimensionless exchange current density *J*<sup>0</sup> is used in the calculations using its inverse *invJ*<sup>0</sup> for the operation and the ohmic resistance *Rohm* is replaced using the ratio of ohmic overpotential at 1 C rate to 1 C rate current *<sup>η</sup>IR*,1*<sup>C</sup> <sup>I</sup>*1*<sup>C</sup>* . Define history optimum *pbestp* as *pp*,*τ*, *pp*,*invJ*0, *pp*,*ηIR*,1*<sup>C</sup>* <sup>∈</sup> *<sup>S</sup>* for the *<sup>p</sup>*th particle, and the global optimum particle *gbest* as *gτ*, *ginvJ*0, *gηIR*,1*<sup>C</sup>* ∈ *S* for the whole population, then the velocity and position update rules for a certain particle are:

$$\begin{cases} \boldsymbol{v}\_{p}^{\text{kg}+1} = \boldsymbol{\omega} \ast \boldsymbol{v}\_{p}^{\text{kg}} + c\_{1} \boldsymbol{r} \text{and} \boldsymbol{d}\_{1} \left( p \text{best}\_{p} - \mathbf{x}\_{p}^{\text{kg}} \right) + c\_{2} \boldsymbol{r} \text{and} \boldsymbol{d}\_{2} \left( \mathbf{g} \text{best} - \mathbf{x}\_{p}^{\text{kg}} \right) \\\ \mathbf{x}\_{p}^{\text{kg}+1} = \mathbf{x}\_{p}^{\text{kg}} + \boldsymbol{v}\_{p}^{\text{kg}+1} \end{cases} \tag{16}$$

where *ω* is the inertia weight, which functions to scale the feasible domain. *c*<sup>1</sup> and *c*<sup>2</sup> are the local learning factor and the global learning factor, respectively. *rand*<sup>1</sup> and *rand*<sup>2</sup> are random numbers uniformly distributed in (0,1). *kg* is the current iteration number. Setting the velocity and position bounds for each particle in the population to ensure that the current velocity and position are restricted to the preset range after each iteration:

$$\upsilon\_p^{k\chi+1} = \begin{cases} \ L B\_\upsilon & \upsilon\_p^{k\chi+1} < L B\_\upsilon \\\ L B\_\upsilon & \upsilon\_p^{k\chi+1} > L B\_\upsilon \end{cases} \tag{17}$$

$$\boldsymbol{\alpha}\_p^{k\text{g}+1} = \begin{cases} \boldsymbol{LB\_x} & \boldsymbol{\alpha}\_p^{k\text{g}+1} < \boldsymbol{LB\_x} \\ \boldsymbol{UB\_x} & \boldsymbol{\alpha}\_p^{k\text{g}+1} > \boldsymbol{LB\_x} \end{cases} \tag{18}$$

When facing constrained optimization problems with multi-peak distribution, the following two improvements are applied to this case to avoid the identification algorithm falling into local optimum, and to improve the global search capability.

(i) The decreasing time-varying inertia weight *ωkg* ∈ [*ωmin*, *ωmax*] is introduced, and *ωkgvkg* represents the momentum of particle motion in the population. In this case, the inertia weight decreases uniformly as the number of iterations increases. The purpose of this is to make sure that the particles have good global search ability at the beginning to avoid falling into a local optimum, and at the end of the iteration to facilitate local search and accelerate convergence, which achieves a good balance between convergence efficiency and global searchability.

$$
\omega\_{\rm kg} = \omega\_{\rm max} + kg \ast \frac{\omega\_{\rm max} - \omega\_{\rm min}}{G} \tag{19}
$$

Furthermore, the local learning factor and global learning factor are set as a function of the time-varying inertia weights, defined as:

$$\begin{cases} c\_{1,k\text{g}} = k\_1 \ast \begin{pmatrix} 1 - \omega\_{k\text{g}} \\ 1 - \omega\_{k\text{g}} \end{pmatrix} \\\ c\_{2,k\text{g}} = k\_2 \ast \begin{pmatrix} 1 - \omega\_{k\text{g}} \end{pmatrix} \end{cases} \tag{20}$$

where *k*<sup>1</sup> and *k*<sup>2</sup> are the learning factor gain, generally take *k*<sup>2</sup> > *k*1. *c*1,*kg* and gradually increase with the increase of the number of iterations. The motion of particles in the early stage is less influenced by the history and other particles to enhance the global searchability, while in the later stage particles are increasingly influenced by the history and other particles to accelerate convergence. Then, update the velocity and position of the *p*th particle by the rules below:

$$\begin{cases} \boldsymbol{\upsilon}\_{p}^{\text{kg}+1} = \boldsymbol{\omega}\_{\text{kg}} \* \boldsymbol{\upsilon}\_{p}^{\text{kg}} + \boldsymbol{c}\_{1,\text{kg}} \boldsymbol{r} \boldsymbol{n} \boldsymbol{d}\_{1} \left( \boldsymbol{p} \text{best}\_{p} - \boldsymbol{\mathsf{x}}\_{p}^{\text{kg}} \right) + \boldsymbol{c}\_{2,\text{kg}} \boldsymbol{r} \boldsymbol{n} \boldsymbol{d}\_{2} \left( \boldsymbol{g} \text{best} - \boldsymbol{\mathsf{x}}\_{p}^{\text{kg}} \right) \\\ \boldsymbol{\mathsf{x}}\_{p}^{\text{kg}+1} = \boldsymbol{\mathsf{x}}\_{p}^{\text{kg}} + \boldsymbol{\upsilon}\_{p}^{\text{kg}+1} \end{cases} \tag{21}$$

(ii) In each iteration, a certain probability of particle mutation occurrence is set, and the position information of the selected particles is reassigned to ensure that some of the particles in the population can jump out of the local optimum trap and continue to search for other possible global optimum solutions.

$$\begin{cases} \begin{array}{c} \text{x}\_{p}^{\text{kg}}(\dim) = LB\_{\text{x}}(\dim) + (\text{LB}\_{\text{x}}(\dim) - LB\_{\text{x}}(\dim)) \ast r\_{1} \\ \text{dim} = \text{ceil}(3 \ast r\_{2}) \end{array} \end{cases} \tag{22}$$

where *x kg <sup>p</sup>* (*dim*) ⊆ [*x kg <sup>p</sup>*,*τ*, *x kg <sup>p</sup>*,*invJ*0, *x kg <sup>p</sup>*,*R*] represent the randomly selected dimension among those three dimensions in the position vector for the *p*th particle, and *r*<sup>1</sup> and *r*<sup>2</sup> are random numbers uniformly distributed on the interval (0,1), respectively. Among the population, update *pbestp* and *gbest* after each iteration by the equations below, until the maximum iteration number is met and the value of *gbest* will be the final solution.

$$\begin{cases} \begin{array}{c} pbest\_p = \min \left[ f\_{p\prime} pbest\_p \right] \\ gbest = \min \left[ f\_{p\prime} gbest \right] \end{array} \end{cases} \tag{23}$$

L-M method [45,46] is a classical numerical solution method for solving the minimum of nonlinear equations, which retains both the stability of the steepest descent method and the fast convergence property of the Gaussian Newton method. In this paper, we set the parameter vector *θ* = [*τ*, *invJ*0, *ηIR*,1*C*], the original data is *M* sets of battery cycling data (*um*, *im*), *m* = 1, 2, ... , *M*, and the output voltage of LDM is *u*ˆ(*im*, *θ*). Then, the error function of the L-M method is expressed as:

$$E(\mathfrak{a}) = \sum\_{m=1}^{M} e^2 = \sum\_{m=1}^{M} \left(\mathfrak{u}\_m - \mathfrak{A}(i\_{m\prime}\theta)\right)^2 \tag{24}$$

where *e* is the voltage error for one set of data. The optimal model parameter vector *θ* is obtained by iteratively solving for the minimum of the above error function, and the iterative expression of the L-M method is:

$$\theta\_{k+1} = \theta\_k + \left[J^T \cdot J\_k + \mu L\right]^{-1} \cdot J\_k \cdot e(k) \tag{25}$$

where *μ* is the damping factor, *L* is the identity matrix, and *k* is the current iteration number. *Jk* is the Jacobian matrix:

$$J\_k = \begin{bmatrix} \frac{\partial c\_1}{\partial \theta [1]} & \frac{\partial c\_1}{\partial \theta [2]} & \frac{\partial c\_1}{\partial \theta [3]} \\ \frac{\partial c\_2}{\partial \theta [1]} & \frac{\partial c\_2}{\partial \theta [2]} & \frac{\partial c\_2}{\partial \theta [3]} \\ \vdots & \vdots & \vdots \\ \frac{\partial c\_M}{\partial \theta [1]} & \frac{\partial c\_M}{\partial \theta [2]} & \frac{\partial c\_M}{\partial \theta [3]} \end{bmatrix} \tag{26}$$

The iteration termination conditions are:

$$(k > k\_{\text{max}})OR(e(k) < \varepsilon) \tag{27}$$

where *kmax* is the maximum iteration number, and *ε* is the tolerance. The procedure of the proposed joint algorithm is summarized in Figure 1.


**Figure 1.** Major steps of the joint parameter identification algorithm.
