Energy Management in Plug-in Hybrid Electric Vehicles: Convex Optimization Algorithms for Model Predictive Control

TL;DR

This paper compares convex optimization algorithms for energy management in hybrid electric vehicles, demonstrating that the projected interior point method outperforms ADMM and CVX in accuracy and speed.

Contribution

It introduces a projected interior point method tailored for energy management optimization, improving computational efficiency for high-accuracy solutions.

Findings

Projected interior point method is faster and more accurate than ADMM and CVX.

ADMM is preferable for solutions with modest accuracy.

Both proposed methods outperform general-purpose convex optimization software.

Abstract

This paper details an investigation into the computational performance of algorithms used for solving a convex formulation of the optimization problem associated with model predictive control for energy management in hybrid electric vehicles with nonlinear losses. A projected interior point method is proposed, where the size and complexity of the Newton step matrix inversion is reduced by applying inequality constraints on the control input as a projection, and its properties are demonstrated through simulation in comparison with an alternating direction method of multipliers (ADMM) algorithm, and general purpose convex optimization software CVX. It is found that the ADMM algorithm has favourable properties when a solution with modest accuracy is required, whereas the projected interior point method is favourable when high accuracy is required, and that both are significantly faster than CVX.

Figures

Tables

Table 1. TABLE I: Summary of the matrix/vector operations present in each equation used in the projected interior point algorithm, where v refers to a vector, and M refers to a non-diagonal matrix.

Equation(s)	M-1	M$\cdot$M	M$\cdot$v	$𝒪\left(2^n{N}^n\right)$
(23a)	Yes	Yes	Yes	$n\le 3$
(23b-23c)	No	No	Yes	$n\le 2$
(24-25)	No	No	No	$n=1$
(26)	No	No	Yes	$n\le 2$

Table 2. TABLE II: Summary of the matrix/vecotr present in the ADMM variable updates, where ’v’ refers to a vector, and ’M’ refers to a non-diagonal matrix other than Ψ Ψ \Psi .

Update	M-1	M$\cdot$M	M$\cdot$v	$𝒪\left(N^n\right)$
$u$	No	No	No	$n=1$
$x$	No	No	No	$n=1$
$\zeta$	No	No	Yes	$n\le 2$
${\lambda }_1,{\lambda }_2$	No	No	No	$n=1$
$r_{\text{P}},r_{\text{D}}$	No	No	No	$n=1$

Equations

$${\dot{m}}_f(t)=f(P_{eng}(t),{\omega }_{eng}(t),\sigma (t),t),$$

$$\sigma (t)=\{\begin{array}{ll}1& \text{engine on, clutch engaged}\\ 0& \text{engine off, clutch disengaged}\end{array}.$$

$$P_b(t)$$

$$P_c(t)$$

$$E(t)=E(0)-{\int }_0^t{P}_b(t)\mathrm{d}t.$$

$$P_{drv}(t)=P_{em}(t)+P_{eng}(t)+P_{brk}(t),$$

$${\omega }_{em}(t)=r(t){\omega }_w(t),{\omega }_{eng}(t)=\sigma (t)r(t){\omega }_w(t)$$

$$P_{eng}(t)\ge {\underline{P}}_{eng}(t)$$

$$P_{eng}(t)\le {\overline{P}}_{eng}(t)$$

$$P_{em}(t)\ge {\underline{P}}_{em}(t)$$

$$P_{em}(t)\le {\overline{P}}_{em}(t)$$

$${\underline{\omega }}_{em}\le {\omega }_{em}(t)\le {\overline{\omega }}_{em}\text{and}{\underline{\omega }}_{eng}\le {\omega }_{eng}(t)\le {\overline{\omega }}_{eng},$$

$$\underline{E}\le E(t)\le \overline{E}\text{and}{\underline{P}}_b\le P_b(t)\le {\overline{P}}_b.$$

$$\underset{\sigma (t),r(t),P_{eng}(t),P_{brk}(t)}{\min }$$

$$(\text{ref}eq{n}_{f}ue{l}_{c}onsumption)-(\text{ref}eq{n}_{b}atter{y}_{l}imits)\forall t\in [0,T]$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{drv}}$$

$$\displaystyle\hat{\dot{m}}_{f,k}=\hat{f}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}},\hat{\omega}_{eng,k},\hat{\sigma}_{k})$$

$$=$$

$$\displaystyle\hat{P}_{c,k}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{h}}(\hat{P}_{em,k},\hat{\omega}_{em,k})$$

$$=$$

$${\hat{P}}_{b,k}=$$

$$=$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{eng,k}}$$

$${\overline{P}}_{eng,k}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{em,k}}$$

$${\overline{P}}_{em,k}$$

$$\displaystyle\min_{\hat{\sigma},\hat{r},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng}},\hat{P}_{brk}}$$

$${\hat{E}}_0=E(t)$$

$$\displaystyle\begin{rcases}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{E}}_{k+1}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{E}_{k}}-\delta\hat{g}(\hat{P}_{em,k},\hat{\omega}_{em,k},k)\\ {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{drv,k}}=\hat{P}_{em,k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}}+\hat{P}_{brk,k}\\ \hat{\omega}_{em,k}=\hat{r}_{k}\hat{\omega}_{w,k}\\ \hat{\omega}_{eng,k}=\hat{\sigma}_{k}\hat{r}_{k}\hat{\omega}_{w,k}\\ \hat{r}_{k}\in\{r_{1},\dots,r_{N_{g}}\}\\ \hat{\sigma}_{k}\in\{0,1\}\\ {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{eng,k}}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}}\leq\overline{P}_{eng,k}\\ {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{em,k}}\leq\hat{P}_{em,k}\leq\overline{P}_{em,k}\\ \underline{\omega}_{em}\leq\hat{\omega}_{em,k}\leq\overline{\omega}_{em}\\ \underline{\omega}_{eng}\leq\hat{\omega}_{eng,k}\leq\overline{\omega}_{eng}\\ \underline{P}_{b}\leq\hat{P}_{b,k}\leq\overline{P}_{b}\\ \underline{E}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{E}}_{k+1}\leq\overline{E}\\ \end{rcases}\forall k$$

