Optimal Charging of an Electric Vehicle Battery Pack: A Real-Time Sensitivity-Based MPC approach
Andrea Pozzi Marcello Torchio Richard D. Braatz Davide M. Raimondo


















[mycorrespondingauthor]Corresponding author
andrea.pozzi03@universitadipavia.it
University of Pavia, Via Ferrata 5, 27100, Pavia

United Technologies Research Centre Ireland Ltd, 2nd Floor Penrose Business Centre, Cork, Ireland

Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA, 02139

University of Pavia, Via Ferrata 5, 27100, Pavia














Lithium-ion battery packs are usually composed of hundreds of cells arranged in series and parallel connections. The proper functioning of these complex devices requires suitable Battery Management Systems (BMSs). Advanced BMSs rely on mathematical models to assure safety and high performance. While many approaches have been proposed for the management of single cells, the control of multiple cells has been less investigated and usually relies on simplified models such as equivalent circuit models. This paper addresses the management of a battery pack in which  each cell is explicitly modelled as the Single Particle Model with electrolyte and thermal dynamics. A nonlinear Model Predictive Control (MPC) is presented for optimally charging the battery pack while taking voltage and temperature limits on each cell into account. Since the computational cost of nonlinear MPC grows significantly with the complexity of the underlying model, a sensitivity-based MPC (sMPC) is proposed, in which the model adopted is obtained by linearizing the dynamics along a nominal trajectory that is updated over time.
The resulting sMPC optimizations are quadratic programs which can be solved in real-time even for large battery packs (e.g. fully electric motorbike with 156 cells) while achieving the same performance of the nonlinear MPC. 













Lithium-ion batteries, Battery management systems, Advanced battery management systems, Model predictive control, Predictive control









§ INTRODUCTION


Electric vehicle battery packs usually consist of several cells which are arranged in series and parallel connections in order to meet power and capacity requirements. The proper functioning of these complex devices requires Battery Management Systems (BMSs) <cit.>. One of the main tasks of a BMS is to safely charge the battery. This objective is usually addressed through standard protocols, such as Constant Current or Constant Current–Constant Voltage  <cit.>. These methods rely on fixed voltage limits that assure reasonable performance for each cell composing the pack over its lifetime.
In practice, these approaches lead to a suboptimal exploitation of the battery, in which some cells are undercharged and some are overcharged. As a consequence, some cells may experience fast degradation, thermal runaway, and, in certain cases, even explosions.  
These problems can be significantly alleviated if advanced BMSs, which rely on mathematical models of the pack, are employed <cit.>. 


In this work, we focus on batteries based on lithium-ion chemistries which, thanks to their unique characteristics, have proven to be the most promising energy accumulators for many applications <cit.>. Within this context, the two main categories of cell models employed in advanced BMSs are: Equivalent Circuit Models (ECMs, <cit.>) and Electrochemical Models (EMs, <cit.>). While the former are simple and intuitive, the latter provide a detailed description of the electrochemical phenomena which occur inside a cell. Among the EMs, the Pseudo-Two-Dimensional (P2D) model <cit.> – also known as Doyle-Fuller-Newman – is the most widely used. This model consists of coupled and nonlinear Partial Differential Algebraic Equations (PDAEs).  Due to its high computational cost, the P2D is more suited for simulation purposes rather than for on-board control applications. Moreover, the use of the P2D model within a control framework is limited by identifiability and observability issues <cit.>.  
For all these reasons, the research community has been interested in the development of simplified electrochemical models, which are faster to simulate, are identifiable and observable, and still provide a reasonable description of the internal cell phenomena <cit.>. Among them, the Single Particle Model (SPM) <cit.>), which is derived from the P2D by modelling the electrodes as single particles, has received a lot of attention. The parameters in the SPM are fit to input-output data from the battery. The parameter identifiability and state observability of the SPM have been analyzed in several works, e.g. <cit.>. While the SPM has much lower computational cost than the P2D, which is an advantage for control purposes, the SPM neglects electrolyte, thermal, and ageing dynamics.
The model fidelity has been increased by extending the SPM to include electrolyte dynamics (SPMe) <cit.> and thermal dynamics (SPMeT) <cit.>.

The model-based charging of lithium-ion batteries has been addressed by many authors over the years. Most of the literature focuses on the control of a single lithium-ion cell. Within this context, different control approaches have been considered, such as fuzzy logic <cit.>, empirical rules <cit.>, and optimization-based strategies <cit.>. 

Model Predictive Control (MPC) <cit.> appears to be the most used optimization-based methodology for charging of cells <cit.>. MPC is particularly appropriate for controlling multivariable nonlinear systems while taking an objective function and constraints on both inputs and states into account. In the context of lithium-ion cells, MPC aims to minimize the charging time while satisfying temperature and voltage constraints.

In particular, the works in <cit.> have proposed  MPC strategies based on ECMs, while <cit.> have suggested the use of electrochemical models together with MPC in order to achieve better performance. In order to alleviate the complexity of physics-based models, input-output descriptions have been used in combination with MPC in <cit.> and <cit.>. In particular,  <cit.> relied on a linear input-output model, while a piecewise Auto Regressive Exogenous (ARX) input-output model was proposed in <cit.>. Finally, in order to achieve a reasonable tradeoff between computational time and accuracy in the prediction, the usage of Linear Time Varying (LTV) models for MPC has been addressed <cit.>. The works listed above focused on the charging of a single lithium-ion cell. 

Batteries composed of several cells are necessary in many applications, such as hybrid electric vehicles. The optimal charging of battery packs has been less investigated than the case of single cells. In particular, most of the available literature relies, as control models, on very simple lumped ECMs (see e.g. <cit.>). Few works tackle the optimal control of lithium-ion batteries by directly modelling each cell individually. This level of detail is necessary in important tasks such as the model-based State of Charge (SOC) balancing of series-connected cells <cit.>. Most of the research produced in this area relies on linear ECMs for each cell (with the exception of <cit.> which is based on the SPMeT). The usage of such models in the battery optimal control context has the advantage of a low computational cost, which allows for real-time implementation even in the case of a high number of cells. On the other side, simple linear models  often fail to grasp the real behaviour of the battery pack. As a consequence, the resulting charging may be suboptimal and may not even satisfy the safety constraints on all the different cells.
The use of electrochemical models within a control framework would be a possibility, but comes at the price of a prohibitive computational cost.
As a compromise, the use of linearized  electrochemical models seems promising in order to achieve high performance with a reasonable computational load.

