An ADMM Algorithm for MPC-based Energy Management in Hybrid Electric Vehicles with Nonlinear Losses

TL;DR

This paper introduces a convex formulation and an ADMM algorithm for energy management in hybrid electric vehicles, achieving significant speed improvements over traditional dynamic programming methods.

Contribution

It presents a novel convex formulation that linearizes nonlinear dynamics and applies ADMM for efficient optimization in hybrid vehicle energy management.

Findings

ADMM outperforms DP in solution time by up to two orders of magnitude

Convex formulation enables efficient optimization with nonlinear dynamics

Simulation results validate the effectiveness of the proposed method

Abstract

In this paper we present a convex formulation of the Model Predictive Control (MPC) optimisation for energy management in hybrid electric vehicles, and an Alternating Direction Method of Multipliers (ADMM) algorithm for its solution. We develop a new proof of convexity for the problem that allows the nonlinear dynamics to be modelled as a linear system, then demonstrate the performance of ADMM in comparison with Dynamic Programming (DP) through simulation. The results demonstrate up to two orders of magnitude improvement in solution time for comparable accuracy against DP.

Figures

Equations

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

$${\omega }_k=k_i{v}_k,\text{if}{\underline{v}}_i\le v_k\le {\overline{v}}_i,i=1,\dots ,N_g,$$

$$\mathcal{P}=\{k|{\hat{P}}_{drv,k}>0,{\hat{\omega }}_k\ge {\underline{\omega }}_{eng}\},$$

$$\mathcal{C}=\{k|{\hat{P}}_{drv,k}>0,{\hat{\omega }}_k<{\underline{\omega }}_{eng}\},\mathcal{B}=\{k|{\hat{P}}_{drv,k}\le 0\},$$

$${\hat{\omega }}_{em,k}={\hat{\omega }}_k\forall k,{\hat{\omega }}_{eng,k}=\{\begin{array}{ll}{\hat{\omega }}_k& k\in \mathcal{P}\\ {\underline{\omega }}_{eng}& k\notin \mathcal{P}\end{array}.$$

$$P_f(P_{eng},{\omega }_{eng})={\alpha }_{2,{\omega }_{eng}}{P}_{eng}^2+{\alpha }_{1,{\omega }_{eng}}{P}_{eng}+{\alpha }_{0,{\omega }_{eng}},$$

$$P_c(P_{em},{\omega }_{em})={\beta }_{2,{\omega }_{em}}{P}_{em}^2+{\beta }_{1,{\omega }_{em}}{P}_{em}+{\beta }_{0,{\omega }_{em}},$$

$${\hat{P}}_{f,k}$$

$${\hat{P}}_{c,k}$$

$${\hat{P}}_{b,k}=g_k({\hat{P}}_{em,k})=\frac{V_{oc}^2}{2R}\left(1-\sqrt{1-\frac{4R}{V_{oc}^2}{h}_k({\hat{P}}_{em,k})}\right)$$

$${\underline{\hat{P}}}_{eng,k}$$

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

$${\hat{P}}_{eng}^{\star }=\arg \underset{{\hat{P}}_{eng}}{\min }$$

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

$$\displaystyle\begin{aligned} &\hat{E}_{k+1}=\hat{E}_{k}-g_{k}(\hat{P}_{em,k})\\ &\underline{E}\leq\hat{E}_{k+1}\leq\overline{E}\end{aligned}\Biggr{\}}$$

$$\displaystyle\begin{aligned} &\hat{P}_{drv,k}=\hat{P}_{eng,k}+\hat{P}_{em,k}\\ &\underline{\hat{P}}_{eng,k}\leq\hat{P}_{eng,k}\leq\overline{\hat{P}}_{eng,k}\end{aligned}\Biggr{\}}$$

$${\hat{P}}_{eng,k}={\underline{\hat{P}}}_{eng,k}$$

$${\hat{P}}_{em,k}={\hat{P}}_{drv,k}$$

$${\hat{P}}_{em,k}=\gamma {\hat{P}}_{drv,k}$$

$${\underline{\hat{P}}}_{eng,k}^{+}$$

$${\underline{\hat{P}}}_{em,k}^{+}$$