$$\mathcal{P}=\{k:\hat{\sigma }(k)=1\},$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}}_{drv,k}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}}_{em,k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}}_{eng,k},$$

$${\hat{\dot{m}}}_f(k)$$

$$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{2,k}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}}^{2}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{1,k}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{0,k}},$$

$${\hat{P}}_c(k)$$

$$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{2,k}}\hat{P}_{em,k}^{2}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{1,k}}\hat{P}_{em,k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{0,k}}.$$

$$\displaystyle\hat{P}_{b,k}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}}(\hat{P}_{em,k})=\frac{V^{2}_{oc}}{2R}\left(1-\sqrt{1-\frac{4R}{V_{oc}^{2}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{h}}_{k}}(\hat{P}_{em,k})}\right).$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{eng,k}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{P}_{em,k}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{P}_{em,k}}$$

Full text

I Introduction

Increased electricification of road vehicles has been identified as a key short term solution to important societal issues including climate change and air pollution [1]. Plug-in hybrid electic vehicles (PHEVs), where an electric propulsion system is complemented with an internal combustion engine, are currently a common configuration. Although the low energy density and lengthy recharge time of lithium ion batteries limits the viability of all-electric powertrains, analysis of daily driving behaviour reveals that 50% of internal combustion powered miles can be powered electrically using a hybrid vehicle with an all-electric range of just 40 miles [2]. The inclusion of an additional power source, however, introduces a challenging problem: at each instant during a given journey, how much power should be delivered from the motor, and how much should be delivered from the engine?

This is known as the energy management problem [3], and a simple heuristic is a charge depleting/charge sustaining strategy, where power is delivered from only the electric motor until the battery is sufficiently depleted, and then the vehicle is operated in a charge sustaining mode until the end of the journey [2]. It has, however, been demonstrated that significant savings in fuel consumption can be made by delivering power from both the motor and engine simultaneously, and modulating the fraction delivered from each throughout the journey in what is known as a ‘blended mode’ [4]. There are several methods for controlling the powertrain in this way, and the globally optimal solution can be obtained for complex, nonlinear models and/or integer conrol decisions (such as gear selection) using Dynamic Programming [5, 6], but this approach is too computationally demanding for an online solution. Other studies have investigated methods based on Pontryagin’s Minimum Principle [7, 8, 9], although it is challenging to enforce complex constraints, such as the general state of charge constraint or engine switching, whilst still guaranteeing optimality.

Model predictive control (MPC) has shown promise in this application due to the inherent robustness to uncertainty in both the vehicle model and prediction of future driving behaviour [10], and nonlinear models of losses can be used throughout the hybrid powertrain for improved performance [11]. The associated online optimization problem is still computationally intractable when gear selection and engine switching are considered, so these elements are commonly removed from the problem or optimized externally [11, 12, 13, 14], and it has been demonstrated that the power balance alone can be formulated, with nonlinear losses and without simplification, as a convex optimization problem with linear state dynamics [15].

For the last 30 years, the most popular algorithms for solving inequality constrained optimization problems have been interior point methods [16]. Originally formulated as ‘primal’ methods by approximating inequality constraints in an optimization problem with logarithmic barrier functions, their inherent ill-conditioning and numerical inefficiency rendered interior point methods ineffective until the publication of Karmakar’s method [17] in 1984. The subsequently developed ‘primal-dual’ interior point methods [18] displayed excellent theoretical and practical properties, including polynomial complexity and a near constant number of iterations with variations in problem size, and today, a large volume of research output in the field of optimization for MPC is dedicated to the development and application of primal-dual interior point methods [19, 20, 21].

In the recently published literature, algorithms for solving convex formulations of the energy management problem have not been investigated and the optimization problem has normally been solved using general purpose convex optimization software [13, 22, 23, 14], with the exception of the alternating direction method of multipliers (ADMM) algorithm presented in [15]. The first contribution of this paper is a projected interior point solver for this problem, where the size of the matrix inversion associated with the primal-dual Newton step is reduced by enforcing the element-wise inequality constraints on the decision variable as a projection, thereby reducing the computational requirement of each iteration. The algorithm is not domain specific, and is applicable to any MPC optimization problem with linear dynamics, a separable convex cost function of the control variable, and upper and lower bounds on the control and state variables.

A primary motivation of this paper was to determine the relative computational benefits of second-order and first-order methods for the convex PHEV energy management formulation, and the second contribution is a set of numerical studies where the performance of the projected interior point algorithm is compared with the ADMM algorithm of [15]. In these studies we demonstrate that the projected interior point method has superior convergence properties, but requires more time to obtain a solution with modest accuracy, and consequently is only suitable for real-time solutions over shorter horizons (in this case fewer than 500 samples). We also demonstrate that both methods are significantly faster than CVX [24], and in ADMM (using an improved implementation from [15]) we demonstrate the first method capable of solving the energy management problem in real time, over long horizons ($\geq$1000 samples) when nonlinear system dynamics are considered and hard limits on both power and state of charge are enforced over the entire horizon.

The paper is organised as follows: in section II the energy management problem, MPC framework, and convex reformulation are defined, and section III details the projected interior point method. The ADMM algorithm of [15] is stated in section IV, numerical experiments are presented in section V, and the paper is concluded in section VI.

II Energy Management Problem & Model Predictive Control Framework

Fig. 1 shows a simplified diagram of a parallel PHEV powertrain, and illustrates the energy transfers that are considered as part of the energy management problem. At a given time, $t$, the mass flow rate of fuel delivered to the engine, $\dot{m}_{f}$, can be described by a time-varying function, $f$, of engine output power, $P_{eng}$, and engine shaft speed, $\omega_{eng}$, as

[TABLE]

where $\sigma$ describes the state of the engine and clutch engagement by

[TABLE]

Similarly, the rate of consumption of the battery’s internal chemical energy, $P_{b}$, can be described by a time varying function, $g$, of battery output power (i.e motor input power), $P_{c}$, which can in turn be described by a time varying function, $h$, of motor output power, $P_{em}$, and motor shaft speed, $\omega_{em}$:

[TABLE]

Therefore, the state of charge of the battery, $E$, is given at time $t$ by

