Adaptive Model Predictive Control-Based Energy Management for Semi-Active Hybrid Energy Storage Systems on Electric Vehicles

This paper deals with the energy management strategy (EMS) for an on-board semi-active hybrid energy storage system (HESS) composed of a Li-ion battery (LiB) and ultracapacitor (UC). Considering both the nonlinearity of the semi-active structure and driving condition uncertainty, while ensuring HESS operation within constraints, an adaptive model predictive control (AMPC) method is adopted to design the EMS. Within AMPC, LiB Ah-throughput is minimized online to extend its life. The proposed AMPC determines the optimal control action by solving a quadratic programming (QP) problem at each control interval, in which the QP solver receives control-oriented model matrices and current states for calculation. The control-oriented model is constructed by linearizing HESS online to approximate the original nonlinear model. Besides, a time-varying Kalman filter (TVKF) is introduced as the estimator to improve the state estimation accuracy. At the same time, sampling time, prediction horizon and scaling factors of AMPC are determined through simulation. Compared with standard MPC, TVKF reduces the estimation error by 1~3 orders of magnitude, and AMPC reduces LiB Ah-throughput by 4.3% under Urban Dynamometer Driving Schedule (UDDS) driving cycle condition, indicating superior model adaptivity. Furthermore, LiB Ah-throughput of AMPC under various classical driving cycles differs from that of dynamic programming by an average of 6.5% and reduces by an average of 10.6% compared to rule-based strategy of LiB Ah-throughput, showing excellent adaptation to driving condition uncertainty.


Introduction
As concerns about energy crisis and environmental issues mount, electric vehicles (EVs) have been considered as the most promising substitutes for internal combustion engine vehicles.EVs' performance is strongly effected by their energy storage sources (ESS) [1,2].ESSs for EVs must simultaneously possess high specific energy, high specific power and long cycle life, which is impossible for a single ESS.Consequently, hybrid energy storage system (HESS) have becomes favorable options.Among all ESSs, Li-ion batteries (LiB) are the most extensively used for their high specific energy, however, low specific power and short cycle lives limit their application [3].Compared with LiB, ultracapacitors (UCs) have longer cycle life and higher specific power [4].Therefore utilizing HESSs that include LiB and UC and possess both the high specific energy of LiB and high specific power of UC is an advisable choice.In this way, the LiB burden would be reduced remarkably by the assistance of UC, and LiB lifetime would be effectively extended.
fluctuation.CVX is an effective tool in dealing with optimization problems with multi-states, whereas, it is only applicable to convex cost functions with convex constraints, and its application in nonconvex models is limited.
From the literature described above, optimization-based EMSs are preferred over RBC-based EMSs for their optimality and adaptability to driving condition uncertainty, but due to the strong nonlinearity of semi-active HESS, the mentioned optimization methods are not applicable.In this paper, an AMPC-based EMS for semi-active HESS is proposed to handle the model nonlinearity.The proposed AMPC minimizes LiB Ah-throughput in the prediction horizon and ensures HESS operation within constraints.Three aspects of research are conducted on the proposed AMPC.Firstly, a module for generating control-oriented model is designed, including linearization, direct feedthrough elimination, discretization and state augmentation.Secondly, a time-varying Kalman filter (TVKF) is introduced as state estimator to improve the state estimation accuracy.In this way, both controloriented model and estimated states improve the control action accuracy calculated by the QP solver.Finally, sampling time, prediction horizon and scaling factors of AMPC are determined through simulation.The proposed EMS is created and verified on the Matlab/Simulink (The MathWorks, Natick, MA, US) platform.Compared with standard MPC (SMPC), TVKF reduces estimation errors by 1~3 order of magnitude, AMPC reduces LiB Ah-throughput by 4.3% under UDDS conditions, indicating superior model adaptivity.Furthermore, simulation results under various driving cycles differ from those of DP by an average of 6.5% and reduce an average of 10.6% compared to RBC, indicating that AMPC is able to minimize the LiB Ah-throughput online effectively when the HESS is nonlinear and the driving conditions are uncertain.
The rest of the paper is organized as follows: in Section 2, the detailed HESS model is introduced.In Section 3, the control-oriented model updating process is performed, and TVKF is specified.In Section 4, AMPC parameters are discussed, and AMPC is compared with SMPC, RBC and DP to verify the model adaptivity and adaptability to driving condition uncertainty.Moreover, the influence of the weights is researched and a suggested range is given.Section 5 presents the conclusions.

The Hybrid Energy Storage System Modeling
The structure of the semi-active HESS is shown in Figure 1, in which UC is connected to the DC bus directly.LiB is connected to the DC bus through a full-size bi-directional DC/DC converter.The LiB power output is determined by EMS, and the UC delivers rest power demand passively.

LiB Model
The LiB cell used in this paper is a Panasonic NCR 18650B (Panasonic Co., Ltd, Osaka, Japan), and the main parameters are listed in Table 1.