In this work, we consider a battery pack composed by series-connected modules, each of which are constituted of parallel-connected cells modelled according to the SPMeT. The resulting model consists of a set of nonlinear Differential Algebraic Equations (DAEs), for which suitable linearization methods are applied. In the following, in order to provide high accuracy for the model used for the  control, we rely on a sensitivity-based linearization method, which has been described in detail in <cit.> for a system of Ordinary Differential Equations (ODE) and here extended to the case of a system of DAEs. Such an approach differs from standard LTV techniques <cit.>, since the sensitivities of the states and outputs to input variations are continuously integrated together with the model equations rather than evaluated only at discrete time steps.

The main contribution of this work is the application of a sensitivity-based linear MPC  to the management of large battery packs, taking into account safety constraints (such as temperature and voltage limits) on all the different cells. To the best of our knowledge, this is the first time that such an approach is used in the context of the control of battery packs. As shown in the Results section, a detailed analysis is conducted by comparing a standard nonlinear MPC and the Constant-Current Constant-Voltage algorithm with the proposed sensitivity-based MPC over an increasing number of series and parallel connected cells. The sensitivity-based MPC has much lower computational cost than nonlinear MPC while achieving comparable performance. Since the choice of the nominal trajectory is fundamental in order to obtain an accurate linearization, we also provide an adaptive method in order to update the nominal trajectory during the charge. In addition, the need for an optimal management of the different cells is made evident by highlighting the disadvantages of standard charging protocols, such as the CC-CV. The value of the proposed methodology is demonstrated by application to the optimal management of an electric motorbike battery pack composed by 156 cells. In this application, the sensitivity-based MPC provides optimal performance with a computational cost compatible with the sampling time. 

The paper is organized as follows. Section <ref> recalls the equations of the model used in the control algorithm, while Section <ref> proposes the sensitivity-based  model predictive control algorithm. The adapting of this MPC formulation for the context of lithium-ion batteries and the main results are provided in Section <ref>. Finally, Section <ref> concludes the paper. 





















§ MODEL


In this section, the main equations of the simplified electrochemical model proposed in <cit.> are recalled.  The latter consist of a reduction of the P2D <cit.> model which is detailed but also suitable for control purposes. In particular, starting from the SPMe <cit.>, the Partial Differential Equation (PDE) which describes the diffusion of ions within the electrolyte is spatially discretized according to the finite volume method <cit.>. Moreover, a polynomial approximation of the ion concentration along the radial axis of each electrode is considered in order to reduce the Fick's laws into Ordinary Differential Equations (ODEs) <cit.>, as previously done in the context of SPMe by <cit.>. 
In addition, the thermal dynamics are described by adapting the equations proposed in <cit.>, where the heat generated by the cell is transferred to a proper cooling system. The cell model is presented in Section <ref>, while the equations of the battery pack model are given in Section <ref>.



 §.§ Model of a Single Cell
 
In this subsection, the index j ∈{p,s,n} refers to all the cell sections, while the index i ∈{p,n} is used in equations valid only for the electrodes. The independent variables x∈ℝ and r ∈ℝ are the axial and radial coordinates respectively, and t ∈ℝ is the time. 
According to the approximation in <cit.>, the ion concentration along the radial axis r of each electrode is described by a fourth-order polynomial function of r. This approach results in coefficients that are functions of the solid average concentration c̅_s,i(t) and the average concentration flux q̅_i(t). In the following, we introduce the dynamics of such variables that are necessary for reconstructing the value of the ion concentration along the radial axis. Consider the average stoichiometry in the electrodes defined by

    θ̅_i (t)=c̅_s,i(t)/c_s,i^max

where c_s,i^max is the maximum solid concentration. 
Relying on the fact that the moles of lithium in the solid phase are preserved <cit.>, the average stoichiometry in the anode can be expressed in terms of the cathode by

    θ̅_n(t)= θ^0%_n+θ̅_p(t)-θ^0%_p/θ_p^100%-θ^0%_p (θ_n^100%-θ^0%_n )

where θ_i^0% and θ_i^100% represent the value of the stoichiometries  respectively when the cell is fully discharged and completely charged. 
From <cit.>, the temporal evolution of the average stoichiometry in the positive electrode can be approximated by

    θ̇̅̇_p(t)=3 I_app(t)/a_p R_p,p L_p F A

where  I_app(t) is the applied current input (defined as being negative during charging by convention), R_p,i is the particle radius, L_j is the thickness of the jth section, F is the Faraday constant, A is the contact area between solid and electrolyte phase, and a_i=3 ϵ_i^act/R_p,i is the specific active surface area, with the active material volume fraction ϵ_i^act defined by


    ϵ^act_p   =-C/Δθ_p A F L_p c_s,p^max
    ϵ^act_n   =C/Δθ_n A F L_n c_s,n^max


in which  C is the cell capacity and Δθ_i= θ_i^100%- θ_i^0%.
In accordance with <cit.>, the volume-averaged concentration fluxes can be described by 


    q̇̅̇_p(t)    = -30 D_s,p(T(t))/R_p,p^2q̅_p(t) + 45/2R_p,p^2 F A L_p a_pI_app(t)
    q̇̅̇_n(t)    = -30 D_s,n(T(t))/R_p,n^2q̅_n(t) - 45/2R_p,n^2 F A L_n a_nI_app(t)
 

where 

    D_s,i(T(t))=D_s,i^0 e^-E_a,D_s,i/RT(t)

is the solid diffusion coefficient for the ith section which depends on the cell temperature T(t) according to the Arrenhius law,
D_s,i^0 is a pre-exponential coefficient which is assumed to be constant, E_a,D_s,i is the activation energy associated with the parameter D_s,i(T(t)), and R is the universal gas constant. 
Then, in accordance with <cit.>, the surface stoichiometries in the positive and negative electrodes are respectively given by the algebraic equations


    θ_p(t)   = θ̅_p(t)+8R_p,pq̅_p(t)/35c_s,p^max+R_p,p I_app(t)/35D_s,p(T(t))F A L_p a_p c_s,p^max
    θ_n(t)   = θ̅_n(t)+8R_p,nq̅_n(t)/35c_s,n^max-R_p,n I_app(t)/35D_s,n(T(t))F A L_n a_n c_s,n^max