[TABLE]

The engine output power that is delivered through the clutch is combined additively with power from the motor through a coupling device to drive the gearbox. Assuming that all drivetrain components are $100\%$ mechanically efficient, the power delivered to the wheels (i.e the power demanded by the driver), $P_{drv}$, is given by

[TABLE]

where $P_{brk}$ is the power extracted from the system by the mechanical brakes. Assuming a discrete variable transmission, the rotational velocities of the engine and motor shafts are given as a function of the rotational velocity of the wheels, $\omega_{w}$, as

[TABLE]

where $r:\mathbb{R}\rightarrow\{r_{1},\dots,r_{N_{g}}\}$, and $N_{g}$ is the number of available gear ratios.

The engine has upper and lower limits on torque, $\overline{T}_{eng}$ and $\underline{T}_{eng}$, that are functions of engine speed, so limits on engine power are given by

[TABLE]

The limits on motor power can be given similarly as

[TABLE]

where $\overline{T}_{em}$ and $\underline{T}_{em}$ are upper and lower limits on motor torque. The engine and motor have static limits on rotational speed, given by

[TABLE]

and the battery has static limits on state-of-charge and rate of charge and discharge, given by

[TABLE]

The above system is under-constrained in three degrees of freedom: the fraction of total driver demand power delivered from motor, engine, and brakes; the engine switching and clutch engagement; and the gear selection. Consequently, the parameters $P_{eng}(t)$, $P_{brk}(t)$, $\sigma(t)$, and $r(t)$ must be actively controlled ($P_{em}(t)$ is given in terms of $P_{eng}(t)$ and $P_{brk}(t)$ by (5)). The energy management problem can therefore be written as an open-loop optimal control problem for a journey of length $T$ as

[TABLE]

Note that $\omega_{w}$ and $P_{drv}$ (which determine the vehicle’s speed and acceleration) are not affected in the control problem; the principle of the controller is to always meet the powertrain output demanded by the driver, and to not affect the overall driving behaviour of the vehicle.

II-A MPC Framework

If implemented in a real vehicle, the solution found from (11) will be suboptimal as it is impossible to model the powertrain components with complete accuracy, and because the problem is dependent on future disturbance variables, $P_{drv}$ and $\omega_{w}$, that are impossible to exactly predict a priori. A MPC framework can be used to reduce these limitations, where instead of solving a single instance of the open-loop control problem, the control variables are repeatedly updated as the journey progresses [11]. This allows the predictions of driver behaviour to be improved as new information becomes available, and provides feedback on the vehicle state (i.e the battery state of charge), thus providing a degree of robustness to modelling and prediction errors. We describe a MPC framework for the energy management problem in this section.

Throughout the following text, the notation $\hat{x}$ is used for a variable, $x(t)$, to refer to the discretely sampled prediction used within the MPC framework, as opposed to the physical signals and states described in the previous section. At each control variable update instant, a discretely sampled prediction of demand power and wheel speed is made as

[TABLE]

for $k=0,\dots,N-1$, where the sampling period $\delta$ is assumed to be constant, and the prediction horizon is given by $N$. Although the prediction aspect of the energy management problem is still very much an open issue [25], the focus of this paper is on the subsequent optimization problem, so it is assumed that an accurate method is available to the controller and this aspect is not addressed further. The engine and motor loss maps can be approximated with quasi-static quadratic functions, $\hat{f}$ and $\hat{h}$, [13, 11, 14] as

[TABLE]

where $\alpha_{2}(\hat{\omega}_{eng,k})>0$ $\forall$ $\hat{\omega}_{eng,k}$, and $\beta_{2}(\hat{\omega}_{em,k})>0$ $\forall$ $\hat{\omega}_{em,k}$. The battery is commonly modelled as an equivalent circuit of internal resistance [26, 27, 13, 11, 15, 14, 12] as:

[TABLE]

where $V_{oc,k}$ and $R_{k}$ are the open circuit voltage and internal resistance. The limits on engine, motor, and discharge power (7,8) are now given in discrete time as

[TABLE]

and the MPC optimization at time $t$ is then given as an approximation of (11) with Euler method integration of the state dynamics (4) as

[TABLE]

At each control variable update instance, the first elements of the vectors of optimization variables, $\hat{\sigma}^{\star}_{0}$, $\hat{r}^{\star}_{0}$, $\hat{P}_{eng,0}^{\star}$, and $\hat{P}_{brk,0}^{\star}$ are implemented as $\sigma(t)$, $r(t)$, $P_{eng}(t)$, and $P_{brk}(t)$, where $\hat{\sigma}^{\star}$, $\hat{r}^{\star}$, $\hat{P}_{eng}^{\star}$, and $\hat{P}_{brk}^{\star}$ are the minimizing arguments of (13).

II-B Convex Reformulation

Whilst the MPC framework provides a degree of robustness to prediction and modelling errors, problem (13) is challenging to solve as the cost function is non-convex, it is subject to nonlinear, non-convex constraints, and it has $2N$ discrete decision variables ($N$ is likely to be in the thousands for journeys longer than 15 mins, assuming an update frequency of approximately 1Hz). The global minimum of the problem can be found using Dynamic Programming, however this is a computationally demanding approach and not suitable for an online solution [15]. Instead, it is possible to reduce the complexity of the problem by determining the three variables $\hat{\sigma}$, $\hat{r}$, and $\hat{P}_{brk}$ using heuristic rules [11] or an external optimization routine [13], leaving the power balance between the engine and motor as the only under-constrained parameter. It has previously been demonstrated in [15] that the resulting problem is convex when using the battery power, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{b}}, as the decision variable. This is the approach that is taken in this paper, as discussed below.

Assuming an external method for estimating the engine switching behaviour, the set of timesteps at which the engine is switched on and the clutch is engaged is given by

[TABLE]

and also assuming a method for pre-determining the use of mechanical brake, we re-define the drive demand power from (5) as

[TABLE]

so that $P_{brk,k}$ is no longer an optimization variable in (13). Finally, by assuming that the gear selection, $\hat{r}$, is also pre-determined, $\hat{\omega}_{eng,k}$ and $\hat{\omega}_{em,k}$ are defined for all $k$, so the engine and motor models (12) can be reduced to time-varying polynomial functions of the form

