Parallel ADMM for robust quadratic optimal resource allocation problems

TL;DR

This paper presents a parallel GPU-based ADMM solver for quadratic resource allocation problems, demonstrating significant speedups in energy management applications for hybrid electric vehicles.

Contribution

It introduces a novel parallel GPU implementation of ADMM tailored for quadratic problems with linear constraints, including a custom quartic minimizer.

Findings

GPU implementation achieves faster computation times than serial versions.

Application to hybrid vehicle energy management shows practical effectiveness.

Numerical simulations validate the efficiency of the parallel approach.

Abstract

An alternating direction method of multipliers (ADMM) solver is described for optimal resource allocation problems with separable convex quadratic costs and constraints and linear coupling constraints. We describe a parallel implementation of the solver on a graphics processing unit (GPU) using a bespoke quartic function minimizer. An application to robust optimal energy management in hybrid electric vehicles is described, and the results of numerical simulations comparing the computation times of the parallel GPU implementation with those of an equivalent serial implementation are presented.

Equations

$$\sum _{i=1}^m{x}_k^{(i)}\ge y_k,k\in {\mathbb{N}}_n$$

$$\sum _{k=1}^n{g}_k^{(i)}(x_k^{(i)})\le c^{(i)},i\in {\mathbb{N}}_m.$$

$$\underset{x_k^{(i)}\in [{\underline{x}}_k^{(i)},{\bar{x}}_k^{(i)}]}{\min }$$

$$\begin{array}{rlrl}& \sum _{i=1}^m{x}_k^{(i)}\ge y_k& & k\in {\mathbb{N}}_n\\ & \sum _{k=1}^n{g}_k^{(i)}(x_k^{(i)})\le c^{(i)}& & i\in {\mathbb{N}}_m.\end{array}$$

$$\underset{x_1^{(i)}\in \mathbb{R},x_k^{(i,j)}\in [{\underline{x}}_k^{(i)},{\bar{x}}_k^{(i)}]}{\min }\frac{1}{q}\sum _{j=1}^q\sum _{i=1}^m\sum _{k=1}^n{f}_k^{(i,j)}(x_k^{(i,j)})$$

$$\text{subject to}\begin{array}{rlr}t]2& \sum _{i=1}^m{x}_k^{(i,j)}\ge y_k^{(j)}& k\in {\mathbb{N}}_n,j\in {\mathbb{N}}_q\\ & \sum _{k=1}^n{g}_k^{(i,j)}(x_k^{(i,j)})\le c^{(i)}& i\in {\mathbb{N}}_m,j\in {\mathbb{N}}_q\\ & x_1^{(i,j)}=x_1^{(i)}& i\in {\mathbb{N}}_m,j\in {\mathbb{N}}_q\end{array}$$

$$f_k^{(i,j)}(x)$$

$$g_k^{(i,j)}(x)$$

$$f(x)=A{x}^4+B{x}^3+C{x}^2+Dx+E$$

$$x^3+b{x}^2+cx+d=0.$$

$$Q=c/3-b^2/9,R=bc/6-b^3/27-d/2$$

$$\Delta =Q^3+R^2.$$

$$\theta=\cos^{-1}\Bigl{(}R/\sqrt{-Q^{3}}\Bigr{)},$$

$$x_a$$

$$x_b$$

$$x_c$$

$$S=\sqrt[3]{R+\sqrt{\Delta }},T=\sqrt[3]{R-\sqrt{\Delta }},$$

$$x_a=S+T-b/3.$$

$$x_a=x_b=x_c=-b/3.$$

$$y_k^{(j)}\le x_k^{(1,j)}+x_k^{(2,j)},j\in {\mathbb{N}}_q,k\in {\mathbb{N}}_n$$

$$\underset{x_1^{(i)}\in \mathbb{R},x_k^{(i,j)}\in [{\underline{x}}_k^{(i)},{\bar{x}}_k^{(i)}]}{\min }\frac{1}{q}\sum _{j=1}^q\sum _{k=1}^n{f}_k^{(1,j)}(x_k^{(1,j)})$$

$$\begin{array}{rlrl}\text{subject to}& x_k^{(1,j)}+x_k^{(2,j)}\ge y_k^{(j)},& & j\in {\mathbb{N}}_q,k\in {\mathbb{N}}_n\\ & \sum _{k=1}^n{g}_k^{(2,j)}(x_k^{(2,j)})\le \Delta E,& & j\in {\mathbb{N}}_q\\ & x_1^{(i,j)}=x_1^{(i)},& & i=1,2,j\in {\mathbb{N}}_q\end{array}$$

