Real-time Nonlinear Model Predictive Control (NMPC) Strategies using Physics-Based Models for Advanced Lithium-ion Battery Management System (BMS)

Optimal operation of lithium-ion batteries requires robust battery models for advanced battery management systems (ABMS). A nonlinear model predictive control strategy is proposed that directly employs the pseudo-two-dimensional (P2D) model for making predictions. Using robust and efficient model simulation algorithms developed previously, the computational time of the nonlinear model predictive control algorithm is quantified, and the ability to use such models for nonlinear model predictive control for ABMS is established. © 2020 The Author(s). Published on behalf of The Electrochemical Society by IOP Publishing Limited. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/ by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/ 1945-7111/ab7bd7]


List of variables for P2D model
Related to the solid-phase potential 2 Related to the liquid-phase potential Lithium-ion batteries are now ubiquitous in applications ranging from cellphones, laptops, electric vehicles, and even electric flights. Safety and long recharging times along with capacity and power fade remain some of the major concerns for lithium-ion batteries. Advanced battery management systems (ABMS) that can counter these issues and implement optimal usage patterns are critical for efficient use of batteries. Various optimal charging strategies have been proposed by researchers in recent times that minimize battery degradation or charge the batteries faster. [1][2][3][4] However, most of these strategies have been derived either using reduced-order physicsbased models, or implemented as open-loop control profiles based on offline calculations. While model order-reduction simplifies the governing model and decreases the numerical stiffness of the underlying full model, it often comes at the cost of simplification of actual physics of the system. Additionally, lithium-ion battery models have uncertainties due to low confidence in estimated system parameters, or parameters that can change with time, which makes open-loop control strategy less effective and necessitates a closedloop (feedback) control for optimal system performance.
Model predictive control (MPC) is an advanced closed-loop control strategy, which due to its characteristics, can be incorporated into ABMS to derive optimal charging protocols. This framework, while satisfying physical and operational constraints, evaluates the control objective based on the future predictions of the plant. Various MPC techniques deriving optimal charging profiles using approximated porous electrode pseudo-2-dimensional (P2D) models have been published in the literature. Xavier et al. proposed MPC strategies for controlling lithium-ion batteries using equivalent circuit models. 5 Torchio et al. proposed a linear MPC strategy based on the input-output approximation of the P2D model. 6 Torchio et al. also proposed health-aware charging protocols for lithium ion batteries using a linear MPC algorithm along with piecewise linear approximation and linear time-varying MPC strategies for lithiumion batteries. 7,8 Klein et al. proposed a nonlinear MPC framework based on a reduced-order P2D model. 9 Lee et al. proposed an MPC algorithm for optimal operation of an energy management system containing a solar photovoltaic panel and batteries connected to a local load in a microgrid. 10 Liu et al. derived nonlinear MPC profiles for optimal health of lithium ion batteries using a full single-particle model. 11 Traditionally, the high computational cost of online calculations has been often cited as one of the main reasons for not using detailed P2D models in MPC formulations.
While nonlinear MPC formulations based on battery models have been developed before, we propose implementation of such strategies using more robust and efficient numerical solvers along with reformulated models, allowing us to significantly reduce the computational time of this technique and enabling their use in real-time ABMS platforms. In this work, we design a nonlinear MPC controller capable of deriving optimal charging profiles using the detailed isothermal P2D model in real-time. The nonlinear model predictive control scheme is summarized in section 2, followed by a discussion on the numerical optimization approach used for solving the optimal control problem within the MPC framework. We then implement the nonlinear model predictive control technique to derive optimal charging protocols for the thin film nickel hydroxide electrode, discussed in section 3, for setpoint tracking objectives. Section 4 demonstrates the nonlinear model predictive controller designed by using the detailed reformulated P2D model. Section 5 analyses the effect of tuning parameters on the performance of the designed controller followed by a description of the computational efficiency achieved by the controller while using the detailed P2D model. Section 6 summarizes and outlines the future directions of the work.