[TABLE]

These functions are strictly convex ({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{2,k}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{2,k}}>0), and we then ensure that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}} is also strictly convex by assuming that $V_{oc}$ and $R$ are independent of state of charge (and in this case, constant), so that the battery model is approximated by

[TABLE]

We update the limits on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{eng,k}} and $\hat{P}_{em,k}$ to ensure that $\hat{f}_{k}$, $\hat{h}_{k}$, and $\hat{g}_{k}$ are all non-decreasing and real-valued as

[TABLE]

This also ensures that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}} is a one-to-one function, so we can define

[TABLE]

where

[TABLE]

Using this definition of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}^{-1}}, it is known that for $k\in\mathcal{P}$,

[TABLE]

where the corresponding limits on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{b}} are

[TABLE]

and it is known that for $k\notin\mathcal{P}$

[TABLE]

In [15], it is shown that given the properties of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}}(\cdot) (strictly convex, twice differentiable, non-decreasing, one-to-one) and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{f}_{k}} (convex and non-decreasing), the function {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{f}_{k}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{drv,k}}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{g}_{k}^{-1}}(\hat{P}_{b,k})) is convex, non-increasing, and twice differentiable, so the MPC problem (13) becomes the convex, linearly constrained optimization problem

[TABLE]

For the sake of clarity in the following sections, we now revert to the commonly used notation for MPC problems, where $u$ is the control input (the predicted battery power, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{b}}), and $x$ is the state variable (the predicted state of charge, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{E}}), i.e

[TABLE]

Then, by defining the non-increasing, separable, strictly convex cost function

[TABLE]

we can express (14) equivalently as

[TABLE]

where $\Phi$ is a vector of $N$ ones, and $\Psi$ is an $N\times N$ lower triangular matrix where every non-zero element is equal to $\delta$.

If we describe $\mathcal{F}_{k}$ as the feasible set of state of charge values at timestep $k$, then problem (14) is feasible if and only if $\mathcal{F}_{k}\neq\emptyset$ for all $k=0,\dots,N$, where $\mathcal{F}_{k}$ is defined recursively by

[TABLE]

and

[TABLE]

As $\mathcal{F}_{k}$ is a one-dimensional convex set, it can be parameterized by $\mathcal{F}_{k}=[\min\mathcal{F}_{k},\max\mathcal{F}_{k}]$, and the sequence $\mathcal{F}_{1},\dots,\mathcal{F}_{N}$ can be obtained iteratively from

[TABLE]

with $\mathcal{F}_{0}=\{x_{0}\}$. The problem is therefore feasible if and only if $\max\mathcal{F}_{k+1}$ and $\min\mathcal{F}_{k+1}$ exist (i.e. $\max\mathcal{F}_{k+1}\geq\min\mathcal{F}_{k+1}$ when calculated from (16), and $\overline{u}_{k}\geq\underline{u}_{k}$) for $k\in\{0,\dots,N-1\}$.

The cost function in (15) is known to be non-increasing in $u$, so by inspection, if

[TABLE]

then $\overline{u}$ is the minimizing argument.

III Projected Interior Point Method

Problem (15) is in the form of a convex nonlinear program with affine equality and inequality constraints. Here we present a projected primal-dual interior point algorithm, where the element-wise bounds on the control variable are applied as a projection. This section begins with the definition of the barrier approximation and optimality conditions, followed by a statement of the initialization algorithm and main projected interior point algorithm with accompanying pseudocode. The section is then is then concluded with an analysis of the computational complexity of each iteration, and a discussion of convergence to the minimizing argument of (15).

Firstly, we introduce a slack variable $s\in\mathbb{R}^{2N}$, and write problem (15) equivalently as

[TABLE]

where

[TABLE]

This can then be approximated with

[TABLE]

where $B(s,\mu)$ is a log barrier function defined by

[TABLE]

and $\mu>0$ can be interpreted as the degree to which the log barrier function approximates the inequality constraint, $s\geq 0$. The Lagrangian function associated with problem (18) is

[TABLE]

where $\theta_{1}\in\mathbb{R}^{2N}$, $\theta_{2}\in\mathbb{R}^{N}$, and $\theta_{3}\in\mathbb{R}^{N}$. By defining the set

[TABLE]

the necessary and sufficient conditions for the optimal values $u^{\circ},s^{\circ}$, and $\theta_{1}^{\circ}$ that minimize (18) are

[TABLE]

where $\nabla_{k}F(u)$ is the $k$th row of the gradient of $F$ at $u$, $A_{k}$ is the $k$th column of $A$, and $S=\text{diag}(s)$.

It can be demonstrated that to obtain the optimality conditions of (17), two changes must be made to (19): firstly, (19f) must become a non-strict inequality condition, and secondly, a vector of zeros must replace $\frac{1}{\mu}\textbf{1}$ in the R.H.S of (19b). Therefore, it can be concluded that $u^{\circ}$ converges asymptotically to $u^{\star}$ as $\mu\to\infty$, where $u^{\star}$ is the minimizing argument of (17). The principle of the algorithm presented in this section is that conditions (19d-19g) hold at all iterations, whilst a projected Newton method [28] is used to obtain an approximation of the solutions of (19a-19c) for a fixed value of $\mu$. This is then repeated for progressively larger values of $\mu$, with progressively higher accuracy, to obtain $u^{\star}$.

III-A Initialization

The projected interior point algorithm is initialized using Algorithm 1 to obtain the tube defined by the sequence $\mathcal{G}_{0},\dots,\mathcal{G}_{N}$ where

[TABLE]

The centerline of this tube is then used to obtain values for $u^{(0)}$, $s^{(0)}$, and $\theta_{1}^{(0)}$ that satisfy (19b-19g).

Proposition III.1**.**

The values of $u^{(0)}$, $s^{(0)}$, and $\theta_{1}^{(0)}$ obtained with Algorithm 1 will satisfy conditions (19b-19g), iff

[TABLE]

for $k\in{1,\dots,N}$.

Proof.

By induction. We start by defining the operation

[TABLE]