$${\hat{P}}_{b,k}=g_k({\hat{P}}_{em,k})\iff {\hat{P}}_{em,k}=g_k^{-1}({\hat{P}}_{b,k}),$$

$$g_k^{-1}({\hat{P}}_{b,k})=-\frac{{\beta }_{1,k}}{2{\beta }_{2,k}}+\sqrt{-\frac{R{\hat{P}}_{b,k}^2}{{\beta }_{2,k}{V}_{oc}^2}+\frac{{\hat{P}}_{b,k}-{\beta }_{0,k}}{{\beta }_{2,k}}+\frac{{\beta }_{1,k}^2}{4{\beta }_{2,k}^2}}.$$

$$\displaystyle\frac{\text{d}y}{\text{d}x}=g^{\prime}_{k}(x),\$$

$$\frac{{\text{d}}^{2}x}{\text{d}{y}^2}$$

$$=\frac{-1}{(g_k^{\prime }(x){)}^2}{g}_k^{\prime \prime }(x)\frac{\text{d}x}{\text{d}y}=\frac{-g_k^{\prime \prime }(g_k^{-1}(y)))}{{\left(g_k^{\prime }(g_k^{-1}(y))\right)}^3}\le 0$$

$$f_k({\hat{P}}_{eng,k})=f_k({\hat{P}}_{drv,k}-g_k^{-1}({\hat{P}}_{b,k})).$$

$$P_{drv,k}-g_k^{-1}\left(\lambda x_1+(1-\lambda )x_2\right)$$

$$\le P_{drv,k}-\lambda g_k^{-1}(x_1)-(1-\lambda )g_k^{-1}(x_2).$$

$$f_k\left(P_{drv,k}-g_k^{-1}\left(\lambda x_1+(1-\lambda )x_2\right)\right)$$

$$\le f_k\left(P_{drv,k}-\lambda g_k^{-1}(x_1)-(1-\lambda )g_k^{-1}(x_2)\right)$$

$$\le \lambda f_k\left(P_{drv,k}-g_k^{-1}(x_1)\right)+(1-\lambda )f_k\left(P_{drv,k}-g_k^{-1}(x_2)\right)$$

$${\hat{P}}_b^{\star }=\arg \underset{{\hat{P}}_b}{\min }$$

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

$$\begin{array}{l}{\hat{E}}_{k+1}={\hat{E}}_k-{\hat{P}}_{b,k}\\ \underline{E}\le {\hat{E}}_{k+1}\le \overline{E}\\ {\underline{\hat{P}}}_{b,k}\le {\hat{P}}_{b,k}\le {\overline{\hat{P}}}_{b,k}\end{array}\}\forall k.$$

$$\begin{array}{l}{\underline{\hat{P}}}_{b,k}=g_k({\underline{\hat{P}}}_{em,k}^{+})\\ {\overline{\hat{P}}}_{b,k}=g_k({\hat{P}}_{drv,k}-{\underline{\hat{P}}}_{eng,k}^{+})\end{array}\}k\in \mathcal{P}$$