Finally, the  state of charge  SOC(t) is defined as

    SOC(t)=100θ̅_n(t)-θ_n^0%/θ_n^100%-θ_n^0%


In this work, the PDEs governing the diffusion of the electrolyte concentration c_e,j(x,t)  <cit.> are discretized according to the FV method, as previously done in this context by the authors in <cit.>, where the spatial domain is divided into P non-overlapping volumes for each section. The kth volume, with k=1, ⋯, P, of the jth section is centered at the spatial coordinate x_j,k and spans the interval Ω_j,k=[x_j,k̅,x_j,k], whose  width is  Δ x_j=L_j/P.   Defining c_e,j^[k](t) as the average electrolyte concentration over the kth volume of jth section gives 



    ϵ_p ∂  c^[k]_e,p(t)/∂ t=   .[D̃_e(x,T(t))/Δ x_p∂ c_e,p(x,t)/∂ x]|_x_p,k^x_p,k̅-1-t_+/FAL_pI_app(t)
    ϵ_s ∂  c^[k]_e,s(t)/∂ t=   .[D̃_e(x,T(t))/Δ x_s∂ c_e,s(x,t)/∂ x]|_x_s,k^x_s,k̅
    ϵ_n ∂  c^[k]_e,n(t)/∂ t=   .[D̃_e(x,T(t))/Δ x_n∂ c_e,n(x,t)/∂ x]|_x_n,k^x_n,k̅ +1-t_+/FAL_nI_app(t)


where t_+ is the transference number, ϵ_j is the material porosity, and the other terms are evaluated as explained in detail in <cit.>. In particular, the electrolyte diffusion coefficients D̃_e(x,T(t)) are computed from

    D̃_e(x,T(t))=
     D_e,1   if   x ∈{x_p,P,x_s,1}
    
    D_e,2     if   x ∈{x_s,P,x_n,1}
    
    D^eff_e,j(T(t))    otherwise

with


    D_e,1   =ℋ(D^eff_e,p(T(t)),D^eff_e,p(T(t)),Δ x_p,Δ x_s)
     
    D_e,2   =ℋ(D^eff_e,s(T(t)),D^eff_e,n(T(t)),Δ x_s,Δ x_n)


where ℋ is the harmonic mean operator, defined as

    ℋ(ρ_1,ρ_2,λ_1,λ_2)=ρ_1 ρ_2 (λ_1 + λ_2)/ρ_1λ_2+ρ_2λ_1
 
and D_e,j^eff(T(t)) = D_e(T(t))ϵ_j^p_j with  p_j being the Bruggeman coefficient and D_e(T(t)) being the diffusion coefficient within the electrolyte which depends on the temperature according to the Arrenhius law as in (<ref>).

The terminal voltage is given by

    V(t)    =  -I_app(t)R_sei+U̅_p(t) - U̅_n(t) + η̅_p(t) - η̅_n(t) +  ΔΦ_e(t)

where R_sei(t) is the SEI resistance, the Open Circuit Potentials (OCPs) in the positive and negative electrodes are given by


    U̅_p(t)=   18.45θ_p^6(t)-40.7 θ_p^5(t)+20.94θ_p^4(t)
       +8.07θ_p^3(t)-7.837θ_p^2(t) + 0.02414θ_p^1(t)+4.571
    
    U̅_n(t)=   0.1261θ_n(t)+0.00694/θ_n^2(t)+0.6995θ_n(t)+0.00405


whose expressions in terms of surface stoichiometries are obtained by fitting the experimental data collected in <cit.> and depend on the considered cell (in this case, Kokam SLPB 75106100), and η̅_p(t) and η̅_n(t)  are the overpotentials for the positive and negative electrodese given by


    η̅_p(t)    = 2RT(t)/Fsinh^-1(-I_app(t)/2 A L_p a_p i̅_0,p(t))  
    η̅_n(t)    = 2RT(t)/Fsinh^-1(I_app(t)/2 A L_n a_n i̅_0,n(t))


with the exchange current density defined as

    i̅_0,i(t) = Fk_i(T(t))√(c̅_e,i(t)θ_i(t)(1-θ_i(t)))

where k_i(T(t)) is the temperature-dependent (according to the Arrenhius law) rate reaction constant and c̅_e,i(t) is the average electrolyte concentration in the ith section, approximated by

    c̅_e,i(t)=1/P∑_k=1^P c_e,i^[k](t).
 
Moreover, ΔΦ_e(t) is computed as 

    ΔΦ_e(t) = Φ_e^drop(t) + 2RT(t)/F(1-t_+)log_e(c_e,p^[1]/c_e,n^[P])

where, assuming that the ionic current i_e(x,t) has a trapezoidal shape over the spatial domain <cit.>, the electrolyte voltage drop Φ_e^drop(t) can be  approximated by

    Φ_e^drop(t)≃  -I_app(t)/2 P( ϕ_p(t) + 2ϕ_s(t) + ϕ_n(t) )

in which



    ϕ_p(t)   =Δ x_p ∑_k=1^P 2k-1/κ(c_e,p^[k](t),T(t))ϵ_p^p_p
    ϕ_s(t)   =Δ x_s ∑_k=1^P1/κ(c_e,s^[k](t),T(t))ϵ_s^p_s
    ϕ_n(t)   =Δ x_n ∑_k=1^P2P-2k+1/κ(c_e,n^[k](t),T(t))ϵ_n^p_n


where κ(c_e,j^[k](t),T(t)) is the temperature-dependent electrolyte conductivity for the kth volume of the jth section, which is usually expressed as an empirically derived nonlinear function of the electrolyte concentration in that volume:

    k(γ_j^[k](t),T(t))=   (0.2667 (γ_j^[k](t))^3 -1.2983  (γ_j^[k](t))^2
        +1.7919 γ_j^[k](t) + 0.1726)e^-E_a,κ/RT(t)

where γ_j^[k](t)=10^-3c_e,j^[k](t). The expression in (<ref>) is still referred to the case of  Kokam SLPB 75106100. 



As mentioned above, this paper considers a lumped thermal model <cit.> in which the temperature dynamics are given by

    C_thṪ(t)= Q(t) - T(t)-T_sink/R_th