for any given set $\mathcal{S}$, and note that for $\max\mathcal{F}_{k}$ and $\min\mathcal{F}_{k}$ to exist in (20) it is implied that the problem is feasible, so $|\mathcal{F}_{k}|\geq 0$ for $k\in\{1,\dots,N\}$. Let $k^{\dagger}$ be the smallest value of $k$ such that $\max\mathcal{F}_{k}=\max\mathcal{X}$ and/or $\min\mathcal{F}_{k}=\min\mathcal{X}$, and suppose that $k^{\dagger}\leq N$. In this case we can use (20) to show that $\max\mathcal{F}_{k}>\min\mathcal{F}_{k}$ (i.e. $|\mathcal{F}_{k}|>0$). Furthermore, we know that $|\mathcal{F}_{k+1}|\geq|\mathcal{F}_{k}|$ for any $k$ where $\max\mathcal{F}_{k+1}\neq\max\mathcal{X}$, $\min\mathcal{F}_{k+1}\neq\min\mathcal{X}$, and $|\mathcal{U}_{k}|\geq 0$. Therefore $|\mathcal{F}_{k}|>0$ for all $k\geq k^{\dagger}$. If we assume that $|\mathcal{F}_{k}|>0$, and $|\mathcal{G}_{k+1}|>0$, we can then show that

[TABLE]

In the case where $k^{\dagger}$ exists we know that $|\mathcal{F}_{N}|>0$, so $|\mathcal{G}_{N}|>0$, and the argument in (21) can then be made recursively to show that $|\mathcal{G}_{k}|>0$ for $k=N-1,\dots,k^{\dagger}$. As $\mathcal{G}_{k}\subseteq\mathcal{F}_{k}\subseteq\mathcal{X}$, we now know that that $x_{k}^{(0)}\in\text{int}(\mathcal{X})$ for $k\in\{k^{\dagger},\dots,N\}$.

For $k\in\{1,\dots,k^{\dagger}-1\}$ we know that $|\mathcal{F}_{k}|\geq 0$ from feasibility, and using a similar method to (21) we can then show that $|\mathcal{F}_{k+1}|\geq 0\Rightarrow|\mathcal{G}_{k}|\geq 0$. We also know that $\max\mathcal{F}_{k}\neq\max\mathcal{X}$ and $\min\mathcal{F}_{k}\neq\min\mathcal{X}$ for $k\in\{1,\dots,k^{\dagger}-1\}$ (from the definition of $k^{\dagger}$), so as $\mathcal{G}_{k}\subseteq\mathcal{F}_{k}\subset\mathcal{X}$, it follows that $x_{k}^{(0)}\in\text{int}(\mathcal{X})$ for $k\in\{1,\dots,k^{\dagger}-1\}$. A similar result can be shown for $k\in\{1,\dots,N\}$ if $k^{\dagger}$ does not exist.

The condition $x_{k}^{(0)}\in\text{int}(\mathcal{X})$ for $k\in\{1,\dots,N\}$ ensures that $Au^{(0)}-b>0$, so $s^{(0)}=Au^{(0)}-b$ ensures (19c) and (19f), and $\theta_{1}^{(0)}=\frac{1}{\mu}(S^{(0)})^{-1}\textbf{1}$ ensures (19b) and (19g). Finally, it can also be shown that

[TABLE]

and we can similarly show that $\min u_{k}^{(0)}=\underline{u}_{k}$, which therefore demonstrates (19d-19e) if (20) is true for $k\in\{1,\dots,N\}$.

Conversely, a subset of the above results are shown to be not true if (20) is not true for some $k\in\{1,\dots,N\}$. ∎

An example of the $x^{(0)}$ and $u^{(0)}$ values obtained by the initialization algorithm for a nominal, randomly generated example is shown in Fig. 2, demonstrating that $\underline{u}_{k}\leq u_{k}^{(0)}\leq\overline{u}_{k}\ \forall k$, and that $\underline{x}<x^{(0)}_{k}<\overline{x}\ \forall k$.

In addition to $u$, $s$, and $\theta_{1}$, there are three further parameters that are initialized at the start of the algorithm: $\mu_{0}>0$, $\overline{\mu}\geq\mu_{0}$, $k_{\mu}>1$, and $\tau\in(0,1)$. These parameters can be assigned any value within the stated ranges, and their significance is discussed in the following subsection.

III-B Algorithm

At each iteration, $j$, the elements $k$ are partitioned into the sets

[TABLE]

then an estimate is made of the solution to the equations (19a-19c) using a projected Newton method. Let

[TABLE]

where $\nabla^{2}_{\mathcal{D}}F(u^{(j)})$ is the $k\in{\mathcal{D}^{(j)}}$ rows and columns of the Hessian of $F$ evaluated at $u^{(j)}$, $\nabla_{\mathcal{D}}F(u^{(j)})$ is the $k\in{\mathcal{D}^{(j)}}$ elements of the gradient of $F$ evaluated at $u^{(j)}$, and $\Theta=\text{diag}(\theta)$. We define $\hat{A}$ as the matrix $A$ with the $k\in\{k:\underline{u}_{k}=u^{(j)}_{k}\ \text{or}\ \overline{u}_{k}=u^{(j)}_{k}\}$ columns set to zero, so that $v$ is the $k\in\mathcal{D}^{(j)}$ columns of $A$, and $w$ is the $k\in\mathcal{D}^{(j)}$ columns of $\hat{A}$. The search directions $\Delta\tilde{u}^{(j)}$, $\Delta\tilde{s}^{(j)}$, and $\Delta\theta_{1}^{(j)}$ are obtained from the reduced equations

[TABLE]

and search lengths $\alpha_{s}$ and $\alpha_{\theta}$ are the determined from the ‘fraction to the boundary’ rule, [29, pp. 567] as

[TABLE]

where $\tau\in(0,1)$ is fixed and arbitrary. The variables $u$, $s$, and $\theta_{1}$ are then updated as

[TABLE]

where $\pi_{k}^{u}(u_{k})=\min\{\overline{u}_{k},\max\{\underline{u}_{k},u\}\}$, and $i(k)$ is a function that returns the index of the element of $\mathcal{D}^{(j)}$ that is equal to $k$, assuming that the elements are ordered in a chronological, increasing sequence (e.g $i(7)=2$ if $\mathcal{D}^{(j)}=\{1,7,12,\dots\}$). This iteration is performed repeatedly until the criterion