Nonlinear Model Predictive Control
Model Predictive Control (MPC) is a multivariable control strategy with an explicit constraint-handling mechanism. This strategy involves generating a sequence of manipulated inputs over a control horizon, which optimizes a defined control objective over a prediction horizon, using an explicit process model. 12,13 If a nonlinear process model is used within the framework, then this strategy is termed as Nonlinear Model Predictive Control (NMPC). 12,14 A nonlinear optimal control formulation 15 related to the NMPC strategy given in literature is Formulation-I: • Equation 2 defines the equality constraints that describe the dynamics of the nonlinear plant denoted by a set of differential algebraic equations (DAEs), where functions f and g describe the differential and algebraic relations, respectively, ( ) x t represents the states of the plant, and ( ) u t represents the input signals to the plant.  The optimal control problem in Formulation I, defined by Eqs. 1-4 is a constrained dynamic optimization which can be solved using direct or indirect methods. 15,16 This work implements a direct method referred to as sequential dynamic optimization. The resulting nonlinear program (NLP) used in this method is discussed in the next subsection.
NMPC optimal control problem using sequential dynamic optimization.-In Formulation-I the decision variable ( ) u t of the optimal control problem is a continuous variable as shown in Fig. 1a.
In sequential dynamic optimization, the infinite-dimensional optimal control problem is reduced to a finite-dimensional NLP through discretization of the input signal ( ) u t to N discrete node points, where N is defined as the total time T f over the sampling time Dt 16 In this method, the input signal ( ) u t is assumed to be a piecewise constant at each sampling time instant Dt as shown in Fig. 1b. To formulate the finite-dimensional NLP, the input signal is discretized as Î  U , p a p-dimensional real-valued vector, where p is the prediction horizon. The reformulated finite-dimensional NLP is Formulation-II: is the objective function, minimizing the cost function j, which is solved for a finite number of optimal input signals, at time instants t k for = ¼ k p 1, , . The cost function j in Eq. 5 for the setpoint tracking objective is written in the discrete-time formulation as where v k denotes the controlled variable at the time instant t , k v set denotes its desired setpoint, U k denotes the predicted optimal manipulated variable at the time instant t , k and Q and R denote weighting parameters for setpoint tracking and input variations, respectively.  Formulation II (Eqs. 5-9) can now be numerically solved using an optimizer along with a robust numerical integrator (DAE solver). In any optimal control problem within the NMPC framework, the optimizer is treated as an "outer-loop" and the DAE solver is treated as an "inner-loop." At each iteration in optimization, the vector of the decision variables U provided by the optimizer is fed to the DAE solver to simulate the model for a finite number of time instants. The state variable trajectories from the DAE solver are then used to evaluate the objective and constraint functions. These functional values are sent to the optimizer, which provides an updated vector of the decision variables for the next optimization calculation. The resulting sequence of simulation and optimization iterations is also referred to as sequential simulation-optimization. 16 Receding horizon approach.-In the MPC framework, after obtaining the "p" optimal inputs, the first optimal input is sent to the plant. The resulting feedback from the plant is incorporated by estimating the states to minimize the plant-model mismatch, upon which the resultant NLP is solved recursively at each sampling Given: Mathematical model f, initial condition ( ) x 0 , prediction horizon p, control horizon m, sampling time Dt, and weighting matrices Q and R Step 1: At the current sampling time t , Step 2: Solve Formulation II for a sequence of m optimal input variables and inject the input to the plant Step 4: At the sampling time instant + t , k 1 obtain the plant measurement y m Step 5: Corresponding to y , m estimate the states ( ) + * x t k 1 (this work assumes full state feedback, for which all the states are measurable) Step 7: Shift the prediction horizon p forward and repeat Step 1 Figure 2. Schematic representation of a model predictive controller.
instant. This recursive method is also termed as 'receding horizon control' 13 which is described by the algorithm in Table I. A pictorial illustration the NMPC algorithm is shown in Fig. 2. The design parameters for the NMPC formulation are, (i) prediction horizon p, (ii) control horizon m (m is specified so that m ⩽ p), (iii) the sampling period Dt, and (iv) weighting parameters [ ] Q R , (in the objective function of Formulation II, Eq. 10) for setpoint tracking and input variations. The weighting parameter R makes the response of NMPC sluggish. In this work, it is taken as zero to enable fast charging strategy.
In real systems, it might not be possible to measure all the states of the system. In that case, the states corresponding to the new plant measurement at sampling instant + t k 1 need to be estimated (Step 5 in Table I). In practice, nonlinear state estimators such as Extended Kalman Filter (EKF) or Moving Horizon Estimator (MHE) are used to estimate the states for the control algorithm. The use of these estimators is under investigation by the authors and will be reported in the future work. Here, the model is differentiated from the plant by introducing model uncertainty by perturbing certain parameters of the system, as described in the Appendix.