$${\rho }_i\leftarrow \tau {\rho }_i\text{if}r/\sigma >1.2\bar{r}/\bar{\sigma }\text{and}{\rho }_i\leftarrow {\rho }_i/\tau \text{if}r/\sigma <0.8\bar{r}/\bar{\sigma }$$

$$\underset{\begin{array}{c}{x}_1^{(i)}\in \mathbb{R},x_k^{(i,j)}\in [{\underline{x}}_k^{(i)},{\bar{x}}_k^{(i)}],z_k^{(i,j)}\in \mathbb{R}\\ s_k^{(j)}\ge 0,h^{(i,j)}\le c^{(i)}\end{array}}{\min }\frac{1}{q}\sum _{i,j,k}{f}_k^{(i,j)}(x_k^{(i,j)})$$

$$\text{subject to}\begin{array}{rlr}t]2& \sum _i{x}_k^{(i,j)}-y_k^{(j)}=s_k^{(j)}& j\in {\mathbb{N}}_q,k\in {\mathbb{N}}_n\\ & g_k^{(i,j)}(x_k^{(i,j)})=z_k^{(i,j)}& i\in {\mathbb{N}}_n,j\in {\mathbb{N}}_q,k\in {\mathbb{N}}_n\\ & \sum _k{z}_k^{(i,j)}=h^{(i,j)}& i\in {\mathbb{N}}_n,j\in {\mathbb{N}}_q\\ & x_1^{(i,j)}=x_1^{(i)}& i\in {\mathbb{N}}_n,j\in {\mathbb{N}}_q\end{array}$$

$$L(x_1,x,z,s,h,\mu ,\lambda ,\nu ,p):=\frac{1}{q}\sum _{i,j,k}{f}_k^{(i,j)}(x_k^{(i,j)})$$

$$\displaystyle+\sum_{i,j,k}\Bigl{[}\mathcal{I}_{[\underline{x}_{k}^{(i)},\hskip 0.90417pt\bar{x}_{k}^{(i)}]}(x_{k}^{(i,j)})\!+\!\tfrac{\rho_{1}}{2}\bigl{(}z_{k}^{(i,j)}\!-\!g_{k}^{(i,j)}(x_{k}^{(i,j)})\!+\!\lambda_{k}^{(i,j)}\bigr{)}^{2}\Bigr{]}$$

$$\displaystyle+\sum_{i,j}\Bigl{[}\mathcal{I}_{(-\infty,\hskip 0.90417ptc^{(i)}]}(h^{(i,j)})+\tfrac{\rho_{2}}{2}\bigl{(}h^{(i,j)}-{\bf 1}^{\top}z^{(i,j)}+p^{(i,j)}\bigr{)}^{2}\Bigr{]}$$

$$\displaystyle+\sum_{j}\Bigl{[}{\bf 1}^{\top}\mathcal{I}_{[0,\infty)}(s^{(j)})+\tfrac{\rho_{3}}{2}\bigl{\|}s^{(j)}\!-\!\sum_{i}x^{(i,j)}\!+\!y^{(j)}\!+\!\mu^{(j)}\bigr{\|}^{2}\Bigr{]}$$

$$\displaystyle+{\sum_{i,j}}\tfrac{\rho_{4}}{2}\bigl{(}x_{1}^{(i)}-x_{1}^{(i,j)}+\nu^{(i,j)}\bigr{)}^{2}$$

$$\displaystyle x_{k}^{(i,j)}\leftarrow\Pi_{[\underline{x}_{k}^{(i)},\,\bar{x}_{k}^{(i)}]}\Bigl{(}\begin{aligned} &\operatorname*{argmin}_{x^{(i,j)}_{k}}\Bigl{\{}\tfrac{1}{q}f_{k}^{(i,j)}(x_{k}^{(i,j)})\\ &+\tfrac{\rho_{1}}{2}\bigl{(}z_{k}^{(i,j)}-g_{k}^{(i,j)}(x_{k}^{(i,j)})+\lambda_{k}^{(i,j)}\bigr{)}^{2}\\ &+\tfrac{\rho_{3}}{2}\bigl{(}s_{k}^{(j)}-\sum_{i}x_{k}^{(i,j)}+y_{k}^{(j)}+\mu_{k}^{(j)}\bigr{)}^{2}\\ &+\delta_{k,1}\tfrac{\rho_{4}}{2}\bigl{(}x_{1}^{(i)}-x_{1}^{(i,j)}+\nu^{(i,j)}\bigr{)}^{2}\Bigr{\}}\Bigr{)}\end{aligned}$$