where C_th is the  thermal capacity of the  cell and R_th is the thermal resistance between the cell and the coolant, whose temperature is here assumed to be constant and equal to T_sink.  The heat generation term Q(t) is due to the cell polarization and is described by

    Q(t)=| I_app(t) |·| V(t) - (U̅_p(t)-U̅_n(t)) | .

In order to simplify the notation in the next sections, define the variables for the average stoichiometry and concentration flux


    θ(t)   =θ_p(t)
    q(t)   =[q_p(t), q_n(t)]^⊤

and express the vector of the electrolyte concentrations in the different finite volumes as

    c^e(t)=   [c_e,p^[1](t), ⋯,  c_e,p^[P](t), c_e,s^[1](t), ⋯,  c_e,s^[P](t), c_e,n^[1](t), ⋯,  c_e,n^[P](t)]^⊤.





 §.§ Model of a Battery Pack

This subsection considers the battery pack configuration in Figure <ref>, which consists of N series-connected modules, each of which is constituted of M parallel-connected cells. The total number of cells is given by N_cells=N M, with all cells modelled according to the equations in Section <ref>. The modelling of a lithium-ion battery pack as a set of connected electrochemical models constitutes a key point of this work with respect to the existing literature, in which the battery is usually described with a very simple equivalent circuit model. The use of a more accurate model is motivated by the fact that lithium-ion battery packs are complex systems which require suitable control strategies in order to guarantee the satisfaction of safety constraints on each single cell. However, the usage of an accurate model in optimal control has the disadvantage of high computational burden. This latter can be reduced by properly linearizing  the dynamics (see Section <ref>).

In the following, the indexes i=1, ⋯, N  and j=1, ⋯, M which refer to the jth cell of the ith module. Moreover, the voltage, state of charge, temperature, and applied current of the cells are indicated respectively by V_i,j(t), SOC_i,j(t), T_i,j(t), and I_i,j(t), while the voltage of the ith module is defined as V_i(t). 

The scheme in Figure <ref> consists also of a supply circuitry which allows   charging of the battery (a battery charger) and balancing the energy stored in the different cells  (current generators). In particular, the battery charger supplies the current I_ch(t)≥0, which can be assumed constant (I_ch(t)=I_ch), while the N generators I_b,i(t)≥0,  i=1, 2, ⋯, N allow draining of the current which flows through the different modules. In the case of constant I_ch, these latter represent the system inputs and are fundamental in order to achieve an optimal current control of the different cells. Such current generators can be realized in practice in several ways. For instance, simple proportional-integral-derivative (PID) controllers can be used to regulate the value of a variable resistor in parallel to each module. As an alternative, a Pulse Width Modulation (PWM) approach or a method based on the cell bypass through active elements <cit.> can be adopted. In this work, the heat generation due to the bypassing process is assumed to be negligible (see <cit.> for more details on the different supply and balancing circuitry schemes).
The battery pack is modelled as a system of DAEs with N_cells algebraic variables (I_i,j(t), i=1, ⋯, N, j=1, ⋯, M) and corresponding N_cells algebraic equations. Each module i, with i=1, ⋯,N, is described by M algebraic equations:


    V_i,1(t)=V_i,2(t)
       V_i,2(t)=V_i,3(t)
             ⋮
       V_i,M-1(t)=V_i,M(t)
       I_ch(t)=-∑_j=1^M I_i,j(t)+I_b,i


These equations correspond to the Kirchhoff's current law at each module and are written according to the conventions that the supply charger provides a positive current in order to charge the battery pack, the battery cells are charged by negative currents, and the bypassing system reduces the battery charging by draining positive currents.
The current generated by the battery charger is assumed to be completely bypassed through the generators I_b,i when the ith module completes its charging procedure (i.e. I_b,i(t)=I_ch, t≥t̅_i, where t̅_i is the time at which the ith module is completely charged). The heat exchange between the cells is also assumed to be negligible.  







§ SENSITIVITY-BASED MODEL PREDICTIVE CONTROL


This section proposes a sensitivity-based linear MPC (sMPC) for the control of general nonlinear continuous-time systems described by semi-explicit DAEs as in Section <ref>. The main advantage of this approach is the significant computational time reduction with respect to  nonlinear MPC (nMPC, see Section <ref>), while having comparable performance.

The sMPC strategy, which is presented in Section <ref>, relies on the model linearization around a nominal input trajectory, and is based on the computation of the sensitivities of states, outputs, and algebraic variables with respect to input perturbations. Such sensitivities are obtained by integrating additional continuous-time differential equations together with the model in (<ref>). Using this approach, the resulting linearized model provides higher accuracy than obtained with standard LTV approaches, in which the sensitivities are evaluated only at discrete time steps.
 
Next, Section <ref> describes how to use the sensitivity-based linearized model in an MPC framework. 
In order to address model mismatches in a practical implementation, the softening of the output constraints in sMPC and nMPC is described in Section <ref>.




 §.§ Nonlinear DAE system

Consider the continuous-time system of semi-explicit nonlinear DAEs described by


    ẋ(t)   =f(x(t),u(t),z(t))
    
    0   =h(x(t),u(t),z(t)) 
    
    y(t)   =g(x(t),u(t),z(t))


where t ∈ℝ is the time, x(t) ∈ℝ^n is the states vector, u(t)∈ℝ^m is the control input, z(t) ∈ℝ^s is the vector of the algebraic variables, y(t) ∈ℝ^p is the output, f:ℝ^n× m × s→ℝ^n and g:ℝ^n× m × s→ℝ^p are the state and output functions respectively, and h:ℝ^n× m × s→ℝ^m specifies the set of algebraic constraints.
The system of DAEs in (<ref>) is assumed to be index-1 <cit.>. Moreover, a digital controller is assumed to apply a piecewise constant input at the discrete times t_k, k∈ℕ with sample time T_s. Within this context, define the generic input sequence applied in the time interval [t_k,t_k+H], with H∈ℕ,  as 

    𝐮_[t_k,t_k+H]=[u^⊤(t_k), u^⊤(t_k+1), ⋯ ,u^⊤(t_k+H-1)]^⊤