Thin Film Nickel Hydroxide Electrode Model
To illustrate the implementation of the control scheme, a twoequation model representing the galvanostatic charge process of a thin film nickel hydroxide electrode 17 is described by the DAE model: where the dependent variable y represents the mole fraction of nickel hydroxide and z represents the potential difference at the solid-liquid interface. The parameters used in the model Eqs. 11-14 are in listed in Table II.
Control objective.-The control objective is defined as a setpoint tracking problem. According to the control objective, an optimal current density profile is computed that drives the mole fraction (the controlled variable) from its initial state to the desired setpoint. While fulfilling the objective, the bounds are simultaneously imposed on the current density.
The defined control objective can be formulated as the NLP (for scalar y): Formulation-III: subject to the constraints: model differential and algebraic Eqs. 11-14 is the setpoint tracking control objective where ( ) y k denotes the nickel hydroxide mole fraction for all k sampling instants over the prediction horizon p, with each sampling instant of time Dt, and y set denotes the desired set point for the nickel hydroxide mole fraction.
• Equation 16 defines the bounds on applied current density ( ) I app over the prediction horizon p, and I app max denotes the upper bound on the applied current density.
Simulation results.-The NLP Formulation III is solved using NMPC algorithm discussed in Table I. The closed-loop trajectories of nickel hydroxide mole fraction, potential difference at the solidliquid interface, and applied current density are shown in Fig. 3. The controller tracks the nickel hydroxide mole fraction (controlled variable) to a set point at 0.9. This case study used = Q 1 and = R 0, and A cm 2 was considered to account physical dissimilarities between different charging units. For satisfying this control objective, the controller is designed with a prediction horizon p of 3 sampling periods, control horizon m of 3 sampling periods, and sampling period Dt of 100 s. To study the robustness of the controller, model-plant mismatch is introduced by increasing the mass of the active material W by 10% in the plant simulation.
The controller validates the observation that a higher maximum input current density results in the mole fraction of the nickel hydroxide electrode reaching its reference value more quickly than for a lower maximum current density.
Bounds on additional state variables (such as voltage (z) in this example) can also be introduced in the NMPC framework. Such bounds will be illustrated in detail in the next section, where the implementation of the NMPC strategy using the pseudo-2-dimensional (P2D) model of a lithium-ion battery is discussed.