$$\displaystyle\begin{aligned} &z^{(i,j)}\leftarrow g^{(i,j)}(x^{(i,j)})-\lambda^{(i,j)}\\ &\hskip 25.00003pt+\!\frac{\rho_{2}}{\rho_{1}\!+\!n\rho_{2}}{\bf 1}\Bigl{(}h^{(i,j)}\!+\!p^{(i,j)}\!-\!{\bf 1}^{\!\top}\!\bigl{(}g^{(i,j)}(x^{(i,j)})\!-\!\lambda^{(i,j)}\bigr{)}\Bigr{)}\end{aligned}$$

$$\displaystyle x_{1}^{(i)}\leftarrow\sum_{j}\bigl{(}x_{1}^{(i,j)}-\nu^{(i,j)}\bigr{)}$$

$$\displaystyle h^{(i,j)}\leftarrow\Pi_{(-\infty,c^{(i)}]}\bigl{(}{\bf 1}^{\top}z^{(i,j)}-p^{(i,j)}\bigr{)}$$

$$\displaystyle s^{(j)}\leftarrow\Pi_{[0,\infty)}\Bigl{(}\,\sum_{i=1}^{m}x^{(i,j)}-y^{(j)}-\mu^{(j)}\Bigr{)}$$

$${\mu }^{(j)}\leftarrow {\mu }^{(j)}+s^{(j)}-\sum _i{x}^{(i,j)}+y^{(j)}$$

$${\lambda }^{(i,j)}\leftarrow {\lambda }^{(i,j)}+z^{(i,j)}-g^{(i,j)}(x^{(i,j)})$$

$${\nu }^{(i,j)}\leftarrow {\nu }^{(i,j)}+x_1^{(i)}-x_1^{(i,j)}$$

$$p^{(i,j)}\leftarrow p^{(i,j)}+h^{(i,j)}-1^{\top }{z}^{(i,j)}$$

$$x_k^{(i,j)}$$

Full text

I Introduction

An important class of control problems can be formulated as optimal resource allocation problems with separable convex cost functions, separable convex constraints and linear coupling constraints. Examples include energy management in hybrid electric vehicles [1, 2], optimal power flow and economic dispatch problems [3], portfolio selection [4], and diverse problems in manufacturing, production planning and finance [5]. The computation required to solve these problems grows with the length of the horizon over which performance is optimized, and the rate of growth becomes more pronounced for scenario programming formulations [6] that provide robustness to demand uncertainty. This limits the applicability of strategies such as model predictive control that require online solutions of resource allocation problems.

The alternating direction method of multipliers (ADMM) is a first order optimization method that can exploit separable cost and constraint functions, thus forming the basis of a solver suitable for parallel implementation [7, 8]. This paper describes a graphics processing unit (GPU) parallel implementation of ADMM that avoids the growth in computation time with horizon that is observed in serial implementations. This is possible because GPUs have large numbers (typically thousands) of processing units and efficient memory architectures. Parallel ADMM implemented on general purpose GPUs therefore provides a low-cost approach for time-critical optimal resource allocation problems.

Under the assumption that the separable cost and constraint functions are defined by quadratic polynomials, we derive in this paper an ADMM iteration in which the most computationally demanding step is the minimization of a set of quartic polynomials. Iterative methods have previously been proposed for solving polynomial functions on the GPU [9, 10]. However, to retain the convergence guarantees of ADMM, we devise an exact quartic minimization algorithm for GPUs. The solver was implemented on the GPU using CUDA [11] and we compare its computation times with an equivalent serial implementation for a range of problem sizes. We also describe the application of ADMM to robust optimal energy management in a plugin hybrid electric vehicle (PHEV). The GPU implementation of this solver was tested in simulations against corresponding serial implementations. The results demonstrate a significant speed-up and indicate that the approach is feasible for robust supervisory control of PHEVs with realistic horizon lengths.

II Problem Statement

Consider a resource allocation problem in which a collection of resources, which could for example represent alternative power sources onboard a hybrid vehicle, is to be used to meet predicted demand at each discrete time instant over a future horizon. For a horizon of $n$ time steps, let $y_{k}$ denote the demand and $x^{(i)}_{k}$ denote the resource allocated from source $i$ at time step $k$, for $k\in\mathbb{N}_{n}:=\{1,\ldots,n\}$ and $i\in\mathbb{N}_{m}:=\{1,\ldots,m\}$. In order to meet the predicted demand, we require

