Balancing-Aware Charging Strategy For Series-Connected Lithium-Ion Cells: A Nonlinear Model Predictive Control Approach

TL;DR

This paper introduces a nonlinear model predictive control strategy for balancing and optimizing the charging of series-connected lithium-ion cells, using an electrochemical model and validated through simulations.

Contribution

It presents a novel nonlinear MPC approach tailored for balancing-aware charging of lithium-ion cells, with practical implementation considerations.

Findings

Effective balancing achieved in simulations

Validated approach on commercial cells

Potential for improved battery lifespan

Abstract

Charge unbalance is one of the key issues for series-connected Lithium-ion cells. Within this context, model-based optimization strategies have proven to be the most effective. In the present paper, an ad-hoc electrochemical model, tailored to control purposes, is firstly presented. Relying on this latter, a general nonlinear MPC for balancing-aware optimal charging is then proposed. In view of the possibility of a practical implementation, the concepts are subsequently specialized for an easily implementable power supply scheme. Finally, the nonlinear MPC approach is validated on commercial cells using a detailed battery simulator, with sound evidence of its effectiveness.

Tables

Table 1. TABLE I: Initial capacity ( C 0 superscript 𝐶 0 C^{0} ), film resistance ( R s ​ e ​ i 0 superscript subscript 𝑅 𝑠 𝑒 𝑖 0 R_{sei}^{0} ) and normalized state of charge ( z 0 superscript 𝑧 0 z^{0} ), values for the different cells of the pack.

Parameter	Cell 1	Cell 2	Cell 3	Cell 4	Cell 5	Cell 6
$C^0$ [Ah]	$6.47$	$7.62$	$8.83$	$8.68$	$8.46$	$7.22$
$R_{sei}^0$ [$m\Omega$]	$2.09$	$2.36$	$2.22$	$2.36$	$2.12$	$1.76$
$z^0$	$0.39$	$0.19$	$0.21$	$0.20$	$0.34$	$0.37$

Table 2. TABLE II: Parameters of the NMPC optimization.

$𝑵$	$𝑯$	${𝑻}_{𝒔}$	$𝒂$	$\overline{𝒛}$	${𝑰}_{𝒎𝒂𝒙}$
$6$	$3$	$10s$	$5$	$1$	$7.5A$
${𝑻}_{𝒎𝒂𝒙}$	${𝑽}_{𝒎𝒂𝒙}$	${𝑷}_{𝒎𝒂𝒙}$	${𝜶}_{\mathrm{𝟏}}$	${𝜶}_{\mathrm{𝟐}}$	${𝜶}_{\mathrm{𝟑}}$
$313.15K$	$4.2V$	$0.75V_{max}{I}_{max}$	${10}^4$	$25$	$0$
${𝜶}_{\mathrm{𝟒}}$	${𝜶}_{\mathrm{𝟓}}$	${𝜶}_{\mathrm{𝟔}}$	${𝜶}_{\mathrm{𝟕}}$	${𝑻}_{𝒆𝒏𝒗}$	${𝜶}_{{𝒔}_{𝒊}},\forall i$
$1$	${10}^{-3}$	${10}^{-3}$	${10}^5$	$298.15$	${10}^{15}$

Equations

$$\psi (T(t))={\psi }_0{e}^{\frac{-E_{a,\psi }}{RT(t)}}$$

$${\bar{\theta }}_i(t)=\frac{{\bar{c}}_{s,i}(t)}{c_{s,i}^{max}}$$

$${\bar{\theta }}_n(t)={\theta }_n^{0\%}+\frac{{\bar{\theta }}_p(t)-{\theta }_p^{0\%}}{{\theta }_p^{100\%}-{\theta }_p^{0\%}}({\theta }_n^{100\%}-{\theta }_n^{0\%})$$

$${\dot{\bar{\theta }}}_p(t)=-\frac{\Delta {\theta }_p}{C(t)}{I}_{app}(t)$$

$${\dot{\bar{q}}}_i(t)=-30\frac{D_{s,i}(T(t))}{R_{p,i}^2}{\bar{q}}_i(t)-\frac{45\Delta {\theta }_i{c}_{s,i}^{max}}{6R_{p,i}C(t)}{I}_{app}(t)$$

$${\theta }_i(t)={\bar{\theta }}_i(t)+\frac{8R_{p,i}{\bar{q}}_i(t)}{35c_{s,i}^{max}}-\frac{R_{p,i}^2\Delta {\theta }_i}{105D_{s,i}(T(t))C(t)}{I}_{app}(t)$$

$$z(t)=\frac{{\bar{\theta }}_n(t)-{\theta }_n^{0\%}}{{\theta }_n^{100\%}-{\theta }_n^{0\%}}$$

$${ϵ}_p^{act}(t)$$

$${ϵ}_n^{act}(t)$$

$$\begin{array}{rl}{ϵ}_p\frac{\partial c_{e,p}(x,t)}{\partial t}=& \frac{\partial }{\partial x}\left[D_{e,p}^{eff}(T(t))\frac{\partial c_{e,p}(x,t)}{\partial x}\right]\\ & -\frac{1-t_{+}}{FA{L}_p}{I}_{app}(t)\end{array}$$

$${ϵ}_s\frac{\partial c_{e,s}(x,t)}{\partial t}=$$

$$\begin{array}{rl}{ϵ}_n\frac{\partial c_{e,n}(x,t)}{\partial t}=& \frac{\partial }{\partial x}\left[D_{e,n}^{eff}(T(t))\frac{\partial c_{e,n}(x,t)}{\partial x}\right]\\ & +\frac{1-t_{+}}{FA{L}_n}{I}_{app}(t)\end{array}$$

$$\frac{\partial c_{e,p}(0_p,t)}{\partial x}=$$

$$D_{e,p}^{eff}(T(t))\frac{\partial c_{e,p}(L_p,t)}{\partial x}=$$

$$D_{e,s}^{eff}(T(t))\frac{\partial c_{e,s}(L_s,t)}{\partial x}=$$

$$c_{e,p}(L_p,t)=$$

$$c_{e,p}(L_n,t)=$$

$$\begin{array}{rl}{ϵ}_p\frac{\partial c_{e,p}^{[k]}(t)}{\partial t}=& {\left[\frac{{\tilde{D}}_e(x,T(t))}{\Delta x_p}\frac{\partial c_{e,p}(x,t)}{\partial x}\right]|}_{x_{j,\underline{k}}}^{x_{j,\bar{k}}}\\ & -\frac{1-t_{+}}{FA{L}_p}{I}_{app}(t)\end{array}$$