Pseudo 2-Dimensional (P2D) Model of a Lithium-Ion Battery
The Pseudo-Two-Dimensional (P2D) model is one of the most widely used physics-based electrochemical models for lithium-ion batteries. 18 The complete set of partial differential algebraic equations (PDAEs) describing the governing equations of the P2D model are given in Table AI  Assuming the battery to be limited by the anode capacity, the bulk SOC is calculated as the average of the volume-averaged solidphase lithium concentration across the negative electrode: where c s n max , denotes the maximum solid-phase concentration of lithium in the negative electrode, c s avg denotes the volume-averaged solid-phase concentration in each solid particle in the negative electrode, and L n denotes the length of the negative electrode of the battery. q min and q max are states of charge at fully discharged and charged states, that depend on the stoichiometric limits of the negative electrode. This choice of controlled variable illustrates the ability and speed of the NMPC algorithm. In general, the batteries are often limited by the lithium concentration in the positive electrode (cathode). Additionally, state variables such as cell voltage or temperature can also be used as controlled variables, as they can be measured directly.
Apart from the main reaction of lithium-ion intercalation, various side reactions occur during charging which may potentially damage the battery. 9,19,20 For example, anodic side reactions may deposit lithium on the surface of the negative electrode (lithium plating) thereby resulting in the subsequent loss of the battery's capacity. 20,21 The lithium plating occurs when the over-potential at the anode becomes negative. 21 As the open-circuit potential of the lithium plating side reaction is taken as 0 V (vs Li/Li+), the over-potential of the lithium plating side-reaction is defined as plating 1 2 It has been previously shown that lithium plating is more likely to occur at the anode-separator interface at high charging rates 19 ; hence we apply constraints only at the anode-separator interface throughout our analysis. As F 1 and F 2 are obtained as internal states of the P2D model, the anode over-potential can be tracked at any time during charging. By constraining the anode overpotential to be non-negative during charging, it is possible to restrict lithium plating side reaction, thereby mitigating battery degradation. The accuracy of the underlying model plays a vital role in predicting and thereby restricting the anode over-potential, as it cannot be directly measured. 20 Therefore, using a detailed physics-based model (P2D model) for BMS helps in minimizing battery degradation, thereby enabling the utilization of the battery to its full potential. where v k denotes the controlled variable at the time instant t , k in which the controlled variable is either SOC or voltage for the system considered; v set denotes the desired set point for SOC or voltage; and I app k , denotes the predicted optimal applied current density (input variable) at the time instant t . k The first term in Eq. 19 describes the setpoint tracking objective and the second term represents the changes in the applied current density. The weighting factor (Q) Figure 3. NMPC time profiles from Formulation III for (a) current density, (b) mole fraction, and (c) potential. The simulations are performed using "ode15s" from the MATLAB solver suite as the DAE solver and "fmincon sqp" as the NLP solver. for setpoint tracking is described by a scalar, due to the presence of a single controlled variable in the electrochemical system under study but can be a vector if there are multiple controlled variables.
For Li-ion batteries, the defined objective can be interpreted as deriving a charge current profile that drives and maintains the controlled variable at a desired operating condition. In doing so, it is desired to simultaneously enforce physical and operational constraints for the safe and optimal charging of a battery. With SOC as the desired controlled variable, the control objective (Eq. 19) is reformulated as the NLP with specific constraints, to obtain the optimal control problem ( ) * I app with Q and R are set as 1 and 0, respectively. The original governing PDAEs are spatially discretized using the strategy described in Northrop et al. 22 and the resulting DAEs of the reformulated P2D model are used as constraints. The convergence analysis on the spatial discretization strategy is discussed in Appendix.
The objective defined in Eq. 19 can also be viewed as a pseudo minimum charging time problem as it brings similar results compared to a battery fast-charge problem (a battery fast-charge problem is defined as finding the optimal charging strategy to charge a battery from an initial SOC to the desired SOC, in the shortest possible time, with given constraints on the voltage, current, temperature, overpotential, or other variables, for the same sample time).
Below is a discussion of the derivation of control profiles for various constraints employed on cell voltage and overpotential at the anode-separator interface. Model-plant mismatch and corresponding uncertainty in the model are introduced by changing the parameter values as shown in the Appendix. The tuning parameters and the bounds used are   • Equation 22 represents bounds on the overall cell voltage.
Imposing bounds on the overall voltage of the battery is essential for its safe operation. Every battery is rated by the battery manufacturer to be operated within a specified voltage window. Hence, for safety (and legal warranty issues imposed by the battery manufacturer in most cases), it is recommended to restrict the battery voltage within a finite window described by (22).
• Equation 23 describes the bounds on the applied current density over the prediction horizon p. This study considers an isothermal model to demonstrate the methodology and the computational time of the algorithm. However, additional constraints on other state variables such as temperature can also be included in the algorithm using a thermal model. Figure 4 shows the comparison between traditional CC-CV profiles and optimal control profiles obtained using Formulation IV. NMPC strategy drives and maintains the SOC at its setpoint of 100% while enforcing the bounds on the applied current density and cell voltage. Once the SOC reaches its setpoint, the controller progressively drops current density to zero as expected, thereby maintaining the desired setpoint conditions. This results in a control profile that qualitatively follows the traditional CC-CV profile until the desired SOC is reached, while essentially charging a battery to the final SOC in the shortest possible time. However, it should be noted that the overpotential for the lithium plating reaction at the anode-separator interface becomes negative (h < 0 plating ) at a certain time while charging. As discussed before, this behavior while charging might lead to the deposition of lithium on the surface of the negative electrode, leading to capacity fade and dendrite formation. Therefore, for ensuring safe operating conditions, constraints are imposed on the plating overpotential to avoid the regimes where h < 0, plating as in Formulation V.
1, , 28 In addition to the constraints described in Formulation IV, Eq. 28 describes the constraints on the lithium plating overpotential over the predictive horizon p. As previously discussed, this constraint mitigates battery degradation due to lithium plating. Figure 5 shows the comparison of the traditional CC-CV profiles and optimal control profiles obtained after adding the constraints on overpotential. The results in this case study show that the proposed manipulated variable profiles drive the controlled variable to a desired set point, in the least time possible, while enforcing constraints on the mechanisms which degrade the battery life. Achieving the same SOC levels using a conventional CC-CV charging profile will lead to negative side overpotential which might potentially degrade the battery performance. In other words, though conventional CC-CV protocols are time tested, the significance of optimal control profiles can be gauged when NMPC strategies are implemented while experimentally cycling the cells.
Servo problem.-The explicit time dependence of the stage cost/ control objective and equality and inequality constraints (comprising model equation constraints, input, and state variable constraints) allow for the incorporation of dynamic setpoint trajectories in NLP defined by (5)- (9). 15 In certain applications, it may be desirable for the batteries to experience specific dynamic voltage profiles. Here, NMPC results are presented for a time-varying setpoint on the voltage. Formulation-VI: where V set is given by the "red" dashed line in Fig. 6d. The controller, in this case, was designed with = = p m 3, 2, and D = t 30 s. Figure 6 shows the control profiles obtained for a dynamic setpoint trajectory.