LiB Model
The LiB cell used in this paper is a Panasonic NCR 18650B (Panasonic Co., Ltd, Osaka, Japan), and the main parameters are listed in Table 1.The Partnership for a New Generation of Vehicles (PNGV) battery model is selected as the LiB equivalent circuit model [31], as shown in Figure 2. The Partnership for a New Generation of Vehicles (PNGV) battery model is selected as the LiB equivalent circuit model [31], as shown in Figure 2. LiB state of charge (BSOC) is defined as: where, BSOCinit is the value of BSOC when the simulation starts, which is set as 0.7.Qb is the LiB cell nominal capacity.
The OCV under different BSOC is obtained by experiment.When BSOC is between 0.2 and 0.9, the OCV increases nearly linearly with BSOC.Limiting BSOC within the range of [0.2, 0.9], and we linearized the relationship of BSOC and OCV to obtain E and Cbb: Thereafter, BSOC = Vbb/b, E = a [32].Vbb and Vbp are two states expressed by: 1 , Terminal voltage Vbl is expressed by: Based on Kirchhoff Voltage Laws, ib is obtained: where, Pb is the only manipulated variable of the HESS.
A hybrid pulse power characterization test is conducted, and the lumped parameters Cbp, Rbb and Rbp are fitted under different BSOC (0.1 as interval).As BSOC is limited to [0.2, 0.9], Rbb, Rbp and Cbp are evaluated by taking an average within the range.The bi-directional DC/DC converter is simplified In the PNGV model, E and V bb approximate the voltage source of the LiB, in which E is a constant evaluated as the open circuit voltage (OCV) value when the LiB is fully discharged.V bb represents a linear relationship between the variation of source voltage and variation of remaining source energy, just like a capacitor approximately, and C bb is corresponding capacity.C bp and R bp reflect the electrode polarization phenomenon; R bb is the equivalent dc resistance.V bl is the cell terminal voltage, P b is pack power, n b is cell number in pack.
LiB state of charge (BSOC) is defined as: where, BSOC init is the value of BSOC when the simulation starts, which is set as 0.7.Q b is the LiB cell nominal capacity.The OCV under different BSOC is obtained by experiment.When BSOC is between 0.2 and 0.9, the OCV increases nearly linearly with BSOC.Limiting BSOC within the range of [0.2, 0.9], and we linearized the relationship of BSOC and OCV to obtain E and C bb : Thereafter, BSOC = V bb /b, E = a [32].V bb and V bp are two states expressed by: Terminal voltage V bl is expressed by: Based on Kirchhoff Voltage Laws, i b is obtained: Energies 2017, 10, 1063 where, P b is the only manipulated variable of the HESS.A hybrid pulse power characterization test is conducted, and the lumped parameters C bp , R bb and R bp are fitted under different BSOC (0.1 as interval).As BSOC is limited to [0.2, 0.9], R bb , R bp and C bp are evaluated by taking an average within the range.The bi-directional DC/DC converter is simplified as an efficiency module with constant efficiency of 95%.

UC Model
As UC is directly connected to the electric machine (EM) controller, the working voltage must be limited to the range of the controller protection voltage which is 220 V~400 V in this paper.The cell voltage upper bound is set as 2.7 V, and based on the bound, the cell number is calculated as 150.According to the protection voltage lower bound as well as calculated cell number, the cell voltage lower bound is obtained as 1.5 V.The cell capacity selection is a sizing problem coupled with EMS optimization.Based on various literatures [33][34][35], pack capacity is chosen as 100 Wh.Based on the pack capacity, the cell capacity is calculated as 0.67 Wh.A Maxwell BCAP0650 is selected whose cell capacity is 0.66 Wh, the main parameters of which are listed in Table 2.The equivalent circuit of the UC is shown in Figure 3 [4].This model consists of two parts: the part for main energy storage (C 0 and C 1 ) and the part reflecting electrode kinetics (R ub , R up and C up ).In the main energy storage part, C 0 and C 1 are connected in parallel.C 0 is characterized by constant capacitance, and C 1 is characterized by capacitance varying linearly with the voltage V ub .They characterize the nonlinearity of the source jointly.In the electrode kinetics part, R ub represents the Ohmic resistance, while R up and C up represent polarization resistance and polarization capacitance separately.

UC Model
As UC is directly connected to the electric machine (EM) controller, the working voltage must be limited to the range of the controller protection voltage which is 220 V~400 V in this paper.The cell voltage upper bound is set as 2.7 V, and based on the bound, the cell number is calculated as 150.According to the protection voltage lower bound as well as calculated cell number, the cell voltage lower bound is obtained as 1.5 V.The cell capacity selection is a sizing problem coupled with EMS optimization.Based on various literatures [33][34][35], pack capacity is chosen as 100 Wh.Based on the pack capacity, the cell capacity is calculated as 0.67 Wh.A Maxwell BCAP0650 is selected whose cell capacity is 0.66 Wh, the main parameters of which are listed in Table 2.The equivalent circuit of the UC is shown in Figure 3 [4].This model consists of two parts: the part for main energy storage (C0 and C1) and the part reflecting electrode kinetics (Rub, Rup and Cup).In the main energy storage part, C0 and C1 are connected in parallel.C0 is characterized by constant capacitance, and C1 is characterized by capacitance varying linearly with the voltage Vub.They characterize the nonlinearity of the source jointly.In the electrode kinetics part, Rub represents the Ohmic resistance, while Rup and Cup represent polarization resistance and polarization capacitance separately.
USOC is obtained by: where, Vub,max is 2.8 V according to the manual.
Terminal voltage Vul is expressed as: Based on Kirchhoff Voltage Laws, iu is represented as: V ub and V up are two states calculated by: USOC is obtained by: where, V ub,max is 2.8 V according to the manual.
Energies 2017, 10, 1063 Terminal voltage V ul is expressed as: Based on Kirchhoff Voltage Laws, i u is represented as: where, P u is replaced by P d − P b , indicating that P b and P d decide the working conditions of UC jointly.
V ub , i u and V ul reveal strong nonlinearity, and standard state space model is insufficient to characterize the model precisely.

AMPC-Based Energy Management Implementation
When the C-rate is within 2C, BSOC is within the range of [0.2, 0.9], Ah-throughput is capable of representing numerically the ageing condition of batteries [36].The purpose of the EMS is: (1) to minimize LiB current (or plus current variation), thereby reduce LiB Ah-throughput and extend battery life; (2) to ensure HESS is operating within the constraints.This is a finite horizon optimization problem with constraints, and MPC is excellent for managing optimization problems like this.In SMPC, the control-oriented model is designed off-line, and remains unchanged during simulation.If the plant is linear or weak nonlinear, the model prediction accuracy is acceptable.However, as the proposed semi-active HESS is highly nonlinear, the constant prediction model isn't accurate, AMPC is used to address the degradation by adapting the prediction model for changing the working conditions.Intrinsically, AMPC has the same structure as SMPC but allows the model parameters and optimization cost function to evolve with time.At each interval, a state space model is obtained by linearizing the HESS model around current working conditions.The original HESS model is: To ensure i b and i u have real solutions, the following two expressions are added to outputs: As AMPC is able to set limits for inputs and outputs only, states V bb and V ub who have engineering constraints are also added to outputs.Besides, as state estimator is required, the model observability must be guaranteed.To this end, V bl and V ul are added to outputs as measured variables.By this means, two current sensors and two voltage sensors are needed.The final model is with inputs Energies 2017, 10, 1063 7 of 21