[TABLE]

where $x_{k}^{(i)}$ is subject to bounds, $x_{k}^{(i)}\in[\underline{x}_{k}^{(i)},\,\bar{x}_{k}^{(i)}]$, $k\in\mathbb{N}_{n}$ as well as constraints on total capacity:

[TABLE]

Capacity constraints represent, for example, the total energy, $\smash{c^{(i)}}$, available from the $i$th power source over the $n$-step horizon, with losses accounted for through the functions ${g_{k}^{(i)}}$.

The cost associated with $x_{k}^{(i)}$ is denoted $f_{k}^{(i)}(x_{k}^{(i)})$. This represents, for example, the operating cost or fuel cost (including losses) associated with the $i$th resource at time $k$. The objective is to minimize total cost subject to constraints:

[TABLE]

This problem is convex if $f_{k}^{(i)}$ and $g_{k}^{(i)}$ are convex functions of their arguments. This is a common assumption in energy management problems with nonlinear losses [12, 13].

Assumption 1

The functions $f_{k}^{(i)}$ and $g_{k}^{(i)}$ are convex, for all $i\in\mathbb{N}_{m}$ and $k\in\mathbb{N}_{n}$.

If the demand sequence ${y=\{y_{k}\}_{k\in\mathbb{N}_{n}}}$ is stochastic, we assume that independent samples $y^{(j)}$, $j\in\mathbb{N}_{q}:=\{1,\ldots,q\}$, can be drawn from the probability distribution of $y$.

Assumption 2

The samples $y^{(j)}=\{y_{k}^{(j)}\}_{k\in\mathbb{N}_{n},j\in\mathbb{N}_{q}}$ are such that $y^{(j)}$ and $y^{(l)}$ are independent for $l\neq j$.

Various robust formulations of the nominal allocation problem (1) are possible using a set of sampled sequences $y^{(j)}$, $j\in\mathbb{N}_{q}$. For example, an open loop formulation optimizes a single sequence $\{x_{k}^{(i)}\}_{k\in\mathbb{N}_{n}}$ for each $i\in\mathbb{N}_{m}$ with the predicted power demand sequence in (1) replaced by $\{\max_{j}y_{k}^{(j)}\}_{k\in\mathbb{N}_{n}}$. Alternatively, a receding (or shinking) horizon approach introduces feedback by solving the resource allocation problem at each sampling instant using current capacity constraints and demand estimates. In this case a sequence $\{x_{k}^{(i,j)}\}_{k\in\mathbb{N}_{n}}$ is computed for each $i\in\mathbb{N}_{m}$ and $j\in\mathbb{N}_{q}$ subject to the additional constraint that, for each $i$, the first element in the sequence should be equal for all $j$:

[TABLE]

Although the objective in (2) is an empirical mean, the approach of this paper allows for alternative definitions such as the maximum over $j\in\mathbb{N}_{q}$ of $\sum_{i,k}f_{k}^{(i,j)}(x_{k}^{(i,j)})$. The probability that the solution of (2) violates the constraints of (1) has a known probability distribution that depends in on the number, $q$ of samples and on the structural properties of (2) – see e.g. [6, 14] for details – moreover these results apply for any demand probability distribution.

An ADMM solver for (2) is summarised in Appendix A. The approach exploits the separable nature of the cost and constraints of (2), which are sums of functions that depend only on variables at individual time-steps, to derive an iteration in which the greatest contribution to the overall computational load is the minimization of scalar polynomials in (7a). Consequently the iteration does not require matrix operations (inversions or factorizations) and, as discussed below, it is suitable for implementation on extremely simple computational hardware. For each $i\in\mathbb{N}_{m}$, the polynomial minimization in step (7a) can be performed in parallel for ${k\in\mathbb{N}_{n}}$ and $j\in\mathbb{N}_{q}$, since the update for $\smash{x_{k}^{(i,j)}}$ does not depend on $\smash{x_{l}^{(i,r)}}$ for $(l,r)\neq(k,j)$. A parallel implementation therefore allows the ADMM iteration to be computed quickly even for long horizons and large numbers of scenarios.