$${ϵ}_s\frac{\partial c_{e,s}^{[k]}(t)}{\partial t}=$$

$$\begin{array}{rl}{ϵ}_n\frac{\partial c_{e,n}^{[k]}(t)}{\partial t}=& {\left[\frac{{\tilde{D}}_e(x,T(t))}{\Delta x_n}\frac{\partial c_{e,n}(x,t)}{\partial x}\right]|}_{x_{j,\underline{k}}}^{x_{j,\bar{k}}}\\ & +\frac{1-t_{+}}{FA{L}_n}{I}_{app}(t)\end{array}$$

$${\tilde{D}}_e(x,T(t))=\{\begin{array}{ll}{D}_{e,1}& \mathrm{if}x\in \{x_{p,\underline{M}},x_{s,\overline{1}}\}\\ D_{e,2}& \mathrm{if}x\in \{x_{s,\underline{M}},x_{n,\overline{1}}\}\\ D_{e,j}^{eff}(T(t))& \mathrm{otherwise}\end{array}$$

$$D_{e,1}$$

$$D_{e,2}$$

$$\mathcal{H}\left({\rho }_1,{\rho }_2,{\lambda }_1,{\lambda }_2\right)=\frac{{\rho }_1{\rho }_2({\lambda }_1+{\lambda }_2)}{{\rho }_1{\lambda }_2+{\rho }_2{\lambda }_1}$$

$$\dot{C}(t)$$

$${\dot{R}}_{sei}(t)$$

$${\bar{j}}_{side}(t)=-\frac{i_{0,side}(t)}{F}{e}^{\left(\frac{0.5F}{RT(t)}{\bar{\eta }}_{side}(t)\right)}$$

$$i_{0,side}(t)=\{\begin{array}{ll}{i}_{0,base}{\left(\frac{|I_{app}(t)|}{I_{1C}}\right)}^w& \mathrm{if}{I}_{app}(t)\le 0\\ 0& \mathrm{otherwise}\end{array}$$

$${\bar{\eta }}_{side}(t)={\bar{\eta }}_n(t)+{\bar{U}}_n(t)-U_{SEI}$$

$${\bar{\eta }}_i(t)=\frac{2RT(t)}{F}{\sinh }^{-1}\left(\frac{\Delta {\theta }_{i}F{R}_{p,i}}{6{\bar{i}}_{0,i}(t)C(t)}{I}_{app}(t)\right)$$

$${\bar{i}}_{0,i}(t)=F{k}_i(T(t))\sqrt{{\bar{c}}_{e,i}(t){\theta }_i(t)(1-{\theta }_i(t))}$$

$${\bar{c}}_{e,i}(t)=\frac{1}{M}\sum _{k=1}^M{c}_{e,i}^{[k]}(t)$$

$$\begin{array}{rl}V(t)=& {\bar{U}}_p(t)-{\bar{U}}_n(t)+{\bar{\eta }}_p(t)-{\bar{\eta }}_n(t)+\\ & +\Delta {\Phi }_e(t)+I_{app}(t)R_{sei}(t)\end{array}$$

$$\Delta {\Phi }_e(t)={\Phi }_e^{drop}(t)+\frac{2RT(t)}{F}(1-t_{+})\log \left(\frac{c_{e,p}^{[1]}}{c_{e,n}^{[M]}}\right)$$

$${\Phi }_e^{drop}(t)\simeq -\frac{I_{app}(t)}{2M}\left({\varphi }_p(t)+2{\varphi }_s(t)+{\varphi }_n(t)\right)$$

$${\varphi }_p(t)$$

Full text

I Introduction

In the last few decades, the ever increasing demand for technologies which rely on portable energy sources (e.g. electric and hybrid vehicles, consumer electronics, microgrids and IoT-related devices) has determined a growth in the request for high-performance energy accumulators. Among the possible manufacturing chemistries, Lithium-ion (Li-ion) batteries have proven to be the most promising [1, 2, 3]. In order to provide high terminal voltage, total capacity and available power, battery packs are usually composed of several cells, arranged in series and parallel connections. Due to unavoidable inconsistencies in the manufacturing process and non-homogeneities in the operating conditions (e.g. temperature and ageing effects, which appear over time), the cells composing a battery pack exhibit slightly different features (e.g. in terms of internal impedance, self-discharge rate and physical volume) [4, 5]. As a consequence, appreciable unbalancing in the stored charge arises already after few charge/discharge cycles if conventional charging methods are employed [6]. This latter is a limiting factor in terms of a complete battery exploitation (undercharging), state of health preservation (premature and non-uniform wear) and safety (overcharging and overdischarging phenomena). In particular, undercharging is a consequence of an early interruption of the charging procedure due to the reaching of the voltage threshold from a subset of cells. Premature and non-uniform degradation, instead, comes mainly from overvoltage exposures and high temperature spikes, which also affect the self-discharge rates. Moreover, during operation, temperature gradients arise between the different cells and even inside of each of them. All of these facts, in turn, contribute to an increase in the differences among the cells characteristics, which eventually lead to greater unbalance.

In order to overcome such issues, a huge research effort has been made towards the development of algorithms implementing the idea of State of Charge (SOC) equalization. Two main categories can be identified for such balancing strategies, namely passive and active methods [7, 8, 9]. The former are based on energy dissipation, while the latter, in order to limit energy losses, rely mainly on charge redistribution . Among all of them, different structures can be adopted, e.g. cell bypass, cell-to-cell, cell-to-pack, pack-to-cell, cell-to-pack-to-cell [10].

At an earlier research stage, equalization was carried out using only model-independent low-level control strategies. Dissipative approaches have been proposed, e.g., in [11, 12, 13], while active voltage-based equalization methods have been presented, e.g., in [14, 15, 16, 17]. Although the voltage can easily be measured in practice, some works have highlighted that a performance improvement can usually be achieved relying on capacity-based algorithms[18, 19], which however require an estimation of the state of charge, e.g. through observations [20, 21, 22, 23]. This approach implies the need for a mathematical model of the cells, whose accuracy directly affects the overall algorithms performance. Several cell models have been proposed in the literature, which can mainly be classified into two categories, namely Equivalent Circuit Models (ECMs) and Electrochemical Models (EMs). While those in the first set give a simplified equivalent description, which only includes lumped states and parameters, the ones in the second group provide an accurate representation of the internal phenomena at the expense of heavier computational burden.

The most used electrochemical model is the well-known Pseudo-Two-Dimensional (P2D) model [24, 25], also referred as Doyle-Fuller-Newman (DFN) model, which consists of tightly coupled and highly nonlinear Partial Differential Algebraic Equations (PDAEs). Note that, due to their complexity, such PDAEs require proper advanced integration methods. Furthermore, it has been shown in the literature that some parameters of the DFN model are unidentifiable without invasive experiments [26, 27, 28] and that some observability issues [29] are present. For these reasons, it appears evident that the P2D model is more suited for simulation purposes rather than control and estimation. In order to preserve physical meaning and a sufficiently high accuracy while providing observability and identifiability, electrochemical reduced order models can be used [30, 31, 32, 33]. For instance, the Single Particle Model (SPM) [34, 35] and the Single Particle Model with Electrolyte dynamics (SPMe) [36] are derived from the P2D considering each electrode as a single spherical particle. The resulting description is greatly simplified but still includes Partial Differential Equations (PDEs) related to the diffusion of the ions in the solid phase (and also in the electrolyte in case of SPMe), which partially preserve the original computational complexity.

Relying on a mathematical representation of the cells behaviour given by one of the available models, optimization-based balancing approaches can be developed. Model Predictive Control (MPC, [37, 38]) is an optimization-based control technique which, due to its ability to deal with complex systems in an optimal and effective way, considering also constraints, appears particularly suitable within the context of Lithium-ion batteries balancing. For instance, in [39] balancing is carried out considering an arbitrary connection topology for the cells, which are modeled as first-order integrators. A linear SOC model is employed by the authors in [40] for minimum-time balancing of a battery pack. In [41] an adaptive MPC is proposed and a slightly more accurate model is considered in order to redistribute energy during the charging procedure. Even if considerable improvements have been shown in the literature, few issues are still evident. As a matter of fact, almost all the proposed works rely on relatively simple cell models. It appears then difficult to assess the validity of those algorithms that are tested only in simulation with such simple models considered also as the real plants. These latter are indeed expected to miss the representation of some key features, at least in particular situations. For instance, temperature as well as ageing dynamics are often neglected, even if they lead to performance degradation over time. In some high-power applications, e.g. electric and hybrid vehicles, cooling systems are arranged in order to dissipate the excessive heat. In such cases, thermal couplings between cells and between cells and coolant play an important role which should be accounted for in the equalization procedure.

The majority of the equalization algorithms are designed to operate offline, at time instants defined by the specific policies, with the aim to mitigate the unbalancing effects of different cells characteristics. This seems hardly applicable in a practical scenario since it usually requires an idle time in which only the charge balancing is performed (except for the cases in which it is performed during the discharge and the load current profile is known in advance, as shown e.g. in [42]). Furthermore, the concept of SOC balancing is effective in maximizing the battery utilization if one supposes that the total capacity and the Solid-Electrolyte Interphase (SEI) resistance are always equal among the cells throughout the battery lifetime. In a real framework, however, due to all the previously mentioned phenomena, this is never the case. Therefore, an offline algorithm such as that usually utilized, not only is unable to maximize battery exploitation, but can also lead to overcharge in the successive charging cycle if a standard method is adopted (e.g. Constant Current - Constant Voltage, CC-CV). In fact, if one balances the SOC at a value lower than the maximum, right before a charge, then the SOCs will diverge one from each other due to the difference in the capacities and the application of a unique current. The only reasonable solution in such cases appears to be the so-called top balancing [43], which consists in balancing the cells SOC to their maximum value at the end of the charge.

In the present paper a general charging strategy based on Nonlinear MPC (NMPC) is proposed with the aim of overcoming the mentioned problems. Then, it is specialized for the case in which each cell is embedded into a power unit (PU) along with two switches, which allow for the cell bypass. Starting from the SPMe, an even more simplified description suitable for control purposes is derived, based only on Ordinary Differential Equations (ODEs). Beside a lumped internal temperature model, already present in the literature [44, 45, 46], the thermal couplings between the cells composing the pack is modelled for an arbitrary topological disposition. Furthermore, in order to enhance the description, coolant temperature as well as ageing dynamics are also considered. This latter accounts for capacity losses and film resistance growth over time, whereas the former can be particularly important in automotive applications. Considering all the features mentioned above in a unique formulation represents a novelty to the best of authors knowledge.

In addition, nullifying the unbalance at the end of each charge cycle, the proposed approach is designed to be able to replace the presence of an equalization system coupled with the battery pack. Due to the simple structure of the power supply control scheme adopted for the practical implementation proposal, which only consists in a variable current generator and a set of PUs, energy losses are maintained at a reasonable level. Furthermore, the predictive nature of the algorithm makes it possible to avoid redistribution, thus increasing energy efficiency and reducing charging time. Suitably designed constraints are enforced for both safety and battery health preservation reasons. The effectiveness of the presented strategy is validated in simulation using an extended version (available on request) of the LIONSIMBA toolbox [47], which is a numerical implementation of the P2D. The simulation environment is set considering a full order model for the solid diffusion law and a lumped model for the thermal dynamics. In order to consider a real world scenario the cells parameters are taken from [48, 49], where a complete parametrization of the Kokam SLPB 75106100 cell is presented.

With the exception of the temperature, the internal states of each cell are not measurable in practice. Nevertheless, in this paper the assumption of availability of all the relevant states is made. The development of a suitable observer and its coupling with the presented control system goes beyond the scope of this paper and will be investigated in the future. This work focuses only on series-connected cells, similarly to what done in the literature. Research effort will be devoted to extend the applicability also to parallel-connected series of cells.

I-A Main Contributions

The main contributions are here summarized for the readers convenience. First of all, a simplified electrochemical model which is suitable for online control but also able to describe all the relevant aspects of a Lithium-ion battery pack is formulated. In particular, beside the internal dynamics of each cell, the introduced model includes coupled thermal dynamics (consisting also of the cooling system) as well as ageing effects. Basing on this description, a general NMPC formulation is proposed for an optimal charging strategy able to achieve balancing at the end of each charge (top balancing). Additionally, the procedure is applied to a proposed supply scheme, easy to implement but more suited for mixed-integer optimization methods. Due to the superior scalability of continuous optimization over the discrete one, particular adaptations are proposed for the coupling of the scheme with the NMPC controller.

I-B Outline

The remaining of the paper is organized as follows. In Section II a detailed explanation of the model components is given. In II-A the simplified version of the SPMe for the single cell electrochemical behaviour is introduced, along with ageing and lumped thermal models. In II-B the thermal couplings are described through a structure able to account for arbitrary relationships between the elements. Then, the model performance is assessed in II-C, through a sinusoidal charge-discharge input current. In Section III, the series balancing-aware optimal charge problem is formally stated (III-A) and applied to the proposed power supply circuit (III-B). Section IV is devoted to the evaluation of the proposed work. In particular, the testbed is described in IV-A, where the characteristics of the considered cells and the battery pack structure are presented. All the cells exhibit different initial conditions in terms of capacity, SEI resistance and SOC, in order to emulate a real scenario in which the non-homogeneities have led to an unbalanced situation. Firstly, in IV-B, a standard charging protocol (namely, the CC-CV) is applied in order to show the effects of such an initial mismatch. Then, a simple method that is based on the proposed circuit but only exploits voltage measurements is tested in IV-C so as to highlight the need for a model-based approach in order to achieve high performance. Eventually, the effectiveness of the proposed strategy is made evident in IV-D. Conclusions are finally drawn in Section V, where possible future enhancements are also discussed.

II Battery Model

In this section, a control-oriented electrochemical model of the battery pack with thermal and ageing dynamics is presented, which is sufficiently accurate but also suitable for online control applications. In II-A the electrochemical behaviour of a single cell is considered along with the respective temperature dynamics and ageing effects. The thermal coupling between cells and coolant is modelled in II-B, where an arbitrary topology of the elements inside the battery pack is considered. At the end (II-C), a sinusoidal charge-discharge cycle is performed to give an idea of the general accuracy of the model with respect to the P2D implemented in LIONSIMBA.

II-A Model of a Single Cell

The SPMe, proposed in [36], provides an approximation of the DFN full order model which still relies on PDEs for describing the diffusion of solid (through the Fick’s law) and electrolyte concentrations. In this paper, we consider a simplified version of the SPMe in which the Fick’s law is reduced for each electrode to an ODE by considering a polynomial approximation of the solid concentration along the particle radius [50]. Furthermore, the electrolyte diffusion PDE is spatially discretized according to the Finite Volume (FV) method [51], as already done in a similar context in, e.g., [47]. In addition, starting from the works in [52, 53] we augment the adopted model with the cell ageing effects in terms of SEI resistance growth and capacity fade, which are among the main degradation phenomena for a Lithium-ion cell [54]. Note that, in the proposed model, some of the parameters exhibit a temperature dependency. In this work, this latter will be modelled according to the Arrenhius law

[TABLE]

where $T(t)$ is the temperature, $E_{a,\psi}$ is the activation energy associated with the generic parameter $\psi(T(t))$ and $R$ is the universal gas constant. Such relationship will be explicitly highlighted for each of the temperature-dependent parameter.

In this Subsection, the index $j\in\{p,s,n\}$ will refer to the three cell sections (respectively, cathode, separator and anode), while $i\in\{p,n\}$ will be used in equations valid only for the electrodes.

II-A1 Solid Concentration Dynamics

According to the authors in [50], the lithium concentration profile along the radial coordinate can be approximated as a fourth-order polynomial function of $r$, $0\leq r\leq R_{p,i}$, with $R_{p,i}$ being the particle radius. The parameters involved in this approximation depend only on the average stoichiometry $\bar{\theta}_{i}(t)$ and the volume-averaged concentration flux $\bar{q}_{i}(t)$. The dynamics of these latter can be described in terms of the ODEs recalled below.

The average stoichiometry $\bar{\theta}_{i}(t)$ is defined, for each electrode, as the ratio between the average solid concentration $\bar{c}_{s,i}(t)$ and the maximum solid concentration $c_{s,i}^{max}$, i.e.

[TABLE]

In order to minimize the number of ODEs required to describe the lithium solid concentration within each cell, we express the anodic stoichiometry in terms of the cathodic one as in [31]

[TABLE]

where $\theta^{100\%}_{i}$ and $\theta^{0\%}_{i}$ are the values of the stoichiometry respectively at the battery charged and discharged states (provided in the cells data-sheet). This relationship is guaranteed by the fact that the moles of lithium in the solid phase are preserved [55]. The dynamics of the average stoichiometry for the cathode can be expressed as

[TABLE]

where $C(t)$ is the available capacity of the cell, $\Delta\theta_{i}=\theta_{i}^{100\%}-\theta_{i}^{0\%}$, and $I_{app}(t)$ is the applied current. Note that, the convention adopted in this work is such that the charging current is negative, i.e. $I_{app}(t)<0$. The volume-averaged concentration flux can be described by

[TABLE]

where $D_{s,i}(T(t))$ is the temperature-dependent solid diffusion coefficient. Then, the positive and negative surface stoichiometries are given respectively by

[TABLE]

Finally, the normalized SOC $z(t)$ is defined as follows

[TABLE]

Note that Equation (7) relies on the widely used convention that the cell SOC corresponds to the anodic one. Moreover, the presented equations are formulated in order to keep (4) intuitive and simplifying a bit the notation. In particular, they are derived from the standard SPM description by taking the active material volume fraction $\epsilon^{act}_{i}(t)$ as a function of the available capacity of the cell $C(t)$, i.e.

[TABLE]

where $A$ is the contact area between solid and electrolyte phase, $F$ is the Faraday constant and $L_{i}$ is the layer thickness. This definition is chosen according to the fact that the SOC-OCV (Open Circuit Voltage) curve variations are small over the battery lifetime [22].

Remark 1

The polynomial approximation used above to describe the lithium solid concentration along the particle radial coordinate can be replaced by a more accurate one, if needed. For instance, a spatial discretization using Chebyshev orthogonal collocation can be employed [23].

II-A2 Electrolyte Concentration Dynamics

The PDEs governing the diffusion of the electrolyte concentration $c_{e,j}(x,t)$ in the SPMe [36] are the following

[TABLE]

where $D_{e,j}^{eff}(T(t))=D_{e}(T(t))\epsilon_{j}^{p_{j}}$ with $\epsilon_{j}$ the material porosity, $p_{j}$ the Bruggeman coefficient, $t_{+}$ the transference number and $D_{e}(T(t))$ the temperature-dependent diffusion coefficient within the electrolyte. The boundary conditions are specified by the following equations [36]

[TABLE]

where $0_{j}$ and $L_{j}$ are, respectively, the beginning and the end of the $j$-th section along the $x$-axis.

In this work, the spatial domain is divided for each section into $M$ non-overlapping volumes with centered nodes. The $k$-th volume, with $k=1,\cdots,M$, of the $j$-th section is centred at the spatial coordinate $x_{j,k}$ and spans the interval $\Omega_{j,k}=\left[x_{j,\bar{k}},x_{j,\underline{k}}\right]$, whose width is $\Delta x_{j}=L_{j}/M$. The PDEs are then discretized according to the FV method. Defining $c_{e,j}^{[k]}(t)$ as the average electrolyte concentration over the $k$-th volume of $j$-th section, one has

[TABLE]

where the terms are evaluated as explained in details in [47]. In particular, the electrolyte diffusion coefficients are computed as follows

[TABLE]

with

[TABLE]

where $\mathcal{H}$ is the harmonic mean operator, defined as

[TABLE]

II-A3 Ageing Effects

The main consequences of ageing in Lithium-ion cells are a loss of capacity and a growth of the SEI film resistance [54]. Such effects are primarily caused by high cycling currents and temperatures. In order to describe these degradation phenomena, we modify the models in [52, 53] to match the SPMe dynamics. In particular, the evolution of the residual available capacity $C(t)$ and the SEI resistance $R_{sei}(t)$ is as follows

[TABLE]

where $M_{w}$ is the molar weight of the negative electrode, $\rho_{n}$ the material density and $\nu$ the admittance of the film. The average side reaction flux is

[TABLE]

while the side reaction exchange current $i_{0,side}(t)$ is

[TABLE]

where $i_{0,base}$ is the base-side reaction current, $I_{1C}$ is the current corresponding to a 1C rate for the considered cell and $w$ is an empirically obtained coefficient. Note that no degradation occurs when the cell is discharged (i.e. for $I_{app}(t)>0$). The side reaction overpotential $\bar{\eta}_{side}(t)$ is given by

[TABLE]

where $U_{SEI}$ is the reference potential for the SEI side reaction, while $\bar{U}_{n}(t)$ and $\bar{\eta}_{n}(t)$ are the negative Open Circuit Potential (OCP) and overpotential, respectively. In particular, the OCPs, for both electrodes, can be represented as nonlinear functions of the surface stoichiometric coefficients with expressions dependent on the particular considered cell (see IV-A). On the other side, $\bar{\eta}_{i}(t)$ is given by

[TABLE]

with

[TABLE]

where $k_{i}(T(t))$ is the temperature-dependent rate reaction constant, and $\bar{c}_{e,i}(t)$ is the average electrolyte concentration in the $i$-th section, approximated as follows

[TABLE]

II-A4 Output Voltage

The terminal voltage is given by

[TABLE]

with elements defined according to equations (16b), (20), and

[TABLE]

in which $\Phi_{e}^{drop}(t)$ is approximated according to the assumption that the ionic current $i_{e}(x,t)$ has a trapezoidal shape over the spatial domain [36] (see Fig. 1 for a schematic representation), i.e.

[TABLE]

In view of the mentioned assumption one can easily derive

[TABLE]

where $\kappa(c_{e,j}^{[k]}(t),T(t))$ is the temperature-dependent electrolyte conductivity for the $k$-th volume of the $j$-th section, which is usually expressed as an empirically derived nonlinear function of the electrolyte concentration in that volume (see IV-A). Such relationship can in principle differ from cell to cell. Note that, the expressions in (26) are computed considering different electrolyte conductivity values for each one of the discrete volumes constituting the cell sections. This approximantion provides better accuracy with respect to e.g. [36], in which the electrolyte conductivity is assumed constant over the spatial domain.

II-B Battery Thermal Model

With a slight abuse of notation, from now on the variables $i$ and $j$ are used to indicate specific cells inside the battery pack. This latter is composed of $N$ cells, which exchange heat through direct contact with other cells (conduction) and possibly with a coolant fluid (by convection). Both conduction and convection phenomena are approximated using lumped thermal resistance parameters $R_{i,j}=R_{j,i}$, which are associated with the couple of cells $(i,j)$, and $R_{i}^{sink}$ that considers the interactions of the $i$-th cell with the coolant. As a matter of fact, any arbitrary disposition of the cells can be fully described from a thermal point of view through the choice of the resistance parameters. Having defined $C_{th,i}$ and $C_{th,sink}$ as the thermal capacities of the $i$-th cell and the sink respectively, the temperature dynamics of the $i$-th cell is given by

[TABLE]

where $T_{i}(t)$ and $T_{sink}(t)$ are the temperatures of the $i$-th cell and the coolant, $\Delta T_{j,i}(t)=T_{j}(t)-T_{i}(t)$ and $\Delta T_{i}^{sink}(t)=T_{i}(t)-T_{sink}(t)$. The nonlinear term $Q_{i}(t)$ is the heat generation due to polarisation inside the cell, defined as

[TABLE]

where the terms in (28) are the instantaneous current, voltage and open circuit potentials associated with the $i$-th cell. The temperature dynamics of the sink is instead given by

[TABLE]

where $\xi(t)$ is the thermal power released by the cooling system to the external environment.

II-C Performance Assessment

Since the model described above is a modified version of the SPMe which includes also the ageing effects, in order to test its accuracy, we performed a comparison with the more detailed P2D model. In particular, for this latter, we rely on the LIONSIMBA toolbox and consider the ageing dynamics as in [35]. In order to evaluate the performance also at high-frequency, the models are excited with the sinusoidal charge-discharge input current depicted in Fig. 2. The resulting voltage, temperature and SOC are plotted respectively in Fig. 3, 4 and 5, while the ageing effects are depicted in Fig. 6. For the sake of clarity, in this latter, only the differences from the initial values are considered. The obtained results highlight the accuracy of the proposed model.

III Balancing Algorithm

In this section, we present a nonlinear MPC scheme suitable for guaranteeing a balancing-aware charging for Lithium-ion cells. In particular, in this work we consider series-connected cells. This configuration suffers from the drawback that, in the absence of a suitable supply scheme and controller, it cannot guarantee any balancing between cells exhibiting differences in terms of e.g. initial SOC, physical parameters, etc. In fact, when an input is applied to the series, the same branch current $I(t)$ flows through all the cells.

The section is organized as follows. The optimization problem, to be solved at each time instant, is introduced in Subsection III-A. The variables to be optimized are the fraction of current which flows through each one of the cells and the branch current itself. The underlying assumption is that it is possible to control each of the former through a suitable supply scheme. In particular, Subsection III-B provides an efficient, easily implementable and reliable design of the latter. Note that this is just one possible implementation which could be replaced by any other feasible one. The method explained below relies on the model described in Section II.

III-A Optimization Problem

In the following we assume that the system is controlled by applying piece-wise constant input actions at discrete times $k$ with a sampling time $T_{s}$. Define as $\mathbf{u}(k)\in\mathbb{R}^{N+1}$

[TABLE]

the set of controllable inputs for the series-connected cells at each $k$, where $\delta_{i}(k)\in\left[0,1\right]$ is the fraction of the branch current $I(k)$ flowing through the $i$-th cell during the time interval $[k,k+1)$. Note that, $I(k)$ is also considered as an optimization variable. In fact, if an adjustable current generator is available, the control of the branch current magnitude over time allows for better performance.

III-A1 Cost Function

In order to consider the control objectives as a scalar index, the following weighted cost function is considered

[TABLE]

with

[TABLE]

where $k_{0}$ is the time instant at which the optimization is performed, $H$ is the prediction horizon, $T_{norm}$ is a normalization temperature which can correspond to the external environment temperature, $I_{max}$ the maximum applicable branch current and $\overline{z}$ the target SOC (usually, $\overline{z}=1$) which the control algorithm aims to reach at the end of the charge. The coefficients $\alpha$ need to be chosen in order to guarantee an optimal trade-off between the different objectives. In particular, the first term (32a) accounts for the objective of charging the cells, i.e. reaching the target state of charge $\overline{z}$. The term (32b) aims to keep the temperature of each cell as low as possible, so as to limit safety risks as well as degradation effects due to ageing phenomena (17). Finally, the control action can be penalized through costs (32c) and (32d), which account for the branch current and its fractions flowing through each of the cells respectively. By limiting the current magnitude, energy losses and heat generation are minimized thus improving overall battery health and efficiency. Costs (32e) and (32f) account for possible penalization of the input ramping rate. The final cost (32g) is introduced with the aim of speeding up the charging process. In fact, it has been noticed in simulation that without such penalty the algorithm often tends to devote most of the charging capability to a subset of the cells in the initial phase, while letting the others at a lower SOC level, so that at the end the total charging time results increased. This could be explained by the fact that the prediction horizon is limited and therefore it is not possible to forecast the long-term future behaviour of the entire system. To this end, the terminal penalty is introduced to consider the fact that additional charging time will be required if differences in the SOCs increase at the end of a certain prediction horizon.

Remark 2

Notice that, since the proposed balancing algorithm is designed with charging capabilities only (i.e. no redistribution is allowed), the reference SOC $\overline{z}$ must be greater than all the cells initial SOC

[TABLE]

III-A2 Constraints

The proper functioning of a charging system requires the consideration of some crucial physical and safety constraints. In particular, the voltage must be limited such that

[TABLE]

III-A3 Problem

The resulting optimization problem, to be solved at each time $k_{0}$ with a prediction horizon of $H$ steps, is then the following

Optimization Problem 1

Find the optimal input sequence $\mathbf{U}^{*}(k_{0})=\left[\mathbf{u}^{*}(k_{0})\cdots\mathbf{u}^{*}(k_{0}+H)\right]$ such that

[TABLE]

with $J$ as in (31), subject to constraints (34).

According to the receding horizon principle, only the first element $\mathbf{u}^{*}(k_{0})$ of the obtained sequence is applied to the system. Then, at the successive time step, the values of the relevant model states are updated according to the obtained measurements and a new optimization is performed.

III-B Practical Implementation

The formulation of Problem 1 is general in the sense that it is designed independently of the underlying power supply circuitry. In this work, with the aim of proposing a possible practical implementation, the following supply scheme is considered.

III-B1 Power Supply System

The elementary unit of the proposed supply system is a PU consisting of the $i$-th Li-Ion cell together with two power switches $S_{i}$ and $\bar{S}_{i}$, as shown in Fig. 7. The two switches work in a complementary way, i.e. when $S_{i}$ is closed $\bar{S}_{i}$ is open and vice versa. Such design allows to fully shunt the cell when necessary. Note that, since the conduction resistance of the switches is assumed to be negligible, the terminal voltage of the $i$-th PU happens to coincide with that of the corresponding cell $V_{i}(t)$. The scheme described above exhibits low cost of implementation, high efficiency and the ability to be easily modularized [10]. On the other hand, a possible drawback is the decay in efficiency when high currents or large number of cells are considered. In these cases, alternative schemes can be adopted [10], which still remain applicable in combination with the presented balancing-aware charging procedure.

According to Fig. 8 the variable input generator provides the branch current $I(t)$ to the whole PUs series. Notice that in such formulation it appears convenient to considered the control variables $\delta_{i}(k)$ as duty cycles, i.e. fractions of the $k$-th time step in which the branch current $I(k)$ flows through the $i$-th cell. The interval $[kT_{s},\,\,(k+1)T_{s})$ is ideally divided into two parts $[kT_{s},\,\,kT_{s}+(1-\delta_{i}(k))T_{s})$ and $[kT_{s}+(1-\delta_{i}(k))T_{s},\,\,(k+1)T_{s})$ as described in Fig. 9, where $I_{app,i}(t)$ is the current actually applied to the $i$-th cell over time. In the first part of the interval, the cell is completely bypassed (i.e. $\bar{S}_{i}$ is open and $S_{i}$ closed in Fig. 7), while in the second part the vice versa occurs so that the cell is charged with a current corresponding to $I(k)$.

Remark 3

Note that to avoid short circuit conditions, one must in practice make the branch current $I(t)$ vanish when all the shunting switches are closed, i.e.

[TABLE]

III-B2 Practical Considerations

The Optimization Problem 1 is designed considering the inputs as constant during the sampling interval. However, the profile of the applied current (Fig. 9) is different from that of a constant current, which in this case can only correspond to the average one $I_{avg,i}(k)=\delta_{i}(k)I(k)$. This brings to wrong predictions of the future states by the MPC algorithm, thus harming optimality and constraints satisfaction if the average value is considered, especially for long prediction horizons. Due to the switched structure of the supply scheme, the problem could be naturally formulated as a mixed integer program [56]. Nevertheless, the online solution of this latter is impracticable as it becomes computationally prohibitive even if only few cells are considered. In order to preserve the possibility of an online implementation, in this work we address the adaptation of the proposed optimization problem to the considered supply circuitry. Since discontinuities in the input are not easily taken into account by optimization procedures, a viable solution is here proposed in order to approximate with sufficient and arbitrary accuracy the step current profile in Fig. 9. In particular, a continuous sigmoid function (see Fig. 10) is considered for any possible control step $k$, such that

[TABLE]

The parameter $a$ in (36) determines the slope of the sigmoid function in $kT_{s}+1-\delta_{i}$ and therefore the accuracy of the approximation, at the expense of a possibly increased stiffness of the resulting ODEs. Simulation results, carried out for an horizon of $200s$ with the input in Fig. 11 are shown in Fig. 12, 13 and 14, where it appears evident that using the average value $I_{avg,i}(k)$ produces non negligible errors, while increasing the parameter $a$ in the sigmoid step approximation (36) enhances considerably the accuracy. Only the evolution of temperature, capacity loss and SEI resistance growth are reported, since they are the unique quantities which exhibit a nonlinear behaviour with respect to the applied current. On the other hand, no benefits in the SOC description can be achieved through the proposed approximation, since its dynamics is linear and sufficiently slow. Notice that the results refer to a simulation time much longer than that expected for an MPC prediction in which the real values of the considered quantities are updated every $T_{s}$ with the measured ones.

Remark 4

The possibility to control the branch current $I(k)$ through a variable current generator allows for better results. In fact, in view of the highly nonlinear relationship between temperature and input current, both the magnitude and the duration of the current applied in the second part of the control step influence the resulting behaviour.

In the following, the possibility of a soft constraints formulation is considered in order to enhance the practical applicability of the presented equalization strategy.

III-B3 Constraints Softening

As for the MPC formulation, a real world application could require the softening of the constraints in (34). In particular, apart from (34c) and (34d), which can not be relaxed outside the ranges of definition, and (34e), which only regards the input, all the other constraints can be relaxed inserting slack variables $s$ such that

[TABLE]

The total cost function $\bar{J}(k_{0})$, at each time instant $k_{0}$, must be then reformulated accordingly as

[TABLE]

with $J(k)$ as in (31) and

[TABLE]

Remark 5

Note that softening the maximum power constraint in (37c) is feasible in practice since generators can deliver a power slightly higher than the nominal one, at least for short periods of time.

III-B4 Practical Problem Formulation

The optimization problem to be solved at each step $k_{0}$ with a prediction horizon of $H$ steps in presence of the power supply scheme described in III-B and the practical considerations therein discussed, becomes then

Optimization Problem 2

Find the optimal input sequence $\mathbf{U}^{*}(k_{0})=\left[\mathbf{u}^{*}(k_{0})\cdots\mathbf{u}^{*}(k_{0}+H)\right]$ such that

[TABLE]

with $\bar{J}$ as in (38), subject to constraints (34), softened as in (37), and (35).

IV Simulation Results

In this section, the effectiveness of the proposed approach is tested in simulation. In order to obtain a realistic scenario, the proposed model is used for the control, while the more accurate LIONSIMBA simulator is considered as the real battery. In IV-A, the (virtual) testbed is described in terms of adopted cells parameters and spatial configuration. At first, a standard charging protocol (namely, the CC-CV) is applied and the results are commented in IV-B. Subsequently, a simple method for equalization at the end of the charge is tested in IV-C. At last, the same battery is charged with the proposed NMPC strategy described in Section III. The improvements obtained with such an algorithm with respect to the previous ones are highlighted in IV-D, where a thorough discussion of the results is carried out. In particular, it is shown that the extractable capacity during the successive discharge cycle is maximized and the safety is guaranteed by constraints satisfaction throughout the charge. It is also evidenced that the ageing effects are reduced.

IV-A Testbed Description

In the simulated battery pack, the considered cells are all Kokam SLPB 75106100. All the parameters, except for the ones related to the ageing dynamics, which are taken from [35], are those experimentally identified in [48, 49]. Notice that this constitutes an appreciable feature, since in this way real commercial cells are accurately modelled. In this work the physical and thermal parameters of the cells are supposed to be equal to the nominal ones given in [48, 49], with the exception of the initial capacity and SEI resistance. In fact, it is known in the literature that these latter exhibit significant variations due to both the manufacturing process and the different ageing exposure. In order to take into account such possible mismatches, the two mentioned parameters are randomly extracted as reported in Table I, where also the initial normalized SOC is listed for each cell.

The topology of the cells series is that depicted in Fig. 15, with two packs of three adjacent cells divided by a space in which the coolant is able to flow (light gray area). It also wraps the rest of the two packs, so that the most external cells have a far higher heat exchange surface compared to the inner ones. A representation of the thermal couplings, described by means of an equivalent electric circuit, is reported in Fig. 16 along with the adopted related values.

Finally, the positive and negative open circuit potentials $\bar{U}_{p}(t)$ and $\bar{U}_{n}(t)$ for the considered cells are given by the following functions of the surface stoichiometries $\theta_{p}(t)$ and $\theta_{n}(t)$ respectively

[TABLE]

Since such functions must be referred to the particular type of cell under investigation, they are obtained by a fitting procedure on the experimental data collected in [48]. The resulting relationships are provided in Fig. 17 for the readers convenience. Finally, the temperature-dependent electrolyte conductivity function is expressed as follows

[TABLE]

where $\gamma_{j}^{[k]}(t)=10^{-3}c_{e,j}^{[k]}(t)$. Without loss of generality, in this work, the thermal power (in [$W$]) that the cooling system releases in the external environment is considered constant in the operating range, for the sake of simplicity. In particular,

[TABLE]

IV-B Standard Charging Method

One of the most used charging protocols in industry is the well-known CC-CV [57, 58]. When this protocol is applied to a set of series-connected cells, these latter are charged while keeping the total branch voltage within a specified threshold. Such limit is computed as the voltage corresponding to SOC=100% for a nominal cell (in this case, taken as $4.15V$) multiplied by the number of cells in the series. In particular, at a first stage a constant current is applied until the limit voltage is reached and, then, such voltage is kept constant reducing the current until a specified threshold is attained. Due to the fact that only the branch voltage is considered, some of the cells usually result in being undercharged and others overcharged. For these latter safety is at risk, while the capacity of the former is obviously not fully exploited. In order to avoid overvoltage exposure, using a voltage sensor for each of the cells, it is possible to interrupt the charging procedure as soon as one of the cells reaches its upper voltage threshold [59]. This, however, reduces even more the overall used battery capacity.

The CC-CV protocol has been applied to the testbed described in IV-A. Fig. 18(a) reports the input branch current applied during the CC-CV. A constant current corresponding to a nominal $I_{1C}$ is firstly imposed until the total voltage of the series reaches $24.9\,V$ (Fig. 18(b)). Due to the differences in terms of total capacity, forcing the same current through all the cells may result in the attainment of different value of SOCs at the end of the charging procedure. This could even result in the amplification of the cell differences over time, as can be seen for the SOC in Fig. 19. Moreover, the total voltage threshold is hit when some cells are well beyond their maximum limit, while others are far below full charge (Fig. 20). While this can constitute a risk for safety (it can even lead to fire or explosions [60]) and a sure harm to the battery health, it is also a waste of capacity, since, as made evident by Fig. 19, the majority of the cells are not even close to their maximum state of charge. The temperature of the various cells, in the first part of the charge, grows mainly according to their characteristics and physical disposition (Fig. 21). In the second part of the charging, the current is decreased in order to maintain the total voltage constant. At the end of the constant-voltage procedure, SOC unbalance is still present, and therefore the total cell capacity is not exploited in a satisfactory way. During this second phase, as expected, the temperatures decrease in response to the lower applied current. The threshold shown in Fig. 21 is drawn for the sake of comparison (it will be used as a constraint in the proposed NMPC algorithm).

IV-C Voltage-Based Method

In order to obtain better performance than the CC-CV protocol, we consider also as a benchmark a simple model-less charging method which exploits the supply circuit described in Section III. Such a strategy only relies on voltage measurements and applies a constant branch current to all the cells until one of them hits its voltage threshold. Then, this latter is fully bypassed, while the charging continues for those with a voltage below the respective limit. Notice that in this case, since the constant-current phase is not followed by a constant-voltage one (a different voltage source for each cell would be required for performing the CV phase on each element of the series independently), we chose as voltage threshold the maximum allowed voltage rather than the one corresponding to $100\%$ of state of charge.

In the carried out simulation, a constant current of $I_{1C}$ has been applied, giving rise to the results reported in Fig. 22, 23 and24. As expected, the charging of each cell is interrupted as soon as the voltage threshold is hit (the maximum voltage $4.2V$, given in [48], is considered). Therefore, the voltage limits are always respected by design. With respect to the CC-CV, it can be noticed that the SOCs are more balanced. However, they do not reach the target value of $100\%$ and even the cell with the lowest capacity is not fully exploited.

IV-D Balancing-Aware NMPC Method

The algorithm proposed in this Section III is now applied to the same battery pack considered above. The values of the parameters adopted in the optimization are listed in Table II. Notice that, due to the operating features of the considered supply scheme, the penalty on the variations of the duty cycles $\delta_{i}(k)$ and the branch current $I(k)$ loses its practical meaning, so that $\alpha_{5}=\alpha_{6}=0$ shall be used. However, in the presented simulation they are taken not null with the unique purpose of enhancing the readability of the plots (for sufficiently low values of such coefficients the overall performance is not affected). The slack variables weights for the soft constraints are taken so high that practically no relaxation takes place, i.e. $\alpha_{s_{1}}=\alpha_{s_{2}}=\alpha_{s_{3}}=10^{15}$.

In the following the obtained results are briefly discussed. In particular, Fig. 25 reports branch current $I_{branch}(t)$ and duty cycles $\delta_{i}(t)$, which produce the average currents plotted in Fig. 26. During the whole process, the cells are charged with different mean currents, with the objective of reaching for each cell the $100\%$ of SOC while reducing the unbalance. These goals are achieved by a suitable tuning of $\alpha_{1}$ and $\alpha_{7}$, respectively. The inputs are obtained taking also into account the power curve of the generator. This translates into a full power supply exploitation, thus allowing for smaller and cheaper generators. As can be evidenced from Fig. 27, the SOC unbalance is reduced over time and the procedure is stopped when all the cells are fully charged. Fig. 28 and 29 highlight the fact that both the voltage and the temperatures remain safely under the imposed limits (only the cell 2 hits and keeps the temperature threshold near the end of the charge). This helps reducing wear and ageing, as well as improving safety.

A key feature of the NMPC algorithm is that, with respect to the other two tested strategies, it allows for a full exploitation of the cells at the end of the charge, i.e. the SOC of all the cells is the maximum possible according to the limits. As the main benefit of this fact, one has that the charge extractable during the successive discharge cycle is maximum. In fact, assuming to discharge the whole battery pack with a constant current of $I_{1C}$ until a cut-off voltage threshold ($2.7V$) is reached by one of the cells, the charge obtainable after the application of the voltage-based method (IV-C) is $6.005Ah$, while after the NMPC it is as high as $6.36Ah$ ($6\%$ more). The same arguments do not hold when taking into account the CC-CV, since it appears evident that the extractable charge is usually higher than that of the NMPC. However, the comparison is not fair since this comes at the cost of an overcharge of some of the cells (which reach a SOC as high as $124.55\%$) with sure harm to safety.

Also as for the ageing effects, the proposed algorithm outperforms the CC-CV, both from the point of view of capacity loss (30) and resistance growth (31). With respect to the voltage-based algorithm, the gain is less noticeable, but still present.

V Conclusions

In the present work a NMPC strategy has been presented for optimally charging series-connected Lithium-ion cells, while avoiding the necessity of periodic offline balancing. A sufficiently accurate electrochemical model, suitable for control, has been formulated in Section II, which takes into account also ageing and thermal effects (comprising the coolant dynamics). This latter has then been used to establish a general NMPC framework (III-A). Finally, a possible practical implementation has been considered in III-B, in which the formulation is adapted to the case of a specific implementable power supply scheme. The effectiveness of the proposed methodology has been validated on a virtual testbed based on the P2D electrochemical model, in which real cells parameters have been considered in order to obtain a realistic scenario. Simulations have shown the ability of the presented algorithm to rapidly achieve state of charge balancing while guaranteeing safety and battery health-related constraints. The comparison with a standard charging protocol (CC-CV) and a simple voltage-based procedure has highlighted that the proposed approach outperforms the other considered methods from all points of view.

In order to expand the possibilities offered by nonlinear model predictive control in the context of advanced battery management systems, future research work could be devoted to solutions able to cope with parallel-of-series of cells, as well as possibly distributed architectures, in order to keep at a reasonable level the computational power requirements also in presence of thousands of cells.