Computational Details
Traditionally, the incorporation of a detailed physics-based model (P2D model) in BMS applications has been said to be computationally expensive due to their large simulation times. 7 Therefore, incorporation of such models for real-time simulation and control applications necessitates efficient, fail-proof and fast solvers. In our previous work, we demonstrated the simulation of the reformulated P2D model with computation time of 15 to 100 ms. 20,[22][23][24] This reduction in the simulation time facilitates the use of P2D model for real-time control applications using NMPC, as demonstrated by the results obtained from this work. All the results reported in this work are obtained using MATLAB. In this environment, the single optimization call to identify optimal current density for single prediction horizon using detailed P2D model was approximately 60 s. The detailed summary of the computation time (using MATLAB) for all cases is given in Table III. The computational time (including single optimization call and single model simulation call) for the NMPC strategy with P2D model will be Figure 6. NMPC time profiles for Formulation VI to identify optimal current density required to match dynamically varying set-points on cell potential. Included constraints on over-potential in Formulation-IV. Compared to Formulation-IV, there is change in current density profile to avoid lithium plating. The optimal control profiles is close to conventional CC-CV charging protocol. But optimal charging is always better strategy as it has the ability to avoid over charging as well as plating compared to conventional charging. VI Isothermal P2D

≈65 ≈0.66
In control theory, it is a servo problem. In practice, some of the applications can have dynamic cell potential profiles. This case studies the ability to implement NMPC strategies that manipulate current density to match dynamic set-points.
lower (≈2 s) when deployed in the C environment. The obtained computational efficiency demonstrates that a detailed P2D model can be used for real-time control applications of BMS. Such detailed models facilitate aggressive and optimal charging protocols, thereby extracting maximum performance from the cell.
Note.-The robustness of this sequential approach relies on the integration solver (odes15s in MATLAB or IDA in C) used in the nonlinear programming problems. In general, the isothermal and thermal battery models pioneered by John Newman are index-1 DAE's. ODE15s is numerical integrator in MATLAB that can handle only index-1 DAEs. There are more robust solvers for index-1 DAE's such as IDA in C developed by SUNDIALs or DASSL/DASSPL. 25 If pressure models are considered in addition to electrochemical models, the resulting DAE's are index-2 DAEs. 26 The best solvers for index-2 DAEs are RADAU. 27 The use of these solver requires the specification of exact initial conditions for the algebraic variables and also requires the identification of index 2 variables.
As of today, for higher index DAEs, the best option is to reformulate and reduce these DAEs to index-1 DAEs and then solve them using Pantelides Algorithm. The difficulty for higher index DAEs are limited to sequential approach. Even for simultaneous approach, there will be reduction in accuracy for higher index DAEs.