In the following sections, we describe a GPU-based parallel implementation of the ADMM iteration defined in (7a-i). We make the assumption that $\smash{f_{k}^{(i,j)}}$ and $\smash{g_{k}^{(i,j)}}$ are convex quadratic functions of their arguments. Under this assumption (2) becomes a quadratically constrained quadratic program, and the iteration (7a) then requires the minimization of a 4th order polynomial, for each $i,j$ and $k$.

Assumption 3

$f_{k}^{(i,j)}$* and $g_{k}^{(i,j)}$ are quadratic functions:*

[TABLE]

for $i\in\mathbb{N}_{m},\,j\in\mathbb{N}_{q},\,k\in\mathbb{N}_{n}$.

III Parallelization

This section discusses parallelization of the ADMM iteration in Appendix A, focussing on the minimization step (7a). Under Assumption 3, this requires the minimization of 4th order polynomials, for which we describe a tailored algorithm in Section III-B. We consider NVIDIA GPUs programmed in CUDA C/C++ [11] based on a single instruction multiple thread (SIMT) execution model. The CPU (host) and GPU (device) have separate memory, making it necessary to copy data between the host and device memory. Copying data back and forth adversely affects computation time, and the computational benefits of parallelization may therefore only become apparent for large problems. We therefore investigate in Section III-C the dependence of the computational speed-up that can be achieved via parallelization on the number of quartic polynomials to be minimized.

III-A CUDA programming model

A CUDA program consists of code running sequentially on the host CPU and code that runs in parallel on the GPU. A function (known as a kernel) involving $N$ parallel executions of an algorithm is implemented by executing $N$ separate CUDA threads. Each CUDA thread resides within a one-, two- or three-dimensional block of threads. A kernel can be executed by multiple equally-shaped thread blocks, which are themselves organized into a grid of thread blocks. Each thread block is required to execute independently, concurrently or sequentially, and this creates a scalable programming model in which blocks can be scheduled on any available multiprocessor within a GPU. During execution, CUDA threads may access data from multiple memory spaces: each thread has private local memory; each thread block has shared memory visible to all threads of the block and with the same lifetime as the block; and all threads have access to the same global device memory. We use the runtime API (Application Programming Interface) to perform memory allocation and deallocation, and to transfer data between the host and device memory.

For the CUDA implementations discussed in Sections III-C and IV-B, individual kernels were devised for specific elementwise operations while cuBLAS library functions were used for vector-matrix multiplications. To compare the performance of non-parallel implementations, serial algorithms were written in C++ with compiler optimizations (e.g. vectorization where possible) set for speed (/Ox) unless stated otherwise. Intel Math Kernel Library BLAS routines were used for vector-matrix multiplications. Both implementations were coded using Microsoft Visual Studio.

III-B Quartic minimization

Iterative solution methods, such as Newton’s method [9] and recursive de Casteljau subdivision [10] have previously been proposed for solving systems of polynomials on a GPU. However, when used to minimize the quartic polynomial in (7a), iterative methods can adversely affect ADMM convergence due to residual errors in their solutions. For parallel GPU implementations of the iteration (7) therefore, a better approach is to compute the closed form solutions of the cubic polynomial obtained after differentiating the quartic in (7a) with respect to $\smash{x_{k}^{(i,j)}}$. Algebraic solutions of cubic equations can be derived in various ways, including Cardano’s method and Vieta’s method e.g. [15, Sec. 14.7]). We use a combination of these methods to first determine whether roots are real and/or repeated before they are calculated. For problems in which real solutions are expected, this allows significant computational savings and avoids the memory requirements of storing complex roots.

Consider the quartic polynomial

[TABLE]

and let $b=3B/4A$, $c=C/2A$ and $d=D/4A$. The extremum points of $f$ are the solutions of the cubic equation

[TABLE]

Let $Q$, $R$ and $\Delta$ be defined by

[TABLE]

If $\Delta\leq 0$, all three roots of (4) are real. In this case let

[TABLE]

then the roots are

[TABLE]

If $\Delta>0$, there is only one real root. In this case, defining

[TABLE]

the real root is given by

[TABLE]

For the special case of $Q=R=0$, the three roots are real and repeated, and Vieta’s formula gives the solution

[TABLE]

We additionally require the solver to determine which root of (4) is the minimizing argument of $f$. This is simplified by sorting the roots and discarding the middle root, which is necessarily a maximizer. The minimizer can then be determined by comparing the value of the objective function at the two remaining roots. The quartic minimization algorithm is summarised in Algorithm 1.

III-C Comparison of CPU/GPU performance