while the corresponding temporal evolution of states, algebraic variables, and outputs is obtained by integrating the equations in (<ref>) over [t_k,t_k+H], with initial condition x(t_k)=x_k, to give


    𝐱_[t_k,t_k+H]   =[x^⊤(t_k), x^⊤(t_k+1), ⋯ ,x^⊤(t_k+H)]^⊤
    𝐳_[t_k,t_k+H]   =[z^⊤(t_k), z^⊤(t_k+1), ⋯ ,z^⊤(t_k+H)]^⊤
    𝐲_[t_k,t_k+H]   =[y^⊤(t_k), y^⊤(t_k+1), ⋯ ,y^⊤(t_k+H)]^⊤







 §.§ Nonlinear MPC
 

This section summarizes the main features of  nonlinear MPC, which is a control technique suitable for multivariable nonlinear systems in the presence of constraints. The formulation considers a continuous-time model as in (<ref>) and a digital controller, which keeps the input constant over each sampling time.
Nonlinear MPC requires the solution of a finite-horizon optimal control problem at each time step  t_k_0, whose solution provides  the optimal control sequence 𝐮^*_[t_k_0,t_k_0+H] over a  prediction horizon of H steps, where 

    𝐮^*_[t_k_0,t_k_0+H]=[u^*^⊤(t_k_0), u^*^⊤(t_k_0+1), ⋯ ,u^*^⊤(t_k_0+H-1)]^⊤.

According to the receding horizon paradigm, only the first element u^*(t_k_0) is applied to the system and the remaining future optimal moves  discarded. The optimization is then repeated at the next time step over a shifted prediction window, with the newly available measurements <cit.>. 
 
The resulting optimization to be solved at each time t_k_0 is described below.

Find the optimal input sequence 𝐮^*_[t_k_0,t_k_0+H] that solves 

    𝐮^*_[t_k_0,t_k_0+H]=_𝐮_[t_k_0,t_k_0+H]J(t_k_0)

for the cost function

    J(t_k_0)=∑_k=k_0^k_0+H‖ y(t_k)-y^ref‖^2_Q+∑_k=k_0^k_0+H-1‖ u(t_k)-u^ref‖^2_R +J_reg(t_k_0)

in which y^ref∈ℝ^p and u^ref∈ℝ^m are the reference vectors for the output and for the input, respectively, and the matrices Q ∈ℝ^p × p and R ∈ℝ^m × m are design parameters which are weights in the MPC cost function, with Q≥0 and R > 0. The additional term J_reg(t_k_0) acts as a small regularization factor in order to avoid abrupt variations and spikes in the input control law and it is given by

    J_reg(t_k_0)=∑_k=k_0^k_0+H-1‖ u(t_k)-u(t_k-1)‖^2_R_reg

where R_reg∈ℝ^m × m.
Note that y(t_k)∈ℝ^p is obtained by evaluating (<ref>) at discrete time instants t_k. 
The optimization is subject to the system dynamic in (<ref>) and the constraints 


    u^lb≤ u(t_k)≤ u^ub,      k=k_0, k_0+1, ⋯, k_0+H-1
    
    y^lb≤ y(t_k)≤ y^ub,      k=k_0, k_0+1, ⋯, k_0+H


where u^lb, u^ub∈ℝ^m and  y^lb, y^ub∈ℝ^p  are the minimum and maximum allowable values of the inputs and outputs, respectively.





 §.§ Sensitivity-Based Linearization Along a Nominal Trajectory
 
The main drawback of the optimal control formulation in Section <ref> is its high computational cost, which comes from the use of a nonlinear model to predict the future system behaviour. This step involves the solution of a  nonlinear optimization at each time step whose complexity can be prohibitive for certain online control applications. 
As an alternative, we propose a sensitivity-based linearization of the system (<ref>), which can be exploited to provide a fast MPC solution. As in the previous sections,  the discussion is carried on by considering a digital controller which applies a  piecewise constant input at the discrete times  t_k.   













Consider a nominal input signal u(t) and the nominal input sequence 𝐮_[t_k,t_k+H]=[u^⊤(t_k), u^⊤(t_k+1), ⋯, u^⊤ (t_k+H-1)]^⊤ over the time window [t_k,t_k+H]. The corresponding  nominal trajectories for states, algebraic variables, and outputs are  𝐱_[t_k,t_k+H],  𝐳_[t_k,t_k+H], and 𝐲_[t_k,t_k+H] (see (<ref>)). Define S_x(t,t_k)=∂ x(t)/∂ u(t_k), S_z(t,t_k)=∂ z(t)/∂ u(t_k), and S_y(t,t_k)=∂ y(t)/∂ u(t_k) which are the sensitivities of the states, algebraic variables, and outputs to a variation in the input u(t_k) with respect to its nominal value u (t_k) at the discrete time t_k. In particular, the matrices S_x(t,t_k)∈ℝ^n× m, S_z(t,t_k)∈ℝ^s× m, and S_y(t,t_k)∈ℝ^p× m are obtained by solving, together with (<ref>), the system of equations


    Ṡ_x(t,t_k)   =F_x(t) S_x(t,t_k) +F_z(t)S_z(t,t_k) + F_u(t)Δ_t_k(t) 
    
    0   =H_x(t)S_x(t,t_k) + H_z(t)S_z(t,t_k) + H_u(t)Δ_t_k(t) 
    
    S_y(t,t_k)   =G_x(t)S_x(t,t_k) + G_z(t)S_z(t,t_k) + G_u(t)Δ_t_k(t)


with initial condition S_x(t_k,t_k)=0_𝐧 × 𝐦 – since we assume that (<ref>) is a causal system – and


    F_ν(t)=∇_ν f(x(t),z(t),u(t)),  ν={x, z, u}
    H_ν(t)=∇_ν h(x(t),z(t),u(t)),  ν={x, z, u}
    G_ν(t)=∇_ν g(x(t),z(t),u(t)),  ν={x, z, u}


where ∇_ν is the Jacobian operator with respect to ν and Δ_t_k(t)=H(t-t_k)-H(t-t_k-T_s), with H(t-t_k) being the unitary Heaviside step function

    H(t-t_k)=
    0,   t<t_k
    
    1,   t≥ t_k