[TABLE]

is met, at which point the value of $\mu$ is updated using

[TABLE]

where $\overline{\mu}>0$ is a pre-determined upper limit on the value of $\mu$, and $k_{\mu}>1$ is a pre-determined, arbitrary constant. The algorithm as a whole is then terminated when both conditions (26) and $\mu=\overline{\mu}$ are met. Algorithm 2 presents a pseudocode implementation of the above description.

III-C Complexity

Lines 6-13 in Algorithm 2 constitute the recursive elements of the projected interior point algorithm, and Table I presents an analysis of the complexity of equations (23-26). The significance of enforcing the bounds on $u$ as a projection is illustrated, as the complexity of the projected interior point operations is a function of the first dimension of $A$, which is $2N$ here. If the bounds on $u$ were applied as log barrier functions (i.e $A=(\Psi,-\Psi,I,-I)$ and $b=(\Phi(x_{0}-\overline{x}),-\Phi(x_{0}-\underline{x}),\underline{u},-\overline{u})$ in (17)) the relevant dimension would instead be $4N$, and the complexity of each update would become $\mathcal{O}(4^{n}N^{n})$. It can be seen that (23a) is the most computationally demanding update due to the presence of both a dense matrix-matrix multiplication and a dense matrix inverse (note that diagonal matrix operations, e.g. $(\nabla^{2}_{\mathcal{D}}F(u^{(j)}))^{-1}$, are omitted from Table I). The computational complexity of each iteration of the projected interior point algorithm is therefore $\mathcal{O}(2^{n}N^{n})$, where $n\leq 3$ is determined by the method used for matrix multiplication and inversion.

III-D Convergence

The algorithm presented in Section III-B can be interpreted as a projected Newton method [28] used to obtain a stationary point, ($u^{\circ}$, $s^{\circ}$, $\theta_{1}^{\circ}$), of the function

[TABLE]

for a given value of $\mu$, subject to the constraint $\underline{u}\leq u\leq\overline{u}$. The strict inequality in the definition of $\mathcal{A}^{\circ}$ means that there is a region of $(u,s,\theta_{1})$-space close to $(u^{\circ},s^{\circ},\theta_{1}^{\circ})$ where $\mathcal{A}^{(j)}=\mathcal{A}^{\circ}$, and a subset of this region will meet the conditions for local quadratic convergence of Newton’s method for nonlinear equations [29, pp. 276]. This means that the R.H.S of (22) converges to 0 as $j\to\infty$, and the termination criterion (26) will be met in a finite number of steps. Global convergence could be ensured by adapting the $\Delta u$ step at each iteration with a line-search [29, pp. 30] of an appropriate merit function, although the merit function from [28] cannot be used as the stationary point ($u^{\circ}$, $s^{\circ}$, $\theta_{1}^{\circ}$) is not a minimum of the function (28) in general. Despite this limitation, convergence was demonstrated for all problem classes in the simulations that follow.

We have previously demonstrated that $u^{\circ}\to u^{\star}$ as $\mu$ is increased towards $\infty$, and the value of $u^{(j)}$ when the criterion (26) is met converges to $u^{\circ}$ as $\mu$ is increased towards $\infty$. Therefore, the value of $u^{(j)}$ when Algorithm 2 terminates can be made arbitrarily close to the minimizing argument of (17) by setting $\overline{\mu}$ arbitrarily high. The algorithm could be further optimised to update the value of $\mu$, possibly at every iteration, to ensure that the iterate remains in, or at least near, to the superlinearly convergent region around $u^{\circ}$, $s^{\circ}$, and $\theta_{1}^{\circ}$, although in the numerical experiments that follow we demonstrate superlinear convergence for a broad class of problems using the simple update (27).

IV Alternating Direction Method of Multipliers

We compare the performance of the projected interior point method in simulation with the ADMM algorithm proposed in [15], which is restated here with a new complexity analysis. We introduce a dummy variable, $\zeta$, and rewrite (15) as

[TABLE]

where the indicator function, $\Lambda$, is defined by

[TABLE]

and the augmented Lagrangian associated with (29) is

[TABLE]

The ADMM algorithm is initialized with the values

[TABLE]

and by defining projection functions

[TABLE]

the iteration is given by

[TABLE]

Problem (29) can be shown to be equivalent to the canonical ADMM form [30, equation (3.1)], for which the iteration is equivalent to (31). As $F(u)+\Lambda(u,x)$ is convex, it can therefore be concluded that iteration (31) will converge to the solution of (29) since the residuals defined by

[TABLE]

necessarily converge to zero. The algorithm is terminated when the conditions $\|r_{\text{P}}^{(j+1)}\|\leq\epsilon$ and $\|r_{\text{D}}^{(j+1)}\|\leq\epsilon$ are met, where $\epsilon$ is a pre-determined threshold.

IV-A Complexity

Algorithm 3 shows a pseudocode implementation of the ADMM algorithm, and the computational complexity of each recursive variable update is presented in Table IV-A. Each $u_{k}$ update is an unconstrained convex optimization problem that we solve here using a Newton method with a backtracking line search, so the $u$ update therefore scales linearly with $N$ if these updates are performed sequentially, or is constant if each $k$ update can be performed in parallel. The matrix inversion in the $\zeta$ update can be computed offline as it involves no decision variables, so only a dense matrix multiplication is required. We note that multiplication by $\Psi$ is not considered a matrix multiplication in the analysis presented in Table II, as it is the equivalent of a linear filtering operation and therefore scales linearly with $N$. This implies that the residual updates also scale linearly with $N$ if they are analytically block multiplied. The complexity of the ADMM iteration is therefore $\mathcal{O}(N^{n})$ where $n\leq 2$ is determined by the method used for matrix multiplication.

V Numerical Experiments