Linearization
At each control interval, the plant is linearized around current operation point (X t , U t ) by first approximation of Taylor expansion.The linearized plant is described as: .
where, the constant matrices f (X 0 ,U 0 ) − AX 0 − BU 0 and g(X 0 ,U 0 ) − CX 0 − DU 0 are set as measured disturbances with constant inputs of 1, and the new inputs are augmented as The state space model is augmented as: Matrix A, B, C and D are calculated as:

Eliminating Direct Feedthrough
Non-empty matrix D which causes direct feedthrough from inputs to outputs is unacceptable by AMPC.To eliminate the feedthrough, the linearized plant in Equation ( 13) is augmented by adding first-order filters with time constant T a to inputs, as shown in Equation ( 15).T a is chosen as one-tenth or smaller of the AMPC sampling time [37,38].With this method, the original inputs are turned into new states and new inputs are identical to old ones except a slight delay: The linearized model is augmented by: Augmented states now are and new inputs are U a = [P ba , P da , M da ].The delay in the inputs intrinsically causes a slight model mismatch, which AMPC is good at dealing with.

Model Discretization
The Zero-Order Hold method is used to discretize the state space model with sample time Ts, the discrete state space is represented as: where: Energies 2017, 10, 1063 The matrix exponential of G (Ts) is calculated using Padé approximant and integral of H (Ts) is implemented using Simpson's Rule.

LiB Current Fluctuation Suppression
To avoid LiB damage caused by enormous current fluctuation, the current variation is also considered as a control target [27].Here we define a state: Based on the discrete state space, and output increments are represented as: ∆i b is expressed by the second row of ∆Y, in which coefficients vector of X a is expressed by G dib , coefficients vector of U a is expressed by H dib .∆i b is expressed as: The discrete state space model is augmented again: The final state space model for AMPC is obtained with internal states as well as control-oriented model matrices G A , H A and C A .