$${\underline{\hat{P}}}_{b,k}={\overline{\hat{P}}}_{b,k}=\{\begin{array}{ll}{g}_k({\hat{P}}_{drv,k})& k\in \mathcal{C}\\ g_k(\gamma {\hat{P}}_{drv,k})& k\in \mathcal{B}\end{array}.$$

$${\hat{P}}_b^{\star }=\arg \underset{{\hat{P}}_b}{\min }$$

$$\zeta =-{\hat{P}}_b$$

Full text

I INTRODUCTION

Automotive powertrains have become increasingly electrified in recent decades, in an effort to reduce both hydrocarbon consumption and tailpipe emissions [1]. Hybrid electric vehicles are a common configuration as they can reduce fuel consumption through regenerative braking, engine shut down, and reduced engine size, but the second power source introduces an additional degree of freedom that makes the powertrain underconstrained, so the fraction of power delivered by the engine and motor must be controlled. The total energy consumed during a journey can be significantly reduced by varying this fraction [2], but a control system that is optimal (in the sense of minimising fuel consumption) whilst also being implementable is still an open problem [3].

Model Predictive Control (MPC) is a control framework where for each control input update, the solution of an open-loop finite-horizon optimal control problem is solved online, and the first element of the optimal predicted control vector is implemented as the current control input. It has been applied to the hybrid vehicle energy management problem as the optimisation can explicitly account for constraints on system variables, and the recursive feedback provides a degree of robustness to modelling, measurement, and disturbance uncertainty [4, 5].

An important issue with MPC is the online solution of the optimal control problem: it must capture the important dynamics of the controlled system, but the optimisation problem may be high-dimensional, constrained, and possibly nonconvex. Dynamic Programming (DP) has been used to provide the globally minimising argument [6], however the approach is too slow to be implemented in real time and is computationally expensive. Conversely, linearisation has been used to reduce the optimisation to a quadratic program [7], which is a well-known problem and can be solved rapidly, however the linearisation process renders the applied control inputs sub-optimal. To alleviate these difficulties, the problem can be reformulated as a convex optimisation [8, 9, 10]. This allows the dynamics of the system to be accurately modelled using convex functions, and enables a guarantee of convergence to the global minimum.

The first contribution of this paper is an alternative proof of convexity for the energy management problem. Previous work has achieved convexity by relaxing the system dynamic constraints to inequalities, and demonstrating that the optimal solution is obtained when the dynamic constraints are satisfied with equality [8]. We instead show that the problem is convex when expressed in terms of battery power output. This modification allows the dynamics to be represented without approximation as a linear system, which simplifies the implementation of convex optimisation algorithms.

The second contribution is an Alternating Direction Method of Multipliers (ADMM) [11] algorithm for the solution of the convex problem. Our work in [4] demonstrated a projected Newton method that only considers the terminal state of charge constraint, and other work in the area has made use of generic convex optimisation solvers; the method presented here enforces state of charge constraints at all times, and is specifically designed to exploit separability of the cost and dynamics. We have previously demonstrated the properties of the algorithm without the assumption of convexity in [12]. Here, we show the comparative performance against DP in a convex formulation through simulation.

The paper is organised as follows: Chapter 2 describes an MPC framework for energy management, and Chapter 3 provides a proof of convexity for the associated optimisation problem. Chapter 4 details the DP and ADMM algorithms, Chapter 5 presents numerical simulations comparing their performance, and Chapter 6 provides conclusions.

II Model Predictive Control Framework

At each time-instant that the control variable is updated, we assume that a prediction is made of future vehicle velocity, $\hat{v}=(\hat{v}_{0},\dots,\hat{v}_{N-1})$, and road gradient, $\hat{\theta}=(\hat{\theta}_{0},\dots,\hat{\theta}_{N-1})$, at a sampling frequency of 1Hz and a horizon length of $N$. The first derivative of $\hat{v}$ is approximated using central difference, and the predicted driver power demand $\hat{P}_{drv}=(\hat{P}_{drv,0},\dots,\hat{P}_{drv,N-1})$ can be calculated using a backwards-causal longitudinal model

[TABLE]

where $m$ is the vehicle mass, $\rho_{a}$ is air density, $c_{d}$ is the drag coefficient, $A$ is the frontal area, $c_{r}$ is the rolling resistance coefficient, and $g$ is the acceleration due to gravity.

Figure 1 shows a simplified diagram of a parallel plug-in hybrid electric vehicle powertrain. Fuel power, $P_{f}$, is fed to an internal combustion engine that delivers mechanical power, $P_{eng}$, through a clutch, which can be disengaged so that the engine can idle (or be turned off) and the vehicle can be driven in an all-electric mode. Chemical power, $P_{b}$, from the battery is delivered as electrical power, $P_{c}$, through a circuit to a motor that delivers mechanical power, $P_{em}$, to the drivetrain. This flow of power can be reversed to charge the battery. It is assumed that all drivetrain components are 100% mechanically efficient, so the motor and engine output powers are combined additively through a mechanical coupling device to provide the driver power demand, $P_{drv}$. We ignore the effect of high speed dynamics such as clutch engagement and engine lag.

A controller is required to regulate the behaviour of this system as it is under-constrained in three degrees of freedom: the fraction of total driver demand power delivered from the engine, the clutch engagement, and the gear selection. Both the clutch disengagement and gear selection (assuming a discrete variable transmission) introduce integer decision variables that make the resulting control problem non-convex (and significantly more complex), so we assume here that these variables are determined by heuristics.

We assume a simple shifting strategy by which the powertrain speed $\omega_{k}$ can be found from

[TABLE]

where $k_{i}$ are constants and $N_{g}$ is the number of gear ratios. We also assume that when the vehicle is braking, the clutch is open, the engine is idle, and the brakes provide a fixed fraction of the braking force $\gamma$. Finally, we assume that when the demand power is positive and the powertrain speed is lower than the minimum engine speed $\underline{\omega}_{eng}$, the clutch is open and the electric motor provides all of the driving power. Therefore, the engine power fraction only needs to be calculated when the driver power demand is positive, and the powertrain speed is greater than the minimum engine speed. These statements define the sets:

[TABLE]

from which predicted motor and engine speeds are given by

[TABLE]

The engine loss map is modelled as a set of convex quadratic functions of $P_{eng}$ of the form

[TABLE]

for $\omega_{eng}$ in the interval $\underline{\omega}_{eng},\dots,\overline{\omega}_{eng}$, and the motor loss map is modelled using convex quadratic functions as

[TABLE]

for $\omega_{em}$ in the interval $\underline{\omega}_{eng},\dots,\overline{\omega}_{eng}$. Loss functions are then found of the form

[TABLE]

by linear interpolation between the loss map coefficients on the basis of $\hat{\omega}_{eng,k}$ and $\hat{\omega}_{em,k}$. The battery dynamics are

[TABLE]

where we assume constant $V_{oc}$ and $R$ for convexity. Limits on engine power are given for $k\in\mathcal{P}$ as

[TABLE]

where $\underline{T}$ and $\overline{T}$ are lower and upper limits on engine and motor torque, $r^{+}_{k}$ is the largest real root of $1-\frac{4R}{V_{oc}^{2}}h_{k}(\hat{P}_{em,k})$, and upper and lower limits on battery state of charge are given by $\overline{E}$ and $\underline{E}$. The MPC optimisation is then given by

[TABLE]

and the first element of $\hat{P}^{\star}_{eng}$ is then applied as the control input to the engine.

III CONVEXITY

We now demonstrate a convex reformulation of (2). Firstly, we restrict the domain of each $f_{k}(\hat{P}_{eng,k})$ and $g_{k}(\hat{P}_{em,k})$ so that they are non-decreasing, i.e impose new lower limits:

[TABLE]

This is expected as an increase in motor or engine output power would require an increase in input power, and it is an assumption also made in [8]. This domain restriction ensures that $g_{k}(\cdotp)$ is a one-to-one function, so

[TABLE]

where $g_{k}^{-1}(\hat{P}_{b,k})$ is given by

[TABLE]

Lemma III.1

Under the assumptions on $g_{k}(\cdotp)$ (convex, twice differentiable, non-decreasing, one-to-one), the inverse map $g_{k}^{-1}(\cdotp)$ is concave, twice differentiable, and increasing.

Proof:

The assumptions on $g_{k}(\cdotp)$ imply that $g^{\prime}_{k}(x)\geq 0,\ g^{\prime\prime}_{k}(x)\geq 0$ (where $g^{\prime}$ and $g^{\prime\prime}$ are the first and second derivatives of $g$), and, if $y=g_{k}(x)$, then $x=g_{k}^{-1}(y)$ is unique for all $x$. Therefore

[TABLE]

for all $y$, which demonstrates that $g_{k}^{-1}(\cdotp)$ is differentiable and monotonically increasing, and it is known that

[TABLE]

for all $y$, which demonstrates that $g_{k}^{-1}(\cdotp)$ is twice differentiable and concave. ∎

Using this definition of $g_{k}^{-1}(\cdotp)$, the cost function in problem (2) can be rewritten as

[TABLE]

Lemma III.2

Under the assumptions on $f_{k}(\cdotp)$ (convex, non-decreasing), the function $f_{k}(\hat{P}_{drv,k}-g_{k}^{-1}(\hat{P}_{b,k}))$ is convex and non-increasing with $\hat{P}_{b,k}$.

Proof:

The composition of a non-decreasing function $f_{k}(\cdotp)$ and decreasing function $\hat{P}_{drv,k}-g_{k}^{-1}(\hat{P}_{b,k})$ is necessarily non-increasing. The concave property of $g^{-1}_{k}(\cdotp)$ implies that for any $x_{1},x_{2}\in\mathcal{X}$ (where $\mathcal{X}$ is a domain such that $g^{-1}_{k}(x)$ is concave), and $\lambda\in[0,1]$,

[TABLE]

Therefore, using the properties that $f_{k}(\cdotp)$ is non-decreasing and convex, we have

[TABLE]

for all $\lambda\in[0,1]$, which demonstrates that $f_{k}(\hat{P}_{drv,k}-g^{-1}(\hat{P}_{b,k}))$ is convex in $\hat{P}_{b,k}$. ∎

Problem (2) can now be re-written as

[TABLE]

This is now a convex optimisation problem, where

[TABLE]

Existing analyses of convexity for similar optimal energy management problems replace the system dynamic equations, which appear in (2) as equality constraints, with inequalities, and then show that for a given solution the dynamic equations will be satisfied with equality [8]. The formulation presented here is superior for this application as it renders the dynamics as a linear system, which simplifies the implementation of convex optimisation algorithms.

IV OPTIMISATION METHODS

From Lemma III.2 it is known that $f_{k}(\hat{P}_{drv,k}-g_{k}^{-1}(\hat{P}_{b,k}))$ is non-increasing in $\hat{P}_{b,k}$, so it can be shown that if $\underline{E}_{k}\leq E_{k}-\overline{\hat{P}}_{b,k}\leq\overline{E}_{k}$ $\forall k$, then the minimising argument of (3) is given trivially by $\hat{P}_{b}^{\star}=\overline{\hat{P}}_{b}$, and $\hat{P}_{eng,k}^{\star}=\hat{P}_{drv,k}-g_{k}^{-1}(\hat{P}_{b,k}^{\star})$ for $k\in\mathcal{P}$. If this is not the case, an optimisation procedure must be applied.

Previous research on convex optimisation for the energy management problem has used general convex optimisation solvers to obtain the control inputs; here, we instead implement an ADMM algorithm similarly to [12]. To demonstrate the relative performance characteristics of the algorithm we compare it in simulation with DP.

IV-A Alternating Direction Method of Multipliers

First, we introduce a dummy variable $\zeta$, and re-write (3) as

[TABLE]

where $\Phi$ is an $N-1$ column of ones, $\Psi$ is an $N-1\times N-1$ lower triangular matrix of ones, and an indicator function is defined by

[TABLE]

The augmented Lagrangian associated with problem (4) is

[TABLE]

where $\rho_{1},\rho_{2}>0$ are constants discussed in Section V-A. We define projection functions for $x=[\,x_{1}\ \cdots\ x_{N}]^{\top}$ as

[TABLE]

then the ADMM iteration is given by

[TABLE]

which is initialised with the values

[TABLE]

The $\zeta$ variable is included in the formulation to ensure that the $\hat{P}_{b}$ update is separable, and the individual $\hat{P}_{b,k}$ updates are convex optimisation problems that we solve using an unconstrained Newton method with a backtracking line-search. The matrix inversion associated with the $\zeta$ update can be computed off-line as it involves no decision variables, and the remaining variable updates are trivial. This iteration is repeated until the residuals $\|r^{j+1}\|_{2}\leq\epsilon$ and $\|s^{j+1}\|_{2}\leq\epsilon$, where $\epsilon$ is the convergence threshold, and

[TABLE]

In [12], we demonstrate that this algorithm will converge to a point satisfying the first order conditions of problem (3), which must be the optimal point in this convex formulation. The solution for problem (2) is then given by

[TABLE]

IV-B Dynamic Programming

Although it is well known that DP is computationally expensive and unlikely to perform as well as a dedicated convex optimisation algorithm, we use it here to demonstrate baseline performance as it is commonly used to solve various forms of the energy management problem and is guaranteed to provide the optimal solution for a sufficiently high mesh density. In the absence of other tailored algorithms for this problem, a comparison with a general purpose convex optimization solver would be similarly flawed, as the results would be obfuscated by the performance of the differing software implementations.

Definition IV.1

Define $J\in\mathbb{R}^{N_{E}\times N+1}$ as a matrix of costs, where $J[E,k]$ is the cost-to-go at state index $E=\underline{E},\dots,\overline{E}$ (we assume the values of $E$ are evenly spaced) and timestep index $k=0,\dots,N$. The function $J(E,k)$ returns the cost-to-go, linearly interpolated from the available values of $E$ and $k$, and returns $\infty$ outside the feasible range of $E$. We use the notation $J[:,k]$ to refer to column $k$ of $J$.

Definition IV.2

Define $U\in\mathbb{R}^{N_{E}\times N+1}$ as a matrix of optimal control inputs, where the value $U[E,k]$ is the optimal control input at state index $E$ and timestep index $k$. The function $U(E,k)$ returns the optimal control input linearly interpolated from the available values of $E$ and $k$, and the notation $U[:,k]$ refers to the column $k$ of $U$.

Definition IV.3

Define $P_{d}\in\mathbb{R}^{N_{P}}$ as the vector with elements equal to the discrete, evenly-spaced values that may be assigned to the control input $P_{eng,k}$, such that $P_{d}=(\underline{P}_{eng,k}^{+},\dots,\overline{P}_{eng,k})=(P_{d,1},\dots,P_{d,N_{P}})$ and $P_{d}[i]=P_{d,i}$.

The algorithm is initialised by setting $U[:,N]=J[:,N]=0$, then for $E=\underline{E},\dots,\overline{E}$ the cost to go is calculated at $k=N-1,N-2,\dots,0$ by solving

[TABLE]

The integer $i$ that minimises this function is found by evaluating all possible values, and selecting the minimum. The same method is used to evaluate

[TABLE]

The optimal trajectory and control inputs can be found from

[TABLE]

for $k=0,\dots,N-1$.

V SIMULATIONS

We are interested in the comparative performance of ADMM and DP and not the properties of the MPC framework itself, so we only consider single-shot instances of the MPC optimisation with assumed driver demand profiles.

We use the FTP-75 cycle to define a vector, $v$, of velocity samples, and generate a vector $\theta$, of gradient samples

[TABLE]

for $t=0,\dots,T$ at 1Hz, where $\theta_{s}$ is a constant defined for each optimisation. The predicted velocity and gradient vectors used for each problem instance are given by

[TABLE]

where $(1-\mu)T$ is rounded to the nearest integer, and $\mu\in[0,1]$. The predicted power demand is then calculated from (1). Essentially, $\theta_{s}$ is a variable used to modulate the magnitude of the predicted power, and $\mu$ is a variable used to modulate the length of the predicted power vector.

We use a vehicle model to represent a 1800kg parallel plug-in hybrid passenger vehicle with a 100kW petrol internal combustion engine, a 50kW electric motor, and a 21.5Ah lithium-ion battery. We assume a fixed open circuit voltage of 350V and internal resistance of $0.1\Omega$ for the battery, and the fixed braking fraction, $\gamma$, is set at 0.4. An initial state of charge $\hat{E}_{0}=(50+10\mu)\%$ is used to represent battery discharge as the drive cycle is shortened, and upper and lower limits on state of charge are set as 62% and 50%.

V-A Optimisation Parameters

For DP we assume that $N_{P}=N_{E}/10$, as the feasible battery storage band is approximately an order of magnitude greater than the maximum power available from the motor, and this ensures that the $P_{d}$ and $J$ meshes have similar numerical density. We then use the value of $N_{E}$ to modulate the desired accuracy of the solution. From numerical experiments we find that $N_{E}=1000$ is sufficient to achieve a high accuracy approximation of the global minimum, and all further references to the global minimum/optimum refer to control inputs and state trajectories that are found this way.

For the ADMM algorithm, $\rho_{1}$ and $\rho_{2}$ must be tuned to match the dynamics of the optimisation problem, and we tune these parameters to the instance where $\theta_{s}=2^{\circ}$ and $\mu=1$. First we obtain the global minimum, then run 300 ADMM iterations and record a cost that measures the deviation of the fuel consumption vector at algorithm termination, $\hat{P}_{f}^{\dagger}$, from the global optimum fuel consumption vector, $\hat{P}_{f}^{\star}$, with

[TABLE]

The results shown in Figure 2 illustrate that for this particular instance of the optimisation problem, variations of $\rho_{1}$ within two orders of magnitude and variations of $\rho_{1}$ within one order of magnitude induce little suboptimality in the resulting optimisation. We use the parameters $\rho_{1}=2.34\times 10^{-4}$ and $\rho_{2}=8.86\times 10^{-9}$ for subsequent ADMM optimisations. The value $\rho_{2}\approx 10^{-8}$ may appear low, however it is appropriate given the relative scaling of the fuel consumption for each timestep ($\approx 10$kJ) and the battery capacity ($\approx 10$MJ) in (5).

V-B Optimisation with variation in power demand

A practically useful algorithm should demonstrate consistent performance across a broad range of optimisation scenarios, and to test this we generate drive cycles with varied magnitudes using $\mu=1$ and $-2^{\circ}\leq\theta_{s}\leq 2^{\circ}$, then solve using DP with $30\leq N_{E}\leq 120$ and ADMM with $10^{5}\leq\epsilon\leq 10^{6}$. These ranges of $N_{E}$ and $\epsilon$ are chosen to produce solutions with comparable accuracy that are computationally tractable. We run each instance of both algorithms 10 times and record the average time for completion and the deviation from the optimum cost measured by (7). The results are shown in Fig. 3.

The ADMM algorithm solves the optimisation fastest with $\theta_{s}=2^{\circ}$ (the cycle it is tuned to) and is slowest to converge for $\theta_{s}=0^{\circ}$. Between these scenarios, all results lie within a linear band with a maximum variation in computation time of less than an order of magnitude, displaying a degree of robustness to variations in drive cycle. Comparatively, the results for DP display more consistency as they are contained within a tighter band, but also show worse performance: for the range of parameters tested, DP takes approximately two orders of magnitude more time for a comparative level of accuracy. For example, using the $\theta_{s}=1^{\circ}$ cycle, ADMM terminates with $J=2.5\times 10^{6}$ after $1.5s$ with $\epsilon=10^{6}$, whereas DP takes $115.8s$ to achieve $J=2.0\times 10^{6}$ using $N_{E}=120$. The bands also diverge, suggesting that ADMM performs even better for higher accuracy, although it is not clear to what extent this relationship can be extrapolated.

V-C Optimisation with variation in prediction horizon

To investigate robustness to variations in horizon length, we generate drive cycles with $\theta_{s}=-2^{\circ}$ and $\mu\in[0.25,1]$, then solve the optimisation and record parameters using the same method as above. The results are shown in Figure 4.

It is shown that the time taken for the ADMM algorithm to terminate decreases dramatically as the length of the horizon is reduced. For example, the time taken with $\epsilon=10^{6}$ is $1.8s$ for $\mu=1$ and $0.9s$ for $\mu=0.75$, then $0.1s$ for $\mu=0.25$. As the time taken for DP decreases in a similar manner, the key observation is that the significant speed improvement with ADMM is sustained as the horizon is reduced. This suggests that if the algorithm can deliver a performance benefit for a broad range of drive cycle characteristics at the start of a journey, it will then also deliver the same performance benefit as the journey is completed.

In Figure 5 we show state trajectories for the $\theta_{s}=-2^{\circ}$ and $\mu=1$ cycle that are obtained with DP and ADMM using a comparative level of computational time, as well as the global optimum. These plots clearly illustrate the superior performance of ADMM for this problem, but also show a significant limitation. Whilst the ADMM state trajectory follows the optimum much more closely, it is also shown to violate the SOC constraints by $\approx 0.02\%$. This is because there is no explicit guarantee in the ADMM termination criteria that the constraints are met: this is only ensured when the primal residual is exactly zero. We do not, however, believe that this is a significant problem with ADMM for this application, as for all optimisations shown in Figures 3 and 4, no state of charge violation was greater that $1\%$, and this level of violation is likely to be lower in magnitude than the uncertainty in the current state of charge.

VI CONCLUSION

A convex formulation of the MPC optimisation associated with energy management in hybrid electric vehicles is presented, with a demonstration of convexity that allows the state dynamics to be modelled as a linear system. We propose an ADMM algorithm for the solution of this problem, and demonstrate approximately two orders of magnitude improvement in computation time compared to DP. A degree of robustness to variations in the magnitude and length of the predicted drive cycle is also demonstrated.