To compare the performance of the algorithms without reference to a particular PHEV powertrain, single-shot instances of problem (14) were created with randomly generated parameters. For each instance of the energy manangement problem, a nominal sampling frequency of 1Hz (i.e $\delta=1$ s) was assumed, and it was also assumed that the engine is always on and the clutch is engaged i.e {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{\sigma}(k)}=1 (for the purposes of these experiments the switching heuristic is arbitrary). Using observations from previous experiments [11], predictions were made of driver power demand as {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hat{P}_{drv,k}}\in[-2.5\times 10^{3},10\times 10^{3}]\,W and hardware parameters were generated from the distributions {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{2,k}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{2,k}}\in[0.5\times 10^{-5},1.5\times 10^{-5}]\,W^{-1}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{1,k}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{1,k}}\in[0.5,1.5], and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{0,k}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta_{0,k}}=0\,W, with $V_{oc}=300\,V$ and $R=0.1\,\Omega$. The limits on state and input were set at $\overline{x}=10^{5}\,J$, $\underline{x}=0\,J$, $\overline{u}=15\times 10^{3}\,W$ and $\underline{u}=-15\times 10^{3}\,W$, and an initial state of charge of $x_{0}=0.9\overline{x}$ was assumed. These limits have little physical significance within the context of this experiment, and were chosen to ensure that the problems are feasible and that both the state and input constraints were active. Finally, the effect of varying $\tau$ was not investigated, and was set at $\tau=0.995$ for all projected interior point solutions. The simulations were implemented in Matlab on a 2.6GHz Intel Core i7-6700HQ CPU.

V-A Optimum

It is demonstrated in Section III that the output of Algorithm 2 can be made arbitrarily close to the solution of (17) by using a sufficiently large value of $\overline{\mu}$, and it is therefore necessary to determine a value of $\overline{\mu}$ that can be used to obtain a sufficiently accurate approximation of $u^{\star}$. 10 problems for each horizon length $N=100,200,300,400$ were generated for a total of 40 problems, and the projected interior point method was used to obtain a solution for each with the parameters $\mu_{0}=\overline{\mu}$, and $\overline{\mu}$ iteratively increased from $10^{2}$ to $10^{5}$ in 20 logarithmically spaced points ($k_{\mu}$ is not required as the algorithm will terminate when condition (26) is first met). For each $\bar{\mu}_{i}$, $i=1,\ldots,20$, the control input vector, $u^{\ast}_{i}$, at termination was recorded, forming the sequence $u^{\ast}_{1},\dots,u^{\ast}_{20}$ for each problem. Fig. 3 shows that as $\overline{\mu}$ was increased, the norm of the difference between the values of $u^{(j)}$ when the algorithm terminated with successive values of $\bar{\mu}$, $\|u^{\ast}_{i}-u^{\ast}_{i-1}\|$, decreased within an inverse band for all cases, and that at $\overline{\mu}=10^{5}$ this metric has reduced to less than $1$ for all 40 problem cases. Given that the decision variable, $u$, can take a range of values in the order of $10^{4}$, it was concluded that $\overline{\mu}=10^{5}$ is therefore sufficiently large to provide a highly accurate solution to problem (17), and all future references to $u^{\star}$ refers to control inputs found using Algorithm 2 with $\mu_{0}=\overline{\mu}=10^{5}$.

V-B Algorithm Tuning

Both algorithms have multiple parameters that must be tuned to provide computationally efficient solutions. For the ADMM algorithm, $\rho_{1}$ and $\rho_{2}$ (which can be loosely interpreted as the step length in a gradient descent algorithm) must be determined, whilst $\mu_{0}$ and $k_{\mu}$ must be determined for the projected interior point method. The energy management MPC framework is commonly implemented with a shrinking horizon, and it is therefore important that the same set of parameters provide a similar level of performance for a broad class of problems over both long and short horizons. This section details the results of investigations to determine the most computationally efficient combination of parameters for each algorithm.

To determine the values of $\rho_{1}$, $\rho_{2}$, $\mu_{0}$, and $k_{\mu}$, that provide optimal convergence for the ADMM and projected interior point algorithms, 20 new problems were generated for each horizon length of of $N=100,200,300,400$. These were solved using the projected interior point algorithm with $\overline{\mu}=10^{5}$, $10^{-5}\leq\mu_{0}\leq 10^{1}$, and $1<k_{\mu}\leq 10^{5}/\mu_{0}$ (any value of $k_{\mu}$ greater than this would ensure that $\mu$ is projected onto $\overline{\mu}$ during the first update step (27)). For each problem instance, the number of iterations required for the algorithm to terminate were recorded, and the average for each combination of parameters is shown in Fig. 4. It can be seen that there is a vertically banded region at $\mu_{0}\approx 10^{-1}$ that requires a minimum number of iterations for all horizon lengths, and that there is a profile to the search space that varies little with changes in horizon length. The values $\mu_{0}=10^{-1}$ and $k_{\mu}=10^{4}$ were therefore selected as the optimal parameters.

For the ADMM algorithm a different approach was taken, as the number of iterations required to achieve the same level of accuracy as the projected interior point algorithm with $\overline{\mu}=10^{5}$ made a similar parameter search intractable. Instead, a total of 100 ADMM iterations were completed for each problem (ignoring the stated termination criteria) with $10^{-6}\leq\rho_{1}\leq 10^{-2}$ and $10^{-8}\leq\rho_{2}\leq 10^{-4}$. The average norm of the difference between the control input at the 100th iteration of ADMM, $u^{(100)}$, and the optimum, $u^{\star}$, was recorded for each case. The results are shown in Fig. 4, and there is a clear region within approximately two orders of magnitude of both $\rho_{1}$ and $\rho_{2}$ where the control vector has a minimum error relative to the optimum, and this region does not change significantly with horizon length. The values of $\rho_{1}=6\times 10^{-5}$ and $\rho_{2}=4\times 10^{-7}$ were therefore selected as the optimal parameters.

V-C Computational Performance

After tuning the parameters of both the projected interior point and ADMM algorithm to the class of problems being investigated, it was possible to analyse their comparative computational performance. This was achieved in two steps: firstly the termination criteria for a ‘sufficiently’ accurate solution was determined, then the variation in computational time with horizon length was investigated.