The results presented in this section and Section IV-B were obtained using an Intel i5-7300HQ CPU with base clock speed 2.5 GHZ, 12GB of RAM, and an Nvidia GTX 1060 3GB GPU. All operations were performed with single-precision floating point numbers.

To test performance we generate $N$ quartic equations with randomly chosen coefficients. For the CPU implementation we solve the $N$ equations sequentially using Algorithm 1. For the GPU implementation, all available threads (up to a limit of 18432 threads for the Nvidia GTX 1060 3GB) are scheduled to perform Algorithm 1 in parallel. The test procedure for the CPU-based algorithm was as follows.

Generate $N$ random sets of quartic coefficients. 2. 2.

Record start time $T_{1}$. 3. 3.

Minimize the $N$ quartics using Algorithm 1. 4. 4.

Record end time $T_{2}$ and elapsed time $T$ = $T_{2}-T_{1}$.

The test procedure for the GPU parallel cubic solver is similar to the CPU procedure with step 3 changed to:

2. (a)

Input data to the GPU. 2. (b)

Minimize the $N$ quartics using Algorithm 1 implemented on the GPU. 3. (c)

Copy results back to the CPU.

The test procedure was repeated at least 10 times and the average value taken for the run-time. Tests were carried out with and without CPU compiler speed optimizations.

Figure 1 compares the GPU and CPU execution speeds as $N$ varies. For $N>10^{6}$, a $30$-fold increase in speed is obtained using the GPU in comparison with the CPU implementation with 4 cores. The speed-up factor rises to $210$ compared with a single-core CPU implementation with speed optimizations off. Although the CPU implementation is faster for $N<10^{3}$, the GPU becomes significantly faster for large $N$. This is expected as the number of GPU threads must be sufficiently large to counter the effects of overheads such as memory transfer times.

IV PHEV Supervisory Control

This section describes an application of the robust resource allocation problem (2) and the parallelization approach of Section III to optimal energy management for a Plugin Hybrid Electric Vehicle (PHEV). PHEVs employ hybrid powertrains that allow the use of more than one energy source. Here we consider a PHEV with two power sources: an electric motor/generator with battery storage, and an internal combustion engine. The supervisory control problem involves determining the optimal power split between the internal combustion engine and the electric motor in order to minimize a performance cost, defined here as the fuel consumption over a prediction horizon, while ensuring that the current and predicted future power demand is met and that battery stored energy remains within allowable limits.

We consider a parallel hybrid architecture [1, 2] in which the propulsive power ($y$) demanded from the vehicle’s powertrain is provided by the sum of power from the combustion engine ($x^{(1)}$) and power from the electric motor $(x^{(2)})$. Independent samples, $\smash{\{y_{k}^{(j)}\}_{j\in\mathbb{N}_{q},k\in\mathbb{N}_{n}}}$, of the future demand sequence are assumed to be available, and the power balance

[TABLE]

is imposed over an $n$-step future horizon. For each demand scenario $j$ and time step $k$ the engine speed $\smash{\omega_{k}^{(j)}}$ is assumed to be known and equal to the electric motor speed whenever the clutch is engaged.

The fuel energy consumed during the $k$th time interval is obtained from a quasi-static map, $\smash{f_{k}^{(1,j)}(x_{k}^{(1,j)})}$, relating engine power output to fuel consumption for the given engine speed $\smash{\omega_{k}^{(j)}}$. Likewise, the battery energy output during the $k$th period, $\smash{g^{(2,j)}_{k}(x_{k}^{(2,j)})}$, is determined for a given $\omega_{k}^{(j)}$ by a quasi-static model of battery and electric motor losses. The energy stored in the battery at the end of the horizon must be no less than a given terminal quantity, $E_{n}$, and hence the fall in battery stored energy over the prediction horizon must not exceed $\Delta E=E_{0}-E_{n}$, where $E_{0}$ is the initial battery state. No constraint is applied to total fuel usage and no cost associated with electrical energy usage (i.e. $g_{k}^{(1,j)}=0$ and $\smash{f_{k}^{(2,j)}}=0$), so the optimization problem (2) becomes

[TABLE]