Summary
This article presents implementation of nonlinear model predictive control on physics-based battery models for deriving optimal charging protocols. We have shown that the designed NMPC controller is efficient in satisfying the given control objectives in the presence of different constraints on the internal state variable and applied current density. It is shown that the proposed controller, through constraining the plating overpotential, can efficiently derive health-conscious charging profiles while still charging the battery to the desired setpoint on SOC. Further, the effectiveness of the controller in tracking a dynamically varying setpoint is also demonstrated. This study demonstrates that a detailed P2D model can be incorporated in the design of ABMS for enabling real-time control of Li-ion batteries. While the objective has been formulated as set-point based on SOC or cell voltage, it can easily be modified to minimize the capacity fade over a charging period (provided that the capacity fade model is incorporated) or minimize the total charging time with constraints on the total charge stored, among others. The formulations discussed in this work are summarized in Table III. For future investigations, we plan to explore implementation of simultaneous numerical optimization strategies instead of sequential strategies for solving the NMPC optimal control problems. Simultaneous strategies, apart from being computationally less expensive, do not depend on a robust DAE solver for evaluating the objective and constraint functions. Further, path constraints through simultaneous strategy can be handled in a more efficient way and need not be approximated as with sequential approach. However, this requires careful and sufficient discretization strategies in time (number of elements, method of discretization, etc.) which will be reported in the future. Future publications will also report on the implementation of an output feedback NMPC, where a nonlinear state estimator is incorporated in the existing framework, for providing the full state information at each sampling instant.   Model Uncertainty (Model-Plant mismatch).-For this work, the plant is simulated by the same model equations. Uncertainty in the model (signifying error in the model), and a corresponding mismatch with the plant, is introduced by perturbing the model parameters compared to the plant parameters. Figure A2 shows the comparison of model vs plant dynamics for a simulation performed at 3C charge rate, using parameters listed in Table AIV. ORCID Suryanarayana Kolluri https://orcid.org/0000-0003-2731-7107 Richard D. Braatz https://orcid.org/0000-0003-4304-3484 Venkat R. Subramanian https://orcid.org/0000-0002-2092-9744 Figure A2. Model uncertainty introduced by changing the solid-phase diffusivity and conductivity kinetic rate constant of the positive electrode as mentioned in Table A