Now consider a modified input sequence ũ_[t_k,t_k+H]=[ũ(t_k), ũ(t_k+1), ⋯, ũ (t_k+H-1)]  which can be obtained from the nominal input sequence by

    ũ_[t_k,t_k+H]=𝐮_[t_k,t_k+H]+δ 𝐮_[t_k,t_k+H]

where δ 𝐮_[t_k,t_k+H]=[δ u^⊤(t_k), δ u^⊤(t_k+1), ⋯, δ u^⊤ (t_k+H-1)]^⊤ is the sequence of the input variations over the time window [t_k,t_k+H]. 
The sensitivity-based approximation of the states, algebraic variables, and output trajectories which correspond to the sequence ũ_[t_k,t_k+H] is given by <cit.> 


    𝐱̂_[t_k,t_k+H]=𝐱_[t_k,t_k+H]+Π^x_[t_k,t_k+H]δ 𝐮_[t_k,t_k+H]
    ẑ_[t_k,t_k+H]=𝐳_[t_k,t_k+H]+Π^z_[t_k,t_k+H]δ 𝐮_[t_k,t_k+H]
    ŷ_[t_k,t_k+H]=𝐲_[t_k,t_k+H]+Π^y_[t_k,t_k+H]δ 𝐮_[t_k,t_k+H]

 
where the matrices Π^ν_[t_k,t_k+H], with ν={x, z, y}, are defined by

    Π^ν_[t_k,t_k+H]=
    [             S_ν(t_k,t_k)                        0                        ⋯                        0;           S_ν(t_k+1,t_k)         S_ν(t_k+1,t_k+1)                        ⋯                        0;           S_ν(t_k+2,t_k)         S_ν(t_k+2,t_k+1)                        ⋯                        0;                        ⋮                        ⋮                        ⋱                        ⋮;       S_ν(t_k+(H-1),t_k)     S_ν(t_k+(H-1),t_k+1)                        ⋯ S_ν(t_k+(H-1),t_k+(H-1));                          ]

in which S_ν, for ν={x, z, y}, is obtained by integrating the continuous-time system in (<ref>) together with (<ref>).
This step constitutes the main difference with respect to a standard LTV approach, for which the model sensitivities are computed only at discrete time instants, resulting in a loss of accuracy.  





 §.§ Sensitivity-Based Linear MPC
 
This section presents a model predictive control approach for the system (<ref>) based on the linearized model  (<ref>). Using the latter, the optimization can be formulated as a Quadratic Program (QP) which significantly reduces the computational cost compared to nonlinear MPC, thus enabling the use of the proposed strategy in on-line control applications. The objective function to be minimized at each time step t_k_0 is

    J_lin(t_k_0)=∑_k=k_0^k_0+H‖ŷ(t_k)-y^ref‖^2_Q+∑_k=k_0^k_0+H-1‖δ u(t_k)+u̅(t_k)-u^ref‖^2_R+J_reg(t_k_0)

and the resulting optimization is formulated below.

Find the optimal  sequence of the input variations  δ 𝐮^*_[t_k_0,t_k_0+H] that solves

    δ 𝐮^*_[t_k_0,t_k_0+H]=_δ 𝐮_[t_k_0,t_k_0+H]J_lin(t_k_0)

for the cost function (<ref>). The optimization is subject to the system dynamics in (<ref>) and the constraints 


    u^lb≤ũ(t_k)≤ u^ub,      k=k_0, k_0+1, ⋯, k_0+H-1
    
    y^lb≤ŷ(t_k)≤ y^ub,      k=k_0, k_0+1, ⋯, k_0+H


  
The optimal control sequence is obtained as

    𝐮^*_[t_k_0,t_k_0+H]=𝐮_[t_k_0,t_k_0+H]+δ 𝐮^*_[t_k_0,t_k_0+H]




















































The performance of the sensitivity-based linearization and the corresponding MPC algorithm are significantly affected by the choice of the nominal input sequence 𝐮_[t_k_0,t_k_0+H]. 
Here we adopt the following choice. Suppose that, at the beginning, the nominal input sequence is not known. Then use the most likely input sequence as the initial guess, since no further information is available. Then, after each iteration of the MPC algorithm, the nominal input sequence is updated as

    𝐮_[t_k_0+1,t_k_0+H+1]=[𝐮^*^⊤_[t_k_0+1,t_k_0+H], 𝐮^*^⊤_[t_k_0+H-1,t_k_0+H]]^⊤

where 𝐮^*_[t_k_0,t_k_0+H] is the optimal solution of the sMPC at the time t_k_0.
In case of highly nonlinear systems, the initialization of the nominal input sequence has to be done carefully to achieve a sufficiently accurate linearized model in the first iteration. 




 §.§ Soft Constraints for Practical Implementation
 
The constraints on the outputs for both the nonlinear (<ref>) and linearized (<ref>) systems can be softened in order to deal with, in a practical implementation, possible model mismatches.  In particular, using a set of slack variables ξ∈ℝ^p, with ξ≥0, the constraint in (<ref>) can be relaxed as


    y^lb≤ y(t_k) +ξ(t_k),      k=k_0, k_0+1, ⋯, k_0+H 
    
    y(t_k)-ξ(t_k)≤ y^ub,      k=k_0, k_0+1, ⋯, k_0+H


and


    y^lb≤ŷ(t_k) +ξ(t_k),      k=k_0, k_0+1, ⋯, k_0+H 
    ŷ(t_k)-ξ(t_k)≤ y^ub,      k=k_0, k_0+1, ⋯, k_0+H


can be used for the constraint in (<ref>).

The cost functions for sMPC and nMPC are reformulated accordingly as 


    J(t_k_0)=   ∑_k=k_0^k_0+H‖  y(t_k)-y^ref‖^2_Q+∑_k=k_0^k_0+H-1‖  u(t_k)-u^ref‖^2_R+J_reg(t_k_0)+J_s(t_k_0)
    
    J_lin(t_k_0)=   ∑_k=k_0^k_0+H‖ŷ(t_k)-y^ref‖^2_Q+∑_k=k_0^k_0+H-1‖δ u(t_k)+u̅(t_k)-u^ref‖^2_R+J_reg(t_k_0)
       +J_s(t_k_0)


with

    J_s(t_k_0)=∑_k=k_0^k_0+H-1c^⊤ξ(t_k)

where c∈ℝ^p is a suitable vector of weights.















§ RESULTS