The updates of $x_{k}^{(i,j)}$ in (7a) therefore require the minimization of a quadratic for $i=1$ and a quartic polynomial for $i=2$. The bounds $\underline{x}_{k}^{(i)}$ are determined by constraints on power and torque output of powertrain components (see e.g. [2]); for example engine braking ($x^{(1,j)}_{k}<0$) is ignored, so ${\underline{x}_{k}^{(1)}=0}$. Assumption 1 (convexity) is satisfied if $f_{k}^{(1,j)}$ and $g_{k}^{(2,j)}$ are convex functions, which is usually the case for real powertrain loss maps ([2]). In addition, cost and constraint functions are approximated by fitting quadratic functions to $\smash{f^{(1,j)}_{k}}$ and $\smash{g^{(1,j)}_{k}}$ at each time-step (see [2] for details), and Assumption 3 is therefore satisfied.

At each sampling instant a supervisory control algorithm based on the online solution of (5) with a shrinking horizon can be summarised as follows.

Algorithm 2** (Supervisory PHEV control)**

**

Draw samples, $\smash{\{y_{k}^{(j)}\}_{j\in\mathbb{N}_{q},\,k\in\mathbb{N}_{n}}}$, from the distribution of the predicted power demand sequence. 2. 2.

For the current battery energy level, $E_{0}$, set $\Delta E=E_{0}-E_{n}$, and compute for each $j\in\mathbb{N}_{q}$ the engine speed sequence $\{\omega_{k}^{(j)}\}_{k\in\mathbb{N}_{n}}$ corresponding to $\{y^{(j)}_{k}\}_{k\in\mathbb{N}_{n}}$, and hence determine $f_{k}^{(1,j)}$ and $g_{k}^{(2,j)}$ for $k\in\mathbb{N}_{n}$. 3. 3.

*Solve (5) by applying the ADMM iteration (7a-i) until the primal and dual residuals satisfy $r<\bar{r}$ and $\sigma<\bar{\sigma}$. * 4. 4.

Apply the engine and electric motor power $x_{1}^{(1)}$, $x_{1}^{(2)}$.

IV-A Modelling assumptions and optimization parameters

We consider a $1900$ kg vehicle with a hybrid powertrain consisting of a $100$ kW gasoline internal combustion engine, $50$ kW electric motor, and a $21.5$ Ah lithium-ion battery. Power demand sequences were generated by adding random sequences (Gaussian white noise of standard deviation $250$ W and $50$ rpm, low-pass filtered with cutoff frequency $0.02$ Hz) to the power demand and engine speed computed for the FTP-75 drive cycle. The sampling interval was fixed at $1$ s. The times at which the vehicle stopped were also modified to simulate variable traffic conditions. This method of generating power demand is used simply to illustrate the computational aspects of the approach; more accurate methods of generating scenarios in this context (e.g. [16, 17]) have been proposed in the literature on stochastic model predictive control.

The initial battery energy level $E_{0}$ was set at $60\%$ of battery capacity, $E_{\max}$, and the terminal level, $E_{n}$, at $50\%$ of capacity. Whenever the power demand $y^{(j)}_{k}$ is negative, $40\%$ is assumed to be recovered through regenerative braking.

Thresholds for the primal and dual residuals were set at $\bar{r}=10^{-6}\Delta E$ and $\bar{\sigma}=10^{-2}$. Step-length parameters $\rho_{1}$, $\rho_{2}$, $\rho_{3}$, $\rho_{4}$ were adapted at each ADMM iteration according to

[TABLE]

with $\tau=1.1$ and initial values $\rho_{1}=10^{-4}$, $\rho_{2}=2\times 10^{-6}$, $\rho_{3}=\rho_{4}=5\times 10^{-6}$. This adaptive approach, which is based on the discussion in Sec. 3.4.1 of [8], was found to reduce the effect of the value of $q$ on the rate of convergence of the ADMM iteration. A typical solution for the battery state of charge for $q=100$ demand scenarios is shown in Fig. 2.

IV-B Comparison of CPU/GPU performance

The CPU-based ADMM implementation was tested using the following procedure:

Perform steps 1 and 2 of Algorithm 2. 2. 2.

Record start time $T_{1}$. 3. 3.

Execute step 3 of Algorithm 2 on the CPU. 4. 4.

Record end time $T_{2}$ and elapsed time $T$ = $T_{2}-T_{1}$.

For the parallel implementation, step 3 was replaced by:

2. (a)

Input data to the GPU. 2. (b)

Execute step 3 of Algorithm 2 on the GPU. 3. (c)

Copy results back to the CPU.