Optimization Problem
The AMPC is designed based on the model shown in Equation (22), in which P ba is the only manipulated variable, P da and M da are measured disturbances, i b , i u , V bl and V ul are measured outputs, V bb , V ub , y b , y u and ∆i b are unmeasured outputs.
The AMPC solves a QP problem at each control interval and determines the control action of next interval.The QP problem here includes outputs reference tracking and constraints violation, the cost function to be minimized is as follows [39,40]: where, n y is the number of outputs, p is prediction horizon, w j is the weight of jth output, ρ ε is the constraints violation penalty weight, and ε is the slack variable at interval k.J involves a trade-off between the output reference tracking and constraint violation by weighting, and corresponding QP decision is [U a (k|k) U a (k + 1|k to be activated, and the controller performance can be substandard.S y,j and s u are scaling factors (SFs) whose roles are to scale inputs and outputs to the same magnitude.ECRs are nonnegative parameters used to soften inputs and outputs constraints, the larger the ECRs are, the greater constraints violation are allowed to obtain optimal solution.Suppose H Au is coefficient matrix of P ba , and H Av is coefficient matrix of [P da , M da ], the model in Equation ( 22) is rewritten as: Assign U au = P ba and U av = [P da , M da ], based on Equation ( 23), the predicted output is: . . . Define Then substitute ∆U au into Equation ( 25), the predicted output is expressed as: where, H Au is coefficient matrix of P ba from H A , H Av is coefficient matrix of [P da , M da ] from H A .
As shown in Equation ( 26), G A , H A (H Au and H Av ) and C A are used by QP solver for calculation.ECRs for inputs and outputs are decided based on whether their upper or lower bounds are allowed to be violated or how great their bounds are allowed to be violated.Constraints for P b are not set for they are limited intrinsically by the limitation of V bb and i b , and their ECRs are set as 1 by default.V bb is limited to the range of 0.18 V~0.81 V based on BSOC constraints, both constraints leave some margin and can be violated, ECR 1,min and ECR 1,max are set as 1.Constraints for i b are −3.2A~6.4A according to the manual, the lower bound is strictly limited by the manufacturer, while the upper bound can be set as high as 10 A in our experiment, ECR 2,min is 0 and ECR 2,max is 1.Constraints for V bl and V ul are not set for they are intrinsically limited by other states, corresponding ECRs are set as 1.Outputs y b and y u must be non-negative, and the constraints are not allowed to be violated.Considering linearization error, the lower bounds for y b and y u are set as 0.5, ECR 3,min and ECR 6,min are 0. V ub is limited to meet the voltage range of DC bus, cell voltage range is thus between 1.5 V~2.7 V.As the energy stored in the UC increases exponentially with V ub , the upper bound would be violated slightly by absorbing much energy, which is acceptable, ECR 4,max is 1.While the lower bound would be violated severely by even delivering little energy, which is unacceptable, ECR 4,min is 0. Constraints for i u are −1200 A~1200 A and can be violated, for huge current do not harm the UC cell obviously, ECR 5,min and ECR 5,max are 10.Constraints for ∆i b are not set for the constraints of i b intrinsically avoid severe variation of ∆i b , meanwhile weights for i b and ∆i b prevent them from wide variation.

Controller State Estimation
QP solver receives states for calculation, as shown in Equation (26).As HESS states can't be measured directly, a state estimator is needed.Currents and terminal voltages of LiB and UC are related to all four states to be estimated, and they are all measured variables.By gathering voltage and current signals, the HESS observability is guaranteed.TVKF is applied to construct the estimator.Different from the standard Kalman filter, TVKF calculates the gain matrix with varying model matrices G A , H A and C A .The model in Equation ( 22) is augmented with noise vector w(k) and v(k) [26]: where, w(k) is process noise with covariance matrix Q, its gaining matrix is the same as inputs, and v(k) is measurement noise with covariance matrix R. Q is n-by-n diagonal matrix, n is equal to the number of states.R is m-by-m diagonal matrix, m is equal to the number of measured outputs.Suppose w(k) and v(k) are driven by white noise with unity gain, Q =H A H T A , and R is diagonal matrix with diagonal elements [SF 2 , SF 3 , SF 6 , SF 7 ] considering influence of scaling.
The estimation process at control interval k is as follows: (1) Prediction updating: (2) Measurement correction: where XA (k − 1|k − 1) is the corrected states prediction of previous interval, XA (k|k − 1) is the uncorrected states prediction of current interval, P(k − 1|k − 1) is the corrected states estimation error covariance matrix of previous interval, P(k|k − 1) is the uncorrected error covariance matrix of current interval, K g (k) is gain matrix for measurement correction, XA (k|k) is the corrected states prediction of current interval, P(k|k) is the corrected error covariance matrix of current interval.It's obvious that time-varying G A , H A and C A influence the value of K g significantly.

Measured Disturbance
The control-oriented model of AMPC has two measured disturbances: P da and M da .M da is a constant set as 1, while P da is HESS power demand determined by vehicle power balance equation shown below: where, v is vehicle velocity from driving cycle, m is vehicle mass, f is rolling resistance coefficient, C D is drag coefficient, A is windward area, δ is correction coefficient of rotating mass, a is vehicle acceleration.
As AMPC solves the QP problem with the assumption that measured disturbances in the prediction horizon are known, the velocity needs to be predicted.Yet, model prediction inevitably causes errors, and the predicted power demand wouldn't improve AMPC efficiency, and sometimes it may even deteriorate the efficiency if the velocity prediction is not that accurate.As a result, current power demand is applied to AMPC, and remains unchanged in prediction horizon.The whole process of the AMPC is shown in Figure 4, in the red box is the AMPC.

Sampling Time, Prediction Horizon and Scaling Factor
As sampling time (Ts) decreases, rejection of unmeasured disturbance as well as model error improves.Once prediction duration (pd) and Ts are determined, prediction horizon (ph) is calculated by ph = pd/Ts.ph decides how many steps the model prediction is conducted, as ph increases, the prediction accuracy gradually decreases.ph also decides the size of the matrices in the QP problem, as ph increases, the matrices become large and the computation burden is heavy.The AMPC with different Ts and ph is conducted under UDDS, and the results are shown in Table 3.Finally, Ts is set as 0.1 s, pd is chosen as 2 s, ph is 20.Scaling factors (SFs) are set to adjust inputs and outputs to the same magnitude, and improper SFs may cause inaccurate states estimation as well as poor QP solutions.As the magnitudes of HESS inputs and outputs vary a lot, it is important to choose proper SFs.SFs are usually determined empirically as their respective spans, whereas the simulation results may not be optimal.Based on spans of inputs and outputs, SFs of inputs [Pb, Pd, Md] are initially selected as [10,000, 10,000, 1], and SFs for outputs [Vbb, ib, Vbl, yb, Vub, iu, Vul, yu, ∆ib] are initially selected as [0.0005, 1, 1, 1, 0.01, 10, 1, 1, 1].Among them, SF of Pb and iu are found to be sensitive to the simulation results.Table 4 shows the results under different SFs:

Sampling Time, Prediction Horizon and Scaling Factor
As sampling time (Ts) decreases, rejection of unmeasured disturbance as well as model error improves.Once prediction duration (p d ) and Ts are determined, prediction horizon (p h ) is calculated by p h = p d /Ts.p h decides how many steps the model prediction is conducted, as p h increases, the prediction accuracy gradually decreases.p h also decides the size of the matrices in the QP problem, as p h increases, the matrices become large and the computation burden is heavy.The AMPC with different Ts and p h is conducted under UDDS, and the results are shown in Table 3.Finally, Ts is set as 0.1 s, p d is chosen as 2 s, p h is 20.Scaling factors (SFs) are set to adjust inputs and outputs to the same magnitude, and improper SFs may cause inaccurate states estimation as well as poor QP solutions.As the magnitudes of HESS inputs and outputs vary a lot, it is important to choose proper SFs.SFs are usually determined empirically as their respective spans, whereas the simulation results may not be optimal.Based on spans of inputs and outputs, SFs of inputs [P b , P d , Md] are initially selected as [10,000, 10,000, 1], and SFs for outputs [V bb , i b , V bl , y b , V ub , i u , V ul , y u , ∆i b ] are initially selected as [0.0005, 1, 1, 1, 0.01, 10, 1, 1, 1].Among them, SF of P b and i u are found to be sensitive to the simulation results.Table 4 shows the results under different SFs: Finally, SFs for P b and i u are reselected as 1000 and 30.

Model Adaptivity Verification
To validate the AMPC adaptivity to a nonlinear model, SMPC is introduced for comparison.SMPC is designed off-line with the same process of AMPC around HESS initial working points, and keeps constant through the simulation.A steady state Kalman filter (SSKF) is used as the estimator of SMPC.
SSKF estimation of SMPC is shown in Figure 5, TVKF estimation of AMPC is shown in Figure 6.We can see from Figure 5 that u ub estimation gradually deviates from the real value, and the error is accumulative.This is because the gain matrix solved by SSKF is constant through the simulation, an the state estimation error is not applied to update the gain matrix.For the same reason, excessive deviation happens in the remaining four estimations.In Figure 6, TVKF exhibits a pretty precise estimation.TVKF of AMPC solves gain matrix with updated G A , H A and C A at each control interval, it is obvious that varying G A , H A and C A could approximate the real HESS model better, TVKF thus significantly outperforms SSKF in estimation accuracy.
Average estimation errors of SSKF and TVKF are listed in Table 5. Estimation errors of TVKF are lower than those of SSKF in order of 1~3 magnitude, quantitatively proving TVKF is more accurate and model adaptive.Finally, SFs for Pb and iu are reselected as 1000 and 30.

Model Adaptivity Verification
To validate the AMPC adaptivity to a nonlinear model, SMPC is introduced for comparison.SMPC is designed off-line with the same process of AMPC around HESS initial working points, and keeps constant through the simulation.A steady state Kalman filter (SSKF) is used as the estimator of SMPC.
SSKF estimation of SMPC is shown in Figure 5, TVKF estimation of AMPC is shown in Figure 6.We can see from Figure 5 that uub estimation gradually deviates from the real value, and the error is accumulative.This is because the gain matrix solved by SSKF is constant through the simulation, an the state estimation error is not applied to update the gain matrix.For the same reason, excessive deviation happens in the remaining four estimations.In Figure 6, TVKF exhibits a pretty precise estimation.TVKF of AMPC solves gain matrix with updated GA, HA and CA at each control interval, it is obvious that varying GA, HA and CA could approximate the real HESS model better, TVKF thus significantly outperforms SSKF in estimation accuracy.Average estimation errors of SSKF and TVKF are listed in Table 5. Estimation errors of TVKF are lower than those of SSKF in order of 1~3 magnitude, quantitatively proving TVKF is more accurate and model adaptive.Figure 8 shows the SMPC and AMPC comparison for ib, Figure 8b,c is zooms of Figure 8a.In Figure 8b, the LiB current of SMPC is charging and discharging dramatically, which is not observed in AMPC.In Figure 8c, the LiB current of SMPC is fluctuating remarkably around 0 when HESS is    Average estimation errors of SSKF and TVKF are listed in Table 5. Estimation errors of TVKF are lower than those of SSKF in order of 1~3 magnitude, quantitatively proving TVKF is more accurate and model adaptive.Figure 8 shows the SMPC and AMPC comparison for ib, Figure 8b,c is zooms of Figure 8a.In Figure 8b, the LiB current of SMPC is charging and discharging dramatically, which is not observed in AMPC.In Figure 8c, the LiB current of SMPC is fluctuating remarkably around 0 when HESS is Figure 8 shows the SMPC and AMPC comparison for i b , Figure 8b,c is zooms of Figure 8a.In Figure 8b, the LiB current of SMPC is charging and discharging dramatically, which is not observed in AMPC.In Figure 8c, the LiB current of SMPC is fluctuating remarkably around 0 when HESS is charging.Furthermore, the LiB current of SMPC distinctly goes beyond that of AMPC.All the situations mentioned above lead to an increase of LiB Ah-throughput, as shown in Figure 8d.The phenomena indicate that time-varying prediction model matrices as well as accurate state estimation by TVKF significantly improve accuracy of control action calculated by QP solver of AMPC.Besides, as SMPC isn't able to estimate states precisely, the constraints on V ub and y u may be violated, resulting in SMPC failure, which will cause damage to HESS.
Energies 2017, 10, 1063 14 of 21 charging.Furthermore, the LiB current of SMPC distinctly goes beyond that of AMPC.All the situations mentioned above lead to an increase of LiB Ah-throughput, as shown in Figure 8d.The phenomena indicate that time-varying prediction model matrices as well as accurate state estimation by TVKF significantly improve accuracy of control action calculated by QP solver of AMPC.Besides, as SMPC isn't able to estimate states precisely, the constraints on Vub and yu may be violated, resulting in SMPC failure, which will cause damage to HESS.By comparison, the superior model adaptivity of AMPC is validated.

Driving Condition Uncertainty Adaptation Verification
To verify the adaptation to uncertain driving conditions, AMPC is simulated under various classic driving cycles, DP and RBC are constructed for comparison. .

DP and RBC Description
As DP solution time increases exponentially with the number of states, to reduce computation burden, USOC is selected as the only state.Within DP, LiB and UC are simplified as resistor capacitor (RC) circuits, USOC discretization step size is 0.001, and the cycle discretization step size is 1 s.To numerically understand the efficiency of different EMSs on extending LiB life, DP cost function includes LiB current absolute only.The cost function is: where, ib is calculated by: where, Pb is calculated by Pd − Pu, in which Pu is obtained from USOC variation.The DP process is shown in Figure 9. Forward iteration calculation is conducted.LiB Ahthroughput from the initial USOC point to different end USOC points are calculated, and the minimal solution as well as corresponding USOC path are recorded.Due to USOC discretization and model simplification, there is slight deviation between the DP result and real optimum value of the semi- By comparison, the superior model adaptivity of AMPC is validated.

Driving Condition Uncertainty Adaptation Verification
To verify the adaptation to uncertain driving conditions, AMPC is simulated under various classic driving cycles, DP and RBC are constructed for comparison.

DP and RBC Description
As DP solution time increases exponentially with the number of states, to reduce computation burden, USOC is selected as the only state.Within DP, LiB and UC are simplified as resistor capacitor (RC) circuits, USOC discretization step size is 0.001, and the cycle discretization step size is 1 s.To numerically understand the efficiency of different EMSs on extending LiB life, DP cost function includes LiB current absolute only.The cost function is: where, i b is calculated by: where, P b is calculated by P d − P u , in which P u is obtained from USOC variation.The DP process is shown in Figure 9. Forward iteration calculation is conducted.LiB Ah-throughput from the initial USOC point to different end USOC points are calculated, and the minimal solution as well as corresponding USOC path are recorded.Due to USOC discretization and model simplification, there is slight deviation between the DP result and real optimum value of the semi-active structure.Yet, considering that more accurate computation companies with much heavier computation burden and this result is accurate enough, it is adopted as comparison.
active structure.Yet, considering that more accurate computation companies with much heavier computation burden and this result is accurate enough, it is adopted as comparison.The RBC process is shown in Figure 10a, in which Pd and USOC are inputs, Pb is the output.Pb,avg represents the average power delivered by LiB, Pb,char represents the complementary power supplied from LiB to UC to maintain adequate energy for coming peak power demand.Both of them are tuned and keep constant for all driving cycles.Hysteresis control is shown in Figure 10b, which is designed for preventing USOC from frequent fluctuation round the bounds.The USOC range for hysteresis control is determined based on the UC capacity.As the capacity increases exponentially with USOC, the upper range could be set smaller compared with the lower range.Besides, considering the DC bus voltage constraint, the upper range is chosen as [0.9, 0.92], the lower range is chosen as [0.54, 0.59].

EMSs Results Comparison
Figure 11 shows the AMPC result under UDDS.From Figure 11a, when the HESS power demand is negative (charging), the UC absorbs almost all the braking energy, which is exactly what the AMPC expect.When the demand is positive (discharging), Vub frequently reaches the lower bound.This situation is dealt with by AMPC very softly, for Vub lower bound isn't allowed to be violated.In addition, the AMPC operational mechanism of dealing with hard constraints is shown in Figure 11b,c, which are a zoom from 450 s to 470 s of Figure 11a,d, respectively.We can see from Figure 11b that at about 453 s, when Vub is around 1.65 V and Pd is 20 kW, Pb starts to increase, while Pu gradually decreases to cope with the upcoming lower bound of Vub.After about 10 s, Vub reaches the lower bound softly without going beyond it.By this means, AMPC avoids any hard constraints violation.This mechanism is essential to Ah-throughput minimization, as each violation of Vub will cause an extra charging process from LiB to UC, which significantly increases LiB Ah-throughput.The RBC process is shown in Figure 10a, in which P d and USOC are inputs, P b is the output.P b,avg represents the average power delivered by LiB, P b,char represents the complementary power supplied from LiB to UC to maintain adequate energy for coming peak power demand.Both of them are tuned and keep constant for all driving cycles.Hysteresis control is shown in Figure 10b, which is designed for preventing USOC from frequent fluctuation round the bounds.The USOC range for hysteresis control is determined based on the UC capacity.As the capacity increases exponentially with USOC, the upper range could be set smaller compared with the lower range.Besides, considering the DC bus voltage constraint, the upper range is chosen as [0.9, 0.92], the lower range is chosen as [0.54, 0.59].
active structure.Yet, considering that more accurate computation companies with much heavier computation burden and this result is accurate enough, it is adopted as comparison.The RBC process is shown in Figure 10a, in which Pd and USOC are inputs, Pb is the output.Pb,avg represents the average power delivered by LiB, Pb,char represents the complementary power supplied from LiB to UC to maintain adequate energy for coming peak power demand.Both of them are tuned and keep constant for all driving cycles.Hysteresis control is shown in Figure 10b, which is designed for preventing USOC from frequent fluctuation round the bounds.The USOC range for hysteresis control is determined based on the UC capacity.As the capacity increases exponentially with USOC, the upper range could be set smaller compared with the lower range.Besides, considering the DC bus voltage constraint, the upper range is chosen as [0.9, 0.92], the lower range is chosen as [0.54, 0.59].

EMSs Results Comparison
Figure 11 shows the AMPC result under UDDS.From Figure 11a, when the HESS power demand is negative (charging), the UC absorbs almost all the braking energy, which is exactly what the AMPC expect.When the demand is positive (discharging), Vub frequently reaches the lower bound.This situation is dealt with by AMPC very softly, for Vub lower bound isn't allowed to be violated.In addition, the AMPC operational mechanism of dealing with hard constraints is shown in Figure 11b,c, which are a zoom from 450 s to 470 s of Figure 11a,d, respectively.We can see from Figure 11b that at about 453 s, when Vub is around 1.65 V and Pd is 20 kW, Pb starts to increase, while Pu gradually decreases to cope with the upcoming lower bound of Vub.After about 10 s, Vub reaches the lower bound softly without going beyond it.By this means, AMPC avoids any hard constraints violation.This mechanism is essential to Ah-throughput minimization, as each violation of Vub will cause an extra charging process from LiB to UC, which significantly increases LiB Ah-throughput.

EMSs Results Comparison
Figure 11 shows the AMPC result under UDDS.From Figure 11a, when the HESS power demand is negative (charging), the UC absorbs almost all the braking energy, which is exactly what the AMPC expect.When the demand is positive (discharging), V ub frequently reaches the lower bound.This situation is dealt with by AMPC very softly, for V ub lower bound isn't allowed to be violated.In addition, the AMPC operational mechanism of dealing with hard constraints is shown in Figure 11b,c, which are a zoom from 450 s to 470 s of Figure 11a,d, respectively.We can see from Figure 11b that at about 453 s, when V ub is around 1.65 V and P d is 20 kW, P b starts to increase, while P u gradually decreases to cope with the upcoming lower bound of V ub .After about 10 s, V ub reaches the lower bound softly without going beyond it.By this means, AMPC avoids any hard constraints violation.This mechanism is essential to Ah-throughput minimization, as each violation of V ub will cause an extra charging process from LiB to UC, which significantly increases LiB Ah-throughput.Besides, pulse power in P u is observed when transient power demand comes across, as mentioned in introduction.This phenomenon happens several times in Figure 11b, and would inevitably cause damage to LiB if it's passive, which is exactly what the HESS tries to avoid.Hence, designing LiB as the controlled component is reasonable and necessary.The results of AMPC, DP and RBC are shown in Figure 12, in which 12b and 12c are zooms of 12a. Figure 12b shows ib of three EMSs under a continuous peak power demand which lasts for about 20 s.With the peak demand, optimal ib of DP changes gradually and magnitude is satisfactory, for peak power demand is known in advance.While, ib of AMPC is almost 0 at first to avoid the increase of LiB Ah-throughput, and gradually increases when USOC is about to fall to the lower bound to avoid violation.Though the ib of AMPC is relatively larger, the duration is shortened significantly by QP solver compared with DP and RBC, and the LiB Ah-throughput is reduced through the simulation greatly, indicating excellent flexibility and adaptation of AMPC.LiB current ib of RBC keeps nearly constant at first for Pb keeps at Pb,avg to maintain USOC.At about 200 s, a sudden change in ib appears when USOC reaches the lower bound, as shown in Figure 12e.This indicates RBC is rigid in handling operation mode transitions caused by constraints, and adjustment of USOC by constant Pb,avg is poor.In Figure 12a, ib of DP and AMPC remain nonnegative, while with RBC, LiB charges several times.One of the phenomena is shown in Figure 12c.Referring to Figure 12e, and we can find that the situation appears due to USOC reaches upper bound, again indicating that constant Pb,avg can't adjust USOC flexibly to leave a margin for the braking energy.Consequently, frequent charge obviously increases LiB Ah-throughput, as shown in Figure 12d.As AMPC minimizes ib, energy in UC is always utilized firstly, and USOC always leaves a margin for braking energy, charging to LiB hardly appear.When the power demand is 0, LiB outputs a very small current to UC with DP, which adjusts USOC to deal with peak power demand.Though RBC designs Pb,avg to simulate this action, Pb,avg of RBC can't be real-time optimized based on future demand as well as control target.In Figure 12e, USOC of AMPC fluctuates mainly at the bottom of the range, whereas USOC of RBC fluctuates mainly at the top of the range, again revealing that with AMPC, HESS can take better advantage of UC as a buffer, and absorb more braking energy to reduce LiB Ah-throughput.As a result, when encountering changing driving conditions, AMPC makes the best use of UC as a buffer, while ensuring UC voltage remains within its range, revealing considerable flexibility and adaptivity.Conversely, RBC is rigid in managing bound problems, meanwhile Pb,avg can't be adjusted for realtime, RBC adaptability is poorer and the result is not optimal compared with AMPC.The results of AMPC, DP and RBC are shown in Figure 12, in which 12b and 12c are zooms of 12a. Figure 12b shows i b of three EMSs under a continuous peak power demand which lasts for about 20 s.With the peak demand, optimal i b of DP changes gradually and magnitude is satisfactory, for peak power demand is known in advance.While, i b of AMPC is almost 0 at first to avoid the increase of LiB Ah-throughput, and gradually increases when USOC is about to fall to the lower bound to avoid violation.Though the i b of AMPC is relatively larger, the duration is shortened significantly by QP solver compared with DP and RBC, and the LiB Ah-throughput is reduced through the simulation greatly, indicating excellent flexibility and adaptation of AMPC.LiB current i b of RBC keeps nearly constant at first for P b keeps at P b,avg to maintain USOC.At about 200 s, a sudden change in i b appears when USOC reaches the lower bound, as shown in Figure 12e.This indicates that RBC is rigid in handling operation mode transitions caused by constraints, and adjustment of USOC by constant P b,avg is poor.In Figure 12a, i b of DP and AMPC remain nonnegative, while with RBC, LiB charges several times.One of the phenomena is shown in Figure 12c.Referring to Figure 12e, and we can find that the situation appears due to USOC reaches upper bound, again indicating that constant P b,avg can't adjust USOC flexibly to leave a margin for the braking energy.Consequently, frequent charge obviously increases LiB Ah-throughput, as shown in Figure 12d.As AMPC minimizes i b , energy in UC is always utilized firstly, and USOC always leaves a margin for braking energy, charging to LiB hardly appear.When the power demand is 0, LiB outputs a very small current to UC with DP, which adjusts USOC to deal with peak power demand.Though RBC designs P b,avg to simulate this action, P b,avg of RBC can't be real-time optimized based on future demand as well as control target.In Figure 12e, USOC of AMPC fluctuates mainly at the bottom of the range, whereas USOC of RBC fluctuates mainly at the top of the range, again revealing that with AMPC, HESS can take better advantage of UC as a buffer, and absorb more braking energy to reduce LiB Ah-throughput.As a result, when encountering changing driving conditions, AMPC makes the best use of UC as a buffer, while ensuring UC voltage remains within its range, revealing considerable flexibility and adaptivity.Conversely, RBC is rigid in managing bound problems, meanwhile P b,avg can't be adjusted for real-time, RBC adaptability is poorer and the result is not optimal compared with AMPC.only an average of 6.5%, and educes by an average of 10.6% compared to that of RBC, further revealing the superior optimality as well as adaptivity to driving condition uncertainty of AMPC.

Weights
When i b and ∆i b are to minimize simultaneously, w 2 and w 7 should be adjusted.Fix w 2 at 1, and set w 7 as 1, 5, 10 and 20, simulation results are shown in Figure 13.From Figure 13a we can see that, weight adjustment significantly influences the AMPC control performance.When w 7 increases, i b changes more and more slowly, LiB Ah-throughput increases for this reason and accumulative ∆i b absolute reduces, as shown in Figure 13b,d, while, when w 7 is 20, the trend has been broken.As shown in Figure 13a,c, i b fluctuates severely at 100 s~500 s, ∆i b absolute inevitably increases.As a result, both Ah-throughput and accumulative ∆i b absolute significantly increase during this period, as shown in Figure 13b,d.The accumulative ∆i b absolute when w 7 is 20 is nearly the same with that when w 7 is 5, which is conflicting with the original target of increasing the weight.The situation arises for two reasons: (1) to suppress i b change; (2) to provide energy to UC so that UC could face huge power demand fluctuation.To satisfy two conditions simultaneously, i b absolute has to retain high magnitude and changes drastically between positive and negative.As a result, when w 7 is too big, it's difficult to obtain optimal solution within all constraints.

Weights
When ib and ∆ib are to minimize simultaneously, w2 and w7 should be adjusted.Fix w2 at 1, and set w7 as 1, 5, 10 and 20, simulation results are shown in Figure 13.From Figure 13a we can see that, weight adjustment significantly influences the AMPC control performance.When w7 increases, ib changes more and more slowly, LiB Ah-throughput increases for this reason and accumulative ∆ib absolute reduces, as shown in Figure 13b,d, while, when w7 is 20, the trend has been broken.As shown in Figure 13a,c, ib fluctuates severely at 100 s~500 s, ∆ib absolute inevitably increases.As a result, both Ah-throughput and accumulative ∆ib absolute significantly increase during this period, as shown in Figure 13b,d.The accumulative ∆ib absolute when w7 is 20 is nearly the same with that when w7 is 5, which is conflicting with the original target of increasing the weight.The situation arises for two reasons: (1) to suppress ib change; (2) to provide energy to UC so that UC could face huge power demand fluctuation.To satisfy two conditions simultaneously, ib absolute has to retain high magnitude and changes drastically between positive and negative.As a result, when w7 is too big, it's difficult to obtain optimal solution within all constraints.Define kw = w7/w2, simulation results under different w2 and kw are listed in Tables 7 and 8.No matter how much w2 is, when kw remains unchanged, LiB Ah-throughput and accumulated ∆ib are always the same, indicating that kw is the only decision variable no matter how w2 and w7 change.Define k w = w 7 /w 2 , simulation results under different w 2 and k w are listed in Tables 7 and 8.No how much w 2 is, when k w remains unchanged, LiB Ah-throughput and accumulated ∆i b are always the same, indicating that k w is the only decision variable no matter how w 2 and w 7 change.Furthermore, when k w increases from 0.1 to 20, LiB Ah-throughput increases gradually as well, while accumulated ∆i b first decreases then increases, and the minimum value appears when k w is 10 or 15.As LiB Ah-throughput when k w is 10 is significantly smaller than that when k w is 15, we can conclude that k w should not exceed 10 when both goals are optimized.Based on analysis above, when setting w 2 as 1, and if the ∆i b is to taken into account, w 7 should be set between 5 and 10 (e.g., k w range is between 5 and 10).The proposed AMPC is capable of optimizing multiple objects simultaneously by appropriate weighting, and it's suitable for designing EMS of HESS which usually includes various targets.

Conclusions
In this paper, an AMPC-based EMS for the semi-active HESS is proposed, which minimizes LiB Ah-throughput to extend its life well as ensures HESS operation within the specified ranges.In contrast to standard MPC, the AMPC solves the control action by calculating a QP problem in which prediction model matrices are updated online and model states are estimated by TVKF.In this way, the control action accuracy is significantly improved and LiB Ah-throughput is effectively minimized.The proposed EMS is implemented and verified on the Matlab/Simulink platform, and the simulation results show that compared with SMPC, TVKF reduces the estimation error by 1~3 orders of magnitude, and AMPC reduces LiB Ah-throughput by 4.3% under UDDS.Furthermore, all HESS states are operating within specified ranges, indicating superior model adaptivity.The LiB Ah-throughput of AMPC under six classical driving cycles differs from that of DP by only an average of 6.5%, and is reduced by an average of 10.6% compared to that of RBC, revealing excellent optimality and adaptation to driving condition uncertainty.This method is ideal for managing nonlinear optimization problems with unpredictable disturbances as well as engineering constraints, and can be extended to EMSs for other nonlinear multi-source systems of new energy vehicles.

Table 5 .
Average states estimation errors of SSKF and TVKF.and CA elements of AMPC under UDDS are shown in Figure7, in which elements with marked change are shown in bold.It's obvious that elements variation in three matrices are nonnegligible, further indicating that a single model is unable to approximate the real HESS, and SSKF is unable to make accurate estimation.

Figure 7 .
Figure 7. GA, HA and CA variation of AMPC.

G
A , H A and C A elements of AMPC under UDDS are shown in Figure7, in which elements with marked change are shown in bold.It's obvious that elements variation in three matrices are non-negligible, further indicating that a single model is unable to approximate the real HESS, and SSKF is unable to make accurate estimation.

Table 5 .
Average states estimation errors of SSKF and TVKF.and CA elements of AMPC under UDDS are shown in Figure7, in which elements with marked change are shown in bold.It's obvious that elements variation in three matrices are nonnegligible, further indicating that a single model is unable to approximate the real HESS, and SSKF is unable to make accurate estimation.

Figure 7 .
Figure 7. GA, HA and CA variation of AMPC.

Figure 7 .
Figure 7. G A , H A and C A variation of AMPC.

Figure 8 .
Figure 8.Simulation results comparison of AMPC and SMPC.(a) LiB current comparison.(b) LiB current comparison between 320 s and 390 s; (c) LiB current comparison between 1170 s and 1220 s; (d) LiB Ah-throughput comparison.

Figure 8 .
Figure 8. Simulation results comparison of AMPC and SMPC.(a) LiB current comparison.(b) LiB current comparison between 320 s and 390 s; (c) LiB current comparison between 1170 s and 1220 s; (d) LiB Ah-throughput comparison.
power in Pu is observed when transient power demand comes across, as mentioned in introduction.This phenomenon happens several times in Figure11b, and would inevitably cause damage to LiB if it's passive, which is exactly what the HESS tries to avoid.Hence, designing LiB as the controlled component is reasonable and necessary.

Table 3 .
Influence of sample time and prediction duration.

Table 3 .
Influence of sample time and prediction duration.

Table 5 .
Average states estimation errors of SSKF and TVKF.

Table 7 .
LiB Ah-throughput under different kw and w2.

Table 8 .
LiB accumulated ∆ib under different kw and w2.

Table 7 .
LiB Ah-throughput under different k w and w 2.

Table 8 .
LiB accumulated ∆i b under different k w and w 2.