This section evaluates the proposed methodology to the control of a lithium-ion battery pack. 
Section <ref> introduces the optimal control problem, while the parameters and simulation settings are given in Section <ref>. The simulation results are provided in Sections <ref>–<ref>. Section <ref> compares the proposed sensitivity-based approach to nonlinear MPC, and a standard charging method, namely the Constant Current-Constant Voltage (CC-CV), is considered as a benchmark in Section <ref>. Then, Section <ref> analyzes how the  computational times of sMPC and nMPC grow with  increasing the number of cells. Finally, the methodologies are tested on a challenging scenario of the battery pack of an electric motorbike. 



 §.§ Optimal Control of a Lithium-Ion Battery

In this section, the optimal control methods in Section <ref> are adapted to the optimal management of a lithium-ion battery, whose cells are arranged as in Figure <ref>, with N series modules of M parallel-connected cells. The total number of cells is given by N_cells=NM. The objective of the control algorithm is to bring the state of charge of each cell as close as possible to  100% while satisfying all the safety constraints.
Each cell is modelled according to the equations described in Section <ref>, while the model of the whole battery pack is given in Section <ref>. The variables are defined by


    u(t_k)   =[I_b,1(t_k), I_b,2(t_k), ⋯, I_b,N(t_k)]^⊤
    
    x(t_k)   =[x^⊤_1,1(t_k), x^⊤_1,2(t_k), ⋯, x^⊤_N,M(t_k)]^⊤
    
    z(t_k)   =[I_1,1(t_k), I_1,2(t_k), ⋯, I_N,M(t_k)]^⊤
    
    y(t_k)   =[y^⊤_1,1(t_k), y^⊤_1,2(t_k), ⋯, y^⊤_N,M(t_k)]^⊤


where x_i,j(t_k)∈ℝ^N_x and y_i,j(t_k)∈ℝ^N_y represent respectively the states and outputs of the jth cell of the ith module. The notation introduced in (<ref>) implies that


    x_i,j(k)   =[θ̅_i,j(k), q̅^⊤_i,j(k), c^e^⊤_i,j(k),  T_i,j(k) ]^⊤
    
    y_i,j(k)   =[V_i,j(k), T_i,j(k), I_i,j(k),  SOC_i,j(k)]^⊤


where θ_i,j(k), q̅_i,j(k), and c^e_i,j(k) refer to the stoichiometry, average concentration flux, and electrolyte concentration of the jth cell of the ith module. The elements of the vector c^e_i,j(k) are the values of the electrolyte concentration along the spatial axis which is discretized according to the finite volume method. The resulting model is a semi-explicit continuous time system of DAEs (see (<ref>)) where  n=N_cells  N_x, m=N, s=N_cells, and p=N_cells  N_y, with N_cells=NM and in which N_x=4+3P (for a P number of finite volumes) and N_y=4 are the number of states and outputs of each single cell, respectively. 
The weighting matrices Q ∈ℝ_≥ 0^p× p and R∈ℝ_> 0^m× m were specified to be diagonal






    Q=[ q 0 ⋯ 0; 0 q ⋯ 0; ⋮ ⋮ ⋱ ⋮; 0 0 ⋯ q;   ],
    
    R=[ r 0 ⋯ 0; 0 r ⋯ 0; ⋮ ⋮ ⋱ ⋮; 0 0 ⋯ r;   ]

where  r ∈ℝ_> 0 and q ∈ℝ_≥ 0^N_y × N_y is defined as 

    q=[   q^V     0     0     0;     0   q^T     0     0;     0     0   q^I     0;     0     0     0 q^SOC ].

The reference vector and the limits for the output vector, respectively, are given by 


    y^ref   =[y^⊤_r, y^⊤_r, ⋯,  y^⊤_r]^⊤
    
    
    
    y^lb   =[y^⊤_min, y^⊤_min, ⋯,  y^⊤_min]^⊤
    
    y^ub   =[y^⊤_max, y^⊤_max, ⋯,  y^⊤_max]^⊤


with  y_r, y_min, y_max∈ℝ^N_y defined by


    y_min   =[V_min,  T_min, I_min,  SOC_min]^⊤
    
    y_max   =[V_max,  T_max, I_max,  SOC_max]^⊤
    
    y_r   =[V_r,  T_r, I_r,  SOC_r]^⊤


where I_min and I_max are the maximum and the minimum values of current that can flow through a single cell, V_min and V_max are the upper and lower bounds for the voltage, and T_min and T_max, and SOC_min and SOC_max, are the limits for the temperature and the state of charge respectively. Since the limits on the current flowing through each cell are explicitly considered in the optimization, lower and upper bounds on the input vector (u^lb and u^ub) are not required. Finally, the reference values for voltage, temperature, current, and state of charge are denoted by V_r, T_r, I_r, and SOC_r respectively. 


Although the internal states of each cell are not measurable in practice, all the relevant states are assumed available here. The use of observers goes beyond the scope of this work whose objective is to highlight the suitability of the sensitivity-based approach to real-time optimal control. For the design of observers for battery states, the interested reader can refer to e.g. <cit.>. 




 §.§ Model Parameters and Simulations Settings

Here the virtual testbed considered in the simulations is described in detail in order to allow the presented results to be reproducible by others. The electrochemical parameters adopted for all the cells are those experimentally measured in <cit.> from a complete electrochemical characterization of a commercial cell (the Kokam SLPB 75106100). The thermal capacity, thermal resistance, and sink temperature are assumed equal to C_th=4186J/K, R_th=169.5K/W, and T_sink=298.15K, and the initial value for the temperature is set to T^0=298.15K for all cells. The initial electrolyte concentration and average concentration flux are assumed to start at equilibrium values 1000mol/m^3 and zero respectively. The initial state of charge of the different cells, as well as the capacity and the SEI resistance, are extracted from a Gaussian distribution as


    SOC^0_i,j   ∈𝒩(SOC^0,σ^2_SOC)
    
    C^0_i,j   ∈𝒩(C^0,σ^2_C)
    
    R_sei,i,j^0    ∈𝒩(R_sei^0,σ^2_R_sei)


with SOC^0=50%, 
R_sei^0=15mΩ, 
and C^0=7.5Ah 
(i.e., I_1C=7.5A), 
while the standard deviations are 
σ_SOC=10%, 
σ_R_sei=0.75mΩ, and 
σ_C=0.375Ah.


