For numbers of demand scenarios $q$ between 5 and 500, Figure 3 shows that GPU runtimes were between 10 and 20 times faster than tests executed on a 4-core CPU with (/Ox) speed optimizations. The number of ADMM iterations required tends to increase with $q$, although for the chosen parameters a significant reduction is observed around $q=20$. The speed-up achieved by the GPU implementation increases with $q$, even though for $q>10$ the number of threads required to execute the ADMM iteration exceeds the limit that can be executed simultaneously on the GPU. To reduce the effect on overall performance of computing primal and dual residuals, the termination conditions were checked only once every 10 iterations in both implementations.

V Conclusions

This paper describes an ADMM algorithm for robust quadratic resource allocation problems. We discuss a parallel implementation of the algorithm and its application to optimal energy management for PHEVs with uncertain power demand. Our results suggest that it is feasible to solve optimal PHEV energy management problems online using inexpensive computing hardware. We also note that the state of the art in low-cost parallel processing hardware is evolving rapidly and that next generation GPUs can provide significant speed improvements over the results reported here.

Appendix A: ADMM iteration

To express (2) in a form suitable for ADMM, we introduce slack variables $s_{k}^{(j)}$ into the demand constraints and introduce dummy variables $\smash{z_{k}^{(i,j)}}$, $\smash{h^{(i,j)}}$ into the capacity constraints:

[TABLE]

This formulation preserves separability of the nonlinear cost and constraints by imposing linear constraints on ${z_{k}^{(i,j)}}$. Let $\mathcal{I}_{\mathbb{I}}(\cdot)$ denote the indicator function of the interval $\mathbb{I}$ (for scalar $x$, $\mathcal{I}_{\mathbb{I}}(x)=0$ if $x\in\mathbb{I}$ and $\mathcal{I}_{\mathbb{I}}(x)=\infty$ if $x\notin\mathbb{I}$; for vector $x$, $\mathcal{I}_{\mathbb{I}}(x)$ is the vector with $k$th element equal to $\mathcal{I}_{\mathbb{I}}(x_{k})$). We define the augmented Lagrangian function

[TABLE]

Here $\mu^{(j)}$, $\lambda^{(i,j)}$, $\nu^{(i,j)}$, $p^{(i,j)}$ and $\nu^{(i,j)}$ are multipliers for the equality constraints in (6); $\smash{y^{(j)}}$, $\smash{x^{(i,j)}}$, $\smash{z^{(i,j)}}$, $\smash{s^{(j)}}$ and $\smash{g^{(i,j)}(x^{(i,j)})}$ are the vectors with $k$th elements equal to $\smash{y_{k}^{(i)}}$, $\smash{x_{k}^{(i,j)}}$, $\smash{z_{k}^{(i,j)}}$, $\smash{s_{k}^{(j)}}$ and $\smash{g_{k}^{(i,j)}(x_{k}^{(i,j)})}$, and ${\bf 1}=[\,1\ \cdots\ 1\,]^{\top}$. Also $\rho_{1},\ldots\rho_{4}$ are non-negative scalars that control the step length of multiplier updates in the ADMM iteration [8, Sec. 3.1]. These design parameters control the rate of convergence and can be tuned to the problem at hand (see e.g. [18]).

At each iteration, $L$ is minimized with respect to primal variables $x^{(i,j)}$, $x_{1}^{(i)}$, $z^{(i,j)}$, $s^{(j)}$, $h^{(i,j)}$. The primal and dual variable updates are given, for $k\in\mathbb{N}_{n}$, $i\in\mathbb{N}_{m}$, $j\in\mathbb{N}_{q}$, by

[TABLE]

where $\Pi_{\mathbb{I}}(x)$ for scalar $x$ denotes projection onto the interval $\mathbb{I}$, while for vector $x$, $\Pi_{\mathbb{I}}(x)$ is the vector obtained by projecting each element of $x$ onto $\mathbb{I}$, and $\delta_{k,1}=1$ if $k=1$, $\delta_{k,1}=0$ if $k\neq 1$. Note that the update (7a) can be written equivalently as

[TABLE]

with $\theta_{k}^{(i,j)}=z_{k}^{(i,j)}+\lambda_{k}^{(i)}$ and $\phi_{k}^{(i,j)}=s_{k}^{(j)}-\sum_{l\neq i}x_{k}^{(l,j)}+y_{k}^{(j)}+\mu_{k}^{(j)}$. The iteration terminates when the primal and dual residuals $r$ and $\sigma$ defined by

[TABLE]

(where $[(\cdot)]_{0}$ denotes the value of an iterate $(\cdot)$ at the previous iteration), fall below predefined thresholds, $\bar{r}$ and $\bar{\sigma}$.