A further 20 test cases were generated consisting of 5 cases for each of $N=100,200,300,400$, and using the values of $\rho_{1}$, $\rho_{2}$, $\mu_{0}$ and $k_{\mu}$ determined during the tuning phase, each problem was solved using ADMM for 100 iterations, and using the projected interior point algorithm with $\overline{\mu}=10^{5}$. The absolute difference between the cost evaluated at iteration $j$ and the optimal cost, $|F(u^{(j)})-F(u^{\star})|$, is shown in Fig. 5. The results clearly demonstrate sublinear convergence for the ADMM algorithm, whilst the projected interior point results show superlinear convergence. Therefore, the projected interior point algorithm can produce an extremely accurate solution within a few tens of iterations, whereas significantly more iterations are required for ADMM.

A threshold of 1% was considered high enough for ‘sufficient’ accuracy, as this is likely to be lower than the level of uncertainty in state measurements used to formulate the problem, and the results shown in Fig. 6 illustrate that the deviation between the control vectors obtained by both the ADMM and projected interior point algorithms and the optimum are almost imperceptible at this level of convergence. A slightly larger deviation is observed between the state trajectories, particularly that obtained with ADMM, however this is because the state trajectory is a function of the integral of the control input, and as the cost is not a function of state-of-charge this does not necessarily indicate greater sub-optimality. A key property of the algorithms is also demonstrated in Fig. 6: the state constraints are only guaranteed for both algorithms when the residuals, $r_{\text{IP}},r_{\text{P}},$ and $r_{\text{D}}$, are precisely zero (this is also why we use the absolute error in Fig. 5, as the cost evaluated for each iteration of can be lower than $F(u^{\star})$). Therefore, the termination criteria do not provide a guarantee of enforcing the state constraints, and we can see that for the final three timesteps the lower state limit is violated by $\approx 2.5\%$ of the feasible state band for the ADMM trajectory. This limitation can be reduced by tightening the algorithms’ convergence thresholds, which makes it more significant for ADMM due to its sublinear rate of convergence.

Based on residuals for the projected interior point and ADMM trajectories shown in Fig. 6, it was assumed that $\epsilon=4\times 10^{3}$ and $\overline{\mu}=1$ enforce a ‘sufficient’ level of convergence. A further 20 problems were generated for horizons $50\leq N\leq 1000$, and the iterations to completion, mean time taken per iteration, and time to completion were recorded for each using both ADMM and the projected interior point algorithm. For comparison, the problems were also solved using CVX with default solver SDPT3 v.4.0 [31] and default error tolerance, for which only the total time was recorded (it is not possible to separate the total time from the individual iterations or the overhead required to parse the problem when using CVX). The results are shown in Fig. 7, where it can be seen that whilst the uncertainty in the number of ADMM iterations is high (from as low as 50 to as high as 400), the band of uncertainty is near constant as the horizon is increased, so it can be assumed that the expected number of iterations is effectively constant with horizon length. The uncertainty for the number of projected interior point iterations is lower, and fewer iterations are required for all horizon lengths, however the number of iterations also increases linearly with horizon length from $\sim 10$ iterations at $N=50$ to $\sim 16$ iterations at $N=1000$.

It was shown in sections III-C and IV-A that, as the horizon length is increased, the computational burden of the projected interior point algorithm is dominated by the $\Delta\theta_{1}$ update and the ADMM algorithm is dominated by the $\zeta$ update, so the methods used to perform these calculations will largely determine the time required per iteration. The Matlab operations x=A\b and x=A∗b were used here, and a quadratic trendline is shown to have an approximate fit for both in Fig. 7, although as $N$ was increased towards $1000$, the projected interior point iterations took over two orders of magnitude longer than the ADMM iterations. It was therefore expected that the total time taken for the ADMM algorithm would scale quadratically with horizon length, whereas the time taken for the projected interior point method would scale cubically with horizon length, and this is supported by the results shown in the bottom plot in Fig. 7. It is also shown that (assuming an interval of 1 second between controller optimizations) the projected interior point is only suitable up to a horizon of $N\approx 500$, whereas even up to the maximum horizon length of $N=1000$, the ADMM algorithm only required $\sim 0.1s$. From the previous scaling properties it can be assumed that ADMM is real time implementable for horizons significantly in excess of $1000$ (the performance of the ADMM algorithm presented here exceeds that presented in [15] due to a fully vectorized software implementation). Therefore, whilst the projected interior point algorithm has been shown to converge to an extremely accurate solution in fewer iterations than the ADMM algorithm, for the hardware used in these experiments, less time is required for a moderate level of accuracy using ADMM, and the ADMM algorithm scales better with horizon length. If the accuracy requirement were tightened, however, it is likely that this performance relationship would change significantly.

In comparison, CVX was unable to obtain solutions in less than 1s for any horizon length, and was at least an order of magnitude slower than both algorithms over all horizon lengths; compared to ADMM it was a factor of 1000 slower for $N=1000$. Although CVX is solving the problem to a different error tolerance, it would be expected that ADMM would still be faster were its termination threshold significantly tightened. This is the first demonstration of a method capable of solving the energy management problem in real time, over long horizons ($\geq$1000 samples) when nonlinear system dynamics are considered and hard limits on both power and state of charge are enforced over the entire horizon.

To conclude the numerical experiments, we note that the Newton method for the projected interior point requires the solution to non-diagonal linear systems of equations, which in turn is typically solved using BLAS [32]. This is not an issue when solving the problems on desktop hardware as demonstrated here, but may not be an option for the embedded hardware used for an online solution in a vehicle. In this case only the ADMM algorithm is suitable, as although it also requires a Newton method for the individual control variable updates, this can be performed element-wise and therefore does not require a matrix inversion step.

VI Conclusion

This paper proposes a projected interior point method for the solution of a convex formulation of the optimization problem associated with nonlinear MPC for energy management in hybrid electric vehicles. The performance w.r.t the tailored ADMM algorithm of [15] is demonstrated through numerical experiments, and the projected interior point algorithm is shown to have faster convergence (superlinear) for the class of problems investigated, although the ADMM algorithm is shown to have superior numerical performance and scaling properties when a modest level of accuracy is required. Both algorithms are also shown to have superior computational performance to general purpose convex optimization software.