The optimization settings are reported in Table <ref>, with the voltage, temperature, and current of each cell appearing in the constraint set but not weighted in the cost function (q^V=0, q^T=0, and q^I=0). The corresponding reference values V_r, T_r, and I_r are set equal to zero. 

All the simulations were performed on a Windows 10 personal computer with 16 Gbytes of RAM and a 2.5 GHz i7vPro processor. The control problems were solved using CasADi <cit.>, a symbolic framework for automatic differentiation. This software offers a Matlab interface for the Interior Point Optimization Method (IPOPT) <cit.> used for solving the optimizations, and for the SUNDIALS suite <cit.> used for integrating the process dynamics. Moreover, CasADi was used for the computation of the sensitivity matrices along the nominal trajectory. In order to provide a fair comparison between sMPC and nMPC, both of the underlying optimizations were solved using IPOPT. 

































 §.§ Comparison Between sMPC and nMPC

This section considers a battery pack composed by N=2 modules with M=2 parallel connected cells for each module. The performance of the proposed sensitivity-based MPC is compared with the nonlinear MPC as the benchmark for the parameters of each cell reported in the previous section. The temporal evolution of the states and outputs obtained by sMPC (dotted line) and nMPC (dashed line) are very similar (see Figures 2a–4a), with nearly complete overlap. For both MPC formulations, the state of charge for all of the cells reach the desired target of 100% within 3500s (Figure <ref>). The constraints on the voltage and temperature for all of the battery cells are satisfied for all time (Figures <ref> and <ref>). The current flowing in the different cells are nearly identical for the two MPC formulations (Figure <ref>),  and the input actions are nearly identical for each module (Figure <ref>), consisting of the drained currents I_b,1 and I_b,2.


Since the main goal of this paper is to develop a control algorithm that achieves high performance while having low enough on-line computational cost being implementable in real time, the mean computational times needed by the two methods to compute the optimal control sequence at each time step are reported in Figure <ref>. While having very similar closed-loop performance, the on-line computational cost of the sensitivity-based MPC is significantly lower than for nMPC, motivating the use of sMPC in the context of real-time control of battery packs. Also, the on-line computational time of nMPC is highly variable while being nearly deterministic for sMPC, as its optimization is a quadratic program whose computational cost for solution is weakly dependent on the values of its parameters. Having low variability in its on-line computational cost is another desirable feature of sMPC.






































































 §.§ Standard CC-CV Method

To demonstrate the need for an optimal management of a lithium-ion battery pack that is able to take into account input and output constraints, a standard charging method – namely the Constant Current-Constant Voltage (CC-CV) – is applied to the same configuration considered in Section <ref>. The CC-CV method is composed of two phases. In the first phase, a constant charging current I_cc is applied to the series connected modules. This phase ends separately for the different modules, as soon as their voltage reaches a predefined  threshold V_th (that in this case is assumed to be equal to 4.15V, which corresponds to the OCPs difference value at 100% of SOC). Note that this value can be lower than the maximum allowed voltage specified by the cell data-sheet. The second phase consists of a separate constant voltage charging of the different modules with a voltage generator of V_cv=V_th. The charging procedure is completed when the maximum current flowing through the different modules achieves the threshold I_th (which in this case is assumed to be equal to 0.1MI_1C). In order to implement such a procedure, the battery pack needs to be equipped with a system of switches which allows to commute from the first phase to the second phase of the CC-CV. 


For the CC-CV charging protocols, the temporal evolution of the state of charge, voltage, temperature, and currents applied to the different cells are shown in Figures <ref>, <ref>, <ref>, and  <ref>. Two scenarios were considered which differ in the value of the current applied during the constant current phase. The constant current of I_cc=MI_1C (dotted line) results in high temperature constraint violation in most of the cells, with a charging time of 3960s. On the other hand, for the lower constant current (dashed line) of I_cc=0.85MI_1C, the temperature constraint is satisfied for each cell, but the charging time increases significantly (4360s). Moreover, the value of the constant current I_cc which guarantees the satisfaction of the temperature constraint must be found experimentally and can change according to the external environment conditions as well as with the battery ageing and degradation. 

Table <ref> compares the charging time, computational time, and maximum temperature and voltage for the sMPC, nMPC, and two CC-CV charging protocols. The MPC algorithms have the same charging times and maximum temperature and voltage, while sMPC required only about 6% of the on-line computational time.


















































 §.§ Scaling of the Computational Time for Increasing Number of Cells

The above simulations showed that the sMPC has much lower on-line computational cost than nMPC for a battery pack of 4 cells. Figure 7 displays the mean computational time for sMPC and nMPC for an increasing the number of series and parallel connections. sMPC is about an order of magnitude faster than nMPC, which is a significant savings in on-line computational cost when dealing with large battery packs. 






















  §.§.§ Optimal Control of an Electric Motorbike Battery Pack

This section demonstrates the applicability of the proposed sensitivity-based MPC to the control of large battery packs, for which nMPC has high on-line computational cost. In this case study, the task is to control the battery of a fully electric motorbike, namely, the electric Vespa Piaggio with a stored energy of 86Ah and a nominal voltage of 48V (i.e., a configuration with 156 cells, arranged in 13 series-connected modules, each with 12 parallel-connected cells, for the Kokam SLPB 75106100). The mean computational time required for each iteration of nMPC was 250s,  which is incompatible with the desired sampling time (T_s=40s). On the other hand, the mean computational time of sMPC was 30s, which is less than the sampling time.

 


§ CONCLUSIONS

This work addresses the optimal charging of a battery pack composed of several cells arranged in series and parallel connections. Each cell is described through an electrochemical model that includes kinetics, mass transport, and thermal effects in order to capture the internal physicochemical phenomena. A nonlinear model predictive control (MPC) formulation is formulated that achieves high performance while ensuring constraint satisfaction. An alternative sensitivity-based MPC formulation is proposed that has very similar closed-loop performance but greatly reduces the online computational cost, which makes optimal model-based control suitable for a real-time implementation on a battery pack composed of dozens of cells. The effectiveness of the strategy is demonstrated. The sensitivity-based MPC is successful in providing real-time optimal charging for a fully electric motorbike composed by 156 cells.






