A Lagrangian Policy for Optimal Energy Storage Control

TL;DR

This paper introduces a millisecond-level look-ahead control algorithm for energy storage that is significantly faster than existing convex optimization methods, maintaining optimality and providing insights for system-wide control design.

Contribution

The paper proposes a novel Lagrangian-based control algorithm that achieves constant space complexity and linear runtime, enabling ultra-fast energy storage management.

Findings

Runs up to 100,000 times faster than convex optimization solvers.

Achieves the same control results as state-of-the-art methods.

Provides theoretical insights for centralized and distributed energy storage control.

Abstract

This paper presents a millisecond-level look-ahead control algorithm for energy storage with constant space complexity and worst-case linear run-time complexity. The algorithm connects the optimal control with the Lagrangian multiplier associated with the state-of-charge constraint. It is compared to solving look-ahead control using a state-of-the-art convex optimization solver. Simulation results show that both methods obtain the same control result, while the proposed algorithm runs up to 100,000 times faster and solves most problems within one millisecond. The theoretical results from developing this algorithm also provide key insights into designing optimal energy storage control schemes at the centralized system level as well as under distributed settings.

Tables

Table 1. TABLE I: Piecewise-linear simulation results.

	CVX+Gurobi			Proposed
Trials	${\theta }_0$	$p_1^{\ast }$	cpu [ms]	${\theta }_0$	$p_1^{\ast }$	cpu [ms]
	$T=10$, $J=100$
1	17.53	0.6	408.8	17.53	0.6	0.1
2	16.82	0.0	278.1	16.82	0.0	0.1
3	14.59	0.3	274.5	14.59	0.3	0.1
4	14.21	0.1	281.7	14.21	0.1	0.1
5	7.27	-0.5	283.7	7.27	-0.5	0.1
	$T=10$, $J=1000$
6	18.67	0.0	1851.5	18.66	0.0	0.1
7	20.28	0.8	1791.5	20.28	0.8	0.2
8	19.50	0.0	1788.2	19.50	0.0	0.1
9	17.44	-0.7	1868.4	17.44	-0.7	0.2
10	19.00	0.0	1868.7	19.00	0.0	0.2
	$T=100$, $J=1000$
11	19.09	0.6	18497.7	19.09	0.6	0.1
12	19.95	0.8	19108.2	19.95	0.8	0.2
13	19.44	0.0	18786.1	19.44	0.0	0.1
14	18.91	-1.0	19263.2	18.91	-1.0	0.2
15	19.55	0.0	19080.4	19.55	0.0	0.2

Table 2. TABLE II: Negative Lagrangian dual results.

	Dual results			Control results
Trials	${\underset{¯}{\theta }}_0$	${\theta }_0^{\ast }$	${\overline{\theta }}_0$	${\underset{¯}{p}}_1$	$p_1^{\ast }$	${\overline{p}}_1$
1	-10.7	-10.6	-10.2	-0.14	0.10	0.11
2	-12.9	-12.5	-11.8	0.49	0.53	0.56
2	-8.5	-8.2	-8.2	-0.17	-0.17	-0.16
4	-4.4	-4.4	-4.0	-1.00	-1.00	-1.00
5	-15.0	-14.6	-14.5	0.89	0.91	0.97

Equations

$$p_t^{\ast }\in \arg \underset{p_t}{\min }\sum _{t=1}^T{O}_t(p_t)+C_T(e_T)$$

$$p_t=p_t^{+}-p_t^{-}$$

$$p_t^{+}=0\text{ or }{p}_t^{-}=0\forall t\in \mathcal{T}$$

$$e_t-e_{t-1}=-p_t^{+}/\eta +p_t^{-}\eta :{\theta }_{t-1}$$

$$\begin{array}{rl}{p}_t^{+}& \ge 0:{\underline{\mu }}_t^{+},\\ p_t^{-}& \ge 0:{\underline{\mu }}_t^{-},\\ e_t& \ge 0:{\underline{\nu }}_t,\end{array}\begin{array}{rl}{p}_t^{+}& \le P:{\overline{\mu }}_t^{+}\\ p_t^{-}& \le P:{\overline{\mu }}_t^{-}\\ e_t& \le E:{\overline{\nu }}_t\end{array}$$

$$o_t(p_t^{\ast })+{\overline{\mu }}_t^{+}-{\underline{\mu }}_t^{+}+{\theta }_{t-1}/\eta$$

$$-o_t(p_t^{\ast })+{\overline{\mu }}_t^{-}-{\underline{\mu }}_t^{-}-{\theta }_{t-1}\eta$$

$${\theta }_{t-1}+{\overline{\nu }}_t-{\underline{\nu }}_t-{\theta }_t$$

$${\theta }_T+c_T(e_T^{\ast })$$

$$\begin{array}{rl}{\underline{\mu }}_t^{+}[p_t^{\ast }{]}^{+}& =0,\\ {\underline{\mu }}_t^{-}[-p_t^{\ast }{]}^{+}& =0,\\ {\underline{\nu }}_t{e}_t^{\ast }& =0,\end{array}\begin{array}{rl}{\overline{\mu }}_t^{+}(P-[p_t^{\ast }{]}^{+})& =0,\\ {\overline{\mu }}_t^{-}(P-[-p_t^{\ast }{]}^{+})& =0,\\ {\overline{\nu }}_t(E-e_t^{\ast })& =0,\end{array}$$

$$\displaystyle p^{+}_{t}(x)=\Big{[}\varphi_{t}(-x/\eta)\Big{]}^{P}_{0},\;p^{-}_{t}(x)=\Big{[}-\varphi_{t}(-x\eta)\Big{]}^{P}_{0}$$

$$p_t^{\pi }(x)=p_t^{+}(x)-p_t^{-}(x)$$

$${\phi }_t(x)$$

$$o_t(p_t^{\ast })+{\theta }_{t-1}/\eta$$

$$p_t^{+}=\max \{0,\min \{P,{\phi }_t(-{\theta }_{t-1}/\eta )\}\}.$$

$${\sigma }_t(x)=$$

$${\phi }_t(-x\eta )<0<{\phi }_t(-x/\eta )$$

$$-x\eta <-x/\eta$$

$$p_t^{+}<=P{v}_t,p_t^{-}<=P(1-v_t),v_t\in \{0,1\}$$

$$p_t^{\pi -}(x)=p_t^{+}(x)1_{[p_t^{-}(x)==0]}-p_t^{-}(x)$$

$$p_t^{\pi +}(x)=p_t^{+}(x)-p_t^{-}(x)1_{[p_t^{+}(x)==0]}$$

$${\underline{p}}_t\le -p_t^{-}({\theta }_0)\le p_t^{+}({\theta }_0)\le {\overline{p}}_t$$

$$o_t(p)=c_j\text{ if }{q}_{j-1}\le p<q_j$$

$$O_t=\frac{{\alpha }_t}{2}({\beta }_t-p{)}^2$$

$${\theta }_{t-1}={\theta }_t$$

$${\theta }_{t-1}<{\theta }_t$$

$${\theta }_{t-1}>{\theta }_t$$

$${\theta }_0={\theta }_T=-c_T(e_T^{\ast })=-c_T({\sigma }_T({\theta }_0))$$

$$x$$

$$p_t^{\pi }(x)$$

$${\sigma }_t(x)$$

$$c_T({\sigma }_t(x))$$

$$c_T({\sigma }_T(x))$$

$$x$$

$$p_t^{\pi }(x)$$

Full text

I Introduction

Energy storage devices such as batteries are key resources in future energy systems due to their flexibility and fast response speed, and their convenient installations as either large-scale bulk units or as distributed resources. Most real-time energy storage operations are optimized using predictive control, with applications such as economic dispatch [1], frequency control [2], voltage control [3], renewable integration [4], energy arbitrage [5], peak shaving [6], electric vehicle charging [7], or a combination of several aforementioned applications [8]. These predictive control strategies solve a multi-period optimization problem over a look-ahead horizon at each control step, obtaining the control and state profile over the entire horizon but only applies the first control result, the problem is then updated with a new horizon and state information for the next control step.

The challenge of using look-ahead control in practice is trading off optimality with computational tractability, as a longer look-ahead horizon incorporates more future information and thus improves solution optimality, but increases the computational challenge significantly. For example, real-time economic dispatches in power systems are typically solved over a single period or with a look-ahead horizon less than one hour [9]. However, power system operations have strong daily patterns due to load and weather variations, such as charging storage from solar power during the day and discharge during the night. Thus, being able to incorporate a look-ahead horizon over one day or even longer is crucial for the future power system, but solving such problems over the scale of a realistic power system is extremely computationally challenging [10], especially binary variables must be introduced in certain application to prevent simultaneous charging and discharging, making the problem non-convex [11]. In addition, future uncertainties in power systems are often modeled with scenarios [12], and modeling uncertainties from multiple sources can easily lead to hundreds of scenarios that makes it almost impossible to solve look-ahead economic dispatch with conventional optimization solvers. While methods such as stochastic dual dynamic programming [13] reduce the solution complexity by introducing inter-temporal and scenario decomposition, the computation is still difficult and requires significant memory usage. On the other hand, the optimal control problem must be solved within a reasonable timescale to fully utilize the fast response speed of energy storage devices. For example, a battery ramps from zero to full discharge power within milliseconds [14] thus, a scheme that takes seconds or even minutes to update the control decision is not appropriate for controlling batteries.

Solving storage control from the dual problem is more effective than dealing with the primal problem directly since the storage has only a single state variable with upper and lower bounds. Cruise et a l. [15] concluded the storage control problem can be solved using a search algorithm based on the binding conditions on the state-of-charge, and Hashimi et al [16] has developed an algorithm for energy storage price arbitrage with quadratic time complexity, based on solving the dual problem. Comparably, the technical contributions of this paper and the main advantages of the proposed algorithm is summarized as follows:

We show that energy storage control with a generalized time-varying objective functions can be solved in worst-case linear time complexity and constant space complexity, with respect to the look-ahead horizon. 2. 2.

We conclude the optimal control condition for energy storage without having to go through the full look-ahead horizon, i.e., the current control is optimal with respect to any future realizations that may not be included in the current look-ahead window. 3. 3.

We derive a sufficient condition for the occurrence of simultaneous charging and discharging, and provide upper and lower bounds for the prime and dual results under such non-convex conditions.

The rest of this paper is organized as follows: Section II formulates the problem; Section III presents main analytical results and the algorithm; Section IV demonstrates numerical results; and Section V concludes the paper.

II Formulation and Preliminaries

II-A Problem Formulation

We consider a time period $t\in\mathcal{T}=\{1,\dotsc,T\}$ where $t=1$ is the current control step and $t=2$ to $T$ is the look-ahead horizon. The optimal control profile $p^{*}_{t}$ is a minimizer to the following multi-period optimization problem

[TABLE]

where

$O_{t}(\cdot)|\mathbb{R}\to\mathbb{R}$ is a scalar time-varying convex objective function. Its derivative is denoted as $o_{t}=\dot{O}_{t}$. 2. 2.

$C_{T}(\cdot)|\mathbb{R}\to\mathbb{R}$ is the terminal cost function of the end state of charge $e_{T}$. $C_{T}$ is also convex and its derivative is denoted as $c_{T}=\dot{C}_{T}$. Note that $C_{T}$ can also be used to model the operation beyond $T$ via dynamic programming [17]. 3. 3.

$p_{t}$ is the control decision variable and it is the energy dispatched from the storage during the time period $t$. 4. 4.

$p^{+}_{t}$ is the positive (discharge) component of $p_{t}$. 5. 5.

$p^{-}_{t}$ is the negative (charge) component of $p_{t}$. 6. 6.

$e_{t}$ is the state of charge (SoC) at the end of time period $t$, subjects to an initial value of $e_{0}$. 7. 7.

$\eta\in(0,1]$ is the storage charge and discharge efficiency. 8. 8.

$P\in\mathbb{R}^{+}$ is the maximum energy that can be charged or discharged into the storage during a single period. 9. 9.

$E\in\mathbb{R}^{+}$ is the maximum energy that can be stored in the storage. 10. 10.

$p^{*}_{t}$ is a set of minimizers to the optimization problem. 11. 11.

$\theta_{t}$ is the Lagrangian multiplier associated with the SoC dynamic, its physical meaning is the marginal value of SoC at the end of time $t$ over the future operation $[t+1,T]$. 12. 12.

$\underline{\mu}^{+}_{t}$, $\overline{\mu}^{+}_{t}$, $\underline{\mu}^{-}_{t}$, $\overline{\mu}^{-}_{t}$, $\underline{\nu}_{t}$, $\overline{\nu}_{t}$ are positive dual variables associated with inequality constraints.

The objective function (1a) minimizes the total operating cost over the period $[1,T]$. Constraint (1b) divides the control $p_{t}$ into a positive component $p^{+}_{t}$ and a negative component $p^{-}_{t}$ in order to model the efficiency difference during charge and discharge in the SoC evolution constraint. (1c) is the non-simultaneous charging and discharging constraint that enforces the storage to only charge or discharge at any given time point. (1d) models the SoC evolution subjects to efficiencies. Power and energy ratings are modeled in (1e).

II-B Karush-Kuhn-Tucker conditions

The results in this paper are primarily based on the use of the Karush-Kuhn-Tucker (KKT) conditions [18], which are listed below for (1) ( the non-simultaneous charging and discharging constraint (1c) is non-convex and is excluded from the KKT condition below, treatment of this constraint will be discussed later):

[TABLE]

and the complimentary slackness conditions associated with the inequality dual variables:

[TABLE]

and $\underline{\mu}^{+}_{t}$, $\overline{\mu}^{+}_{t}$, $\underline{\mu}^{-}_{t}$, $\overline{\mu}^{-}_{t}$, $\underline{\nu}_{t}$, $\overline{\nu}_{t}\geq 0$. Note that we replaced the use of $p^{+}_{t}$ and $p^{-}_{t}$ with $p_{t}$ since $\partial p^{+}_{t}=\partial p_{t}$, $\partial p^{-}_{t}=-\partial p_{t}$, $p^{+}_{t}=[p_{t}]^{+}$ and $p^{-}_{t}=[-p_{t}]^{+}$, where $[x]^{+}=\max\{0,x\}$ is the positive value function.

III Main Results

We start by relaxing constraint (1c) so that the rest of the problem is convex. Thus we establish a closed-form connection between the primal and dual problem in Proposition 1. We then present Theorem 3 on identifying the equality relationship between $\theta_{0}$ and any real number $x\in\mathbb{R}$ using numerical simulation, and develop a binary search algorithm that finds the dual result $\theta_{0}$ and , thus, the primal result $p^{*}_{1}$. Then we discuss how we can bound the result when it is necessary to incorporate the non-convex non-simultaneous charging and discharging constraint (1c) using the proposed algorithm.

III-A Optimal Control Policy

We define the policy $\pi$ that calculates a storage control decision $p^{\mathrm{\pi}}_{t}(x)$ for time $t$ from an input $x\in\mathbb{R}$ as

[TABLE]

where $[x]^{y}_{z}=\max\{\min\{x,y\},z\}$ saturates $x$ between $y$ and $z$ ($z\leq y$), and $\varphi_{t}(x):\mathbb{R}\to\mathbb{R}$ is the inverse function of $o_{t}$ (derivative of $O_{t}$) as

[TABLE]

Note that $\varphi_{t}$ is an alternative definition of the inverse function to $o_{t}$ while compatible with a piecewise linear $O_{t}$.

The following proposition states that we can obtain the optimal control $p^{*}_{t}$ by using the Lagrangian multiplier $\theta_{t-1}$ as input to policy $\pi$:

Proposition 1**.**

Policy $\pi$ is a minimizer to problem (1) when using the Lagrangian $\theta_{t-1}$ as the input, i.e., $p^{*}_{t}=p^{\pi}_{t}(\theta_{t-1})$ .

Proof.

We start by rewriting the KKT condition associated with $p^{+}_{t}$ as

[TABLE]

where we substitute the complementary slackness condition into (2a) that replaces $\overline{\mu}^{+}_{t}$ and $\underline{\mu}^{+}_{t}$. It is now trivial to see that we can calculate $p^{*}_{t}$ as $\varphi_{t}(-\theta_{t-1}/\eta)$ and limiting the result between 0 and $P$, hence

[TABLE]

We repeat the similar process for $p^{-}_{t}$ with (2b), and use $p_{t}=p^{+}_{t}-p^{-}_{t}$ which gives us the result in Proposition 1. ∎

The following corollary supplements that with $\theta_{0}$ we can obtain $p^{*}_{1}$ as well as a series of consecutive optimal control decisions by recording the accumulated sum of the control results $\sigma_{t}$ defined as

[TABLE]

with the initial value $\sigma_{0}(x)=e_{0}$, where $[x]^{+}=\max\{0,x\}$ is the positive value function, and $[x]^{-}=\min\{0,x\}$ is the negative value function. $\sigma_{t}(x)$ emulates the SoC evolution but using the control result $p^{\mathrm{\pi}}_{t}(x)$ which may not be optimal. Another difference is that $\sigma_{t}(x)$ is not limited between $[0,E]$, instead, whether any $\sigma_{t}(x)$ falls above $E$ or below [math] is an indicator on the optimality of $p^{\mathrm{\pi}}_{t}(x)$, as defined by the following corollary:

Corollary 2**.**

$p^{*}_{t}=p^{\mathrm{\pi}}_{t}(\theta_{0})$ if $0\leq\sigma_{\tau}(\theta_{0})\leq E$ $\forall$ $\tau\in[1,t]$.

Corollary 2 means that we can maintain optimal control by using $\theta_{0}$ as the input to (5) for control steps beyond $t=1$ if all previous $\sigma_{t}(\theta_{0})$ are within the SoC constraint. Corollary 2 is based on Proposition 1 and the KKT condition associated with $e_{t}$ in (2c) that the $\theta_{t}$ value will not change if both $\overline{\nu}_{t}$ and $\underline{\nu}_{t}$ are zeros, indicating $0\leq e^{*}_{t}\leq E$ and $\sigma_{t}(x)=e^{*}_{t}$. This corollary is thus proved.

III-B Main Theorem on Finding Lagrangian Dual

Theorem 3**.**

Given $x\in\mathbb{R}$, its equality relationship with respect to $\theta_{0}$ can be determined as

If a) $\sigma_{t}(x)$ reaches upper bound first, i.e., $\exists\tau\in\mathcal{T}\text{ s.t. }\sigma_{\tau}(x)>E\text{ and }0\leq\sigma_{\gamma}(x)\leq E$ $\forall\gamma\in[1,\tau)$; or b) $\sigma_{t}(x)$ reached neither bound and $x>-c_{T}(\sigma_{T}(x))$; then $x\geq\theta_{0}$; 2. 2.

If a) $\sigma_{t}(x)$ reaches lower bound first, i.e., $\exists\tau\in\mathcal{T}\text{ s.t. }\sigma_{\tau}(x)<0\text{ and }0\leq\sigma_{\gamma}(x)\leq E$ $\forall\gamma\in[1,\tau)$; or b) $\sigma_{t}(x)$ reached neither bound and $x<-c_{T}(\sigma_{T}(x))$; then $x\leq\theta_{0}$; 3. 3.

If $\sigma_{t}(x)$ reached neither bound and $x=-c_{T}(\sigma_{T}(x))$, then $x=\theta_{0}$.

Proof of this theorem is deferred to Appendix. The intuition is that the Lagrangian dual is the price of the stored energy, its value does not change despite the change with the SoC evolution, except reaching either the upper or lower SoC bound. The SoC series driven by the optimal dual value should never exceed the SoC bounds as the dual value itself reflects the constrained storage capacity. If the SoC exceed the upper SoC bound, it means SoC value is over estimated as the storage does not have enough capacity to store the excessive energy, hence we picked an $x$ that is higher than the optimal Lagrangian dual value. Vice versa, if the SoC exceed the lower bound, meaning we under estimated the dual value.

III-C Solution Algorithm

We design a binary search algorithm that finds $\theta_{0}$ according to Theorem 3, thus we find $p^{*}_{1}=p^{\mathrm{\pi}}_{1}(\theta_{0})$ (Proposition 1) as well as some consecutive optimal control actions (Corollary 2) without needing to explicitly solve Problem (1). The algorithm requires a preset search accuracy $\epsilon$ and is described as follows:

Initialize a search range $L$ and $R$ with which we are confident that $L\leq\theta_{0}\leq R$; 2. 2.

Set $x$ to $(L+R)/2$. If $R-L<\epsilon$, return $x$ as the optimal Lagrangian dual value and $p^{\mathrm{\pi}}_{t}(x)$ as the optimal storage control up to time step $t$; 3. 3.

Run the following iterative simulation

1. (a)

Set $1\to t$ and $e_{0}\to\sigma_{0}(x)$; 2. (b)

Calculate $p^{\mathrm{\pi}}_{t}(x)$ using Eq. (5); 3. (c)

Calculate $\sigma_{t}(x)$ using Eq. (9); 4. (d)

If $\sigma_{t}(x)\geq E$, set $x\to R$, go to Step 2); 5. (e)

If $\sigma_{t}(x)\leq 0$, set $x\to L$, go to Step 2); 6. (f)

If $t<T$, set $t+1\to t$ and go to Step b); 7. (g)

If $x\geq-c_{T}(\sigma_{T}(x))$, set $x\to R$; 8. (h)

If $x<-c_{T}(\sigma_{T}(x))$, set $x\to L$; 4. 4.

Go to Step 2).

An example of a confident search range is that we can assume stored energy always has a positive value and choose $L=0$ and $R=\max\{o_{t}(p)/\eta|t\in\mathcal{T},p\in[-P,P]\}$.

This algorithm achieves the following complexity results:

Constant space complexity: The algorithm achieves $O(1)$ space complexity with respect to the search range and the look-ahead duration $T$, because the equality relationship between $x$ and $\theta_{0}$ can be identified using only the current simulation result $\sigma_{t}(x)$ so that previous simulation results are not required to be stored. 2. 2.

Worst-case linear run-time complexity: The algorithm achieves a worst-case $O(n)$ complexity with respect to the look-ahead horizon $T$ since the worst-case scenario is to simulate all operations steps from $t=1$ to $T$ during each search, but may terminate before reaching $T$ as stated in step 3-d and 3-e. It also achieves $O(\log n)$ time complexity with respect to the search range for using a binary search algorithm.

III-D Non-simultaneous Charging and Discharging

By far we have concluded the optimal storage control when relaxing the non-simultaneous charging and discharging constraint (1b). A sufficient condition for relaxing this constraint without sacrificing result optimality is illustrated in the following proposition:

Proposition 4**.**

A sufficient condition for simultaneous charging and discharging to happen is the Lagrangian dual being negative, i.e., if $\theta_{t-1}\geq 0$ then $p^{+}_{t}p^{-}_{t}=0$.

Proof.

From (5), it is trivial to see that for simultaneous charging and discharge to happen, both terms in (5) must take non-zero values, hence

[TABLE]

and since $\varphi_{t}$ is a monotonic increasing function and $0<\eta\leq 1$, it follows

[TABLE]

hence $x$ must be less than zero. ∎

Recall that $\theta_{t}$ is the value of the stored energy, hence $\theta_{t}$ being negative indicates the stored energy has a negative value, i.e., we have an intention to store as less energy as possible. Thus, when (1b), the storage can charging and discharging at the same time and use round-trip efficiency loss to consume excessive energy even when the storage is full and has no more storage space. This intuition is useful when deciding whether simultaneous charging and discharging should be considered when formulating the problem, for example, this constraint should be considered in price arbitrage for markets with frequent negative prices.

A common method for enforcing non-simultaneous charging and discharging is to add auxiliary binary variables $\bm{v}=\{v_{t}\}$ such that the storage can only charge or discharge at one time, as

[TABLE]

making (1) a mixed-integer programming problem, the Lagrangian dual can thus be calculated given a set of fixed $\bm{v}$.

We assume set $\mathcal{V}$ as the set of all reasonable charging status results, in which $v_{t}$ may be either 0 or 1 during all periods when simultaneous charging and discharge occur, i.e., when (1c) must be enforced. It is worth noting that the optimal result for $v\bm{v}$ must be in $\mathcal{V}$. We denote $\theta_{0}(\bm{v})$ as the resulting Lagrangian dual associated with the charging status $\bm{v}$, then the following proposition stands:

Proposition 5**.**

For all $\bm{v}\in\mathcal{V}$, $\underline{\theta}_{0}\leq\theta_{0}(\bm{v})\leq\overline{\theta}_{0}$, where $\underline{\theta}_{0}$ is calculated using Algorithm 1 by replacing (5) with

[TABLE]

where \mathbf{1}_{[x]}=\{\text{ 0 if x is true and 1 otherwise}\} is the indicator function.

Proof.

First note that compared to (5), (13a) enforces the storage to charge whenever charging and discharging components are both non-zero. Thus when using (13a) to simulate the battery operation in Algorithm 1, the resulting SoC must always be lower than using any $\bm{v}\in\mathcal{V}$. Thus according to Theorem 3, the resulting dual $\theta_{0}$ must be no smaller than any dual $\theta_{0}(\bm{v})$ using charging status $\bm{v}\in\mathcal{V}$. Vice versa, when using (13b), the battery prefer discharging over charging, resulting in a lower bound for $\theta_{0}$. Note that in (13a) and (13b), $p^{\pi}_{t}(x)$ is still a monotonic decreasing function to $x$ since $p^{+}_{t}(x)\geq 0$ and $p^{-}_{t}(x)\geq 0$ for all $x$, $p^{+}_{t}(x)$ is an decreasing function and $p^{-}_{t}(x)$ is an increasing function. Hence the SoC series is still monotonic increasing with respect to $x$, and the convergence optimality of Algorithm 1 will not be effected. ∎

And the result on the dual binding can be extended to bind the primal control results:

Proposition 6**.**

Let $\underline{p}_{t}=p^{\pi-}_{1}(\underline{\theta}_{0})$ and $\overline{p}_{t}=p^{\pi+}_{t}(\overline{\theta}_{0})$, and $p^{*}_{t}$ is the optimal control, then $\underline{p}_{t}\leq p^{*}_{t}\leq\overline{p}_{t}$ for all $t\in[1,\;\tau]$.

Proof.

First given the dual $\theta_{0}$, the optimal control $p^{*}_{1}$ is either $p^{+}_{t}(\theta_{0})$ or $-p^{-}_{t}(\theta_{0})$ as in (5), since the battery either charges or discharges. Then according to Proposition 5, we know the dual must be within $[\underline{\theta},\overline{\theta}]$. Then we have

[TABLE]

which proved this Proposition. ∎

Hence Proposition 5 and Proposition 6 provide bounds for the primal and dual result of the non-convex battery control problem using the proposed algorithm that has worst-case linear time complexity.

IV Numerical Simulation

We use randomly generated data sets to compare the proposed algorithm to solving Problem 1 with different objectives using Gurobi [19] (model generated using CVX [20]). All simulations are performed in Matlab [21] on a 2.3 GHz machine with 16GB memory.The storage parameter is set as $P=1$ p.u., $E=4$ p.u., $e_{0}=2$ p.u., $\eta=0.92$, and the terminal cost function is set to $C_{T}(e_{T})=(E-e_{T})^{2}/2$. The accuracy of the search algorithm is set to $\epsilon=10^{-3}$.

IV-A Piece-wise linear objectives

In this section the proposed algorithm is compared with Gurobi using piecewise linear objective function inspired by the supply curves in power system economic dispatches [22]. The derivative of the objective function $o_{t}$ is written as

[TABLE]

where $j\in[1,J]$ is the piecewise segment index, $J$ is the number of segments, $c_{j}$ is the marginal cost (derivative) of the system when $p$ is between quantities $q_{j-1}$ and $q_{j}$, and the objective is convex if $c_{i}\leq c_{j}$ and $q_{i}\leq q_{j}$ for all $i<j$, $i,j\in[1,J]$. Some examples of the generated cost curve are plotted in Fig. 2.

Similar to the quadratic results, we test the proposed algorithm and Gurobi using different settings and the results are demonstrated in Table I, where trials 1–5 have 10 time steps and 100 cost segments $T=10$, $J=100$, trials 6–10 have 10 time steps and 1,000 cost segments $T=10$, $J=1000$, trials 11–15 have 100 time steps and 1,000 cost segments $T=100$, $J=1000$. The result shows the proposed algorithm obtains the same results in all trials compared to Gurobi, while being hundreds or even thousands of times faster. In particular, in trials 11-15 Gurobi needs around 18 seconds to complete the computation, while the proposed algorithm finishes below 1ms.

In Fig. 3, we further test the computation speed of both methods in solving look-ahead economic dispatches over the size of realistic power systems with 5000 cost segments per five minute dispatch interval. The result shows that the computation time of Gurobi increases significantly with respect to the look-ahead horizon, and in particular at the 6 hour look-ahead, the problem takes more than 5 minutes to solve which is not feasible since the economic dispatch must be calculated within 5 minutes. In contrast, with our proposed algorithm, the average solution speed is below 0.5 milliseconds for look-ahead horizons up to 8 hours, providing a speed-up up to 100,000 times.

IV-B Negative Lagrangian dual example

We consider the following quadratic objective function

[TABLE]

where $\alpha_{t}$ are randomly generated between $[0,10]$, and $\beta_{t}$ between $[-10,0]$. Recall that negative sign is for charging the battery, hence this is a generation tracking problem where the storage wishes to absorb as much energy as possible. This will result in negative dual prices for the storage and simultaneous charging and discharging will occur if not constrained. Table II shows the simulation result for five trails including upper and lower bounds on the dual and control, where ${\theta}^{*}_{0}$ and $p^{*}_{1}$ is calculated using Gurobi with mixed-integer quadratic programming under default settings, solving (1) using the integer constraint (12) for enforcing constraint (1c). In all test trails, the primal and dual result fall within the calculated range. In terms of computation speed, the proposed method all solves less than 1 millisecond, while the benchmark method using Gurobi may need up to several minutes to solve the problem depending on the problem size (thousands of steps), due to solving a mixed-integer quadratic programming problem.

V Conclusion

This paper proposed a novel algorithm for solving look-ahead control for energy storage. The numerical results illustrate that the algorithm provides computation speed in milliseconds for controlling a single energy storage device over an extended planning period. In future research, we plan on expanding this method to controlling multiple energy storage devices subject to network constraints. Moreover, using the generalized terminal state function we plan on incorporating this algorithm into scenario-based stochastic programming or dynamic programming. In addition, our results connects the optimal control with the Lagrangian multiplier associated with the state-of-charge constraint, which we will further explore to provide key insights into designing future electricity pricing and distributed control schemes.

Appendix A Proof of Theorem 1

We start by showing that when moving to the next control step, the value of the Lagrangian $\theta_{t}$ only changes after $e^{*}_{t}$ is reaching the upper or lower SoC bound, or more specifically:

[TABLE]

Hence, it is trivial to see that if $0\leq\sigma_{\tau}(\theta_{0})\leq E$ $\forall$ $\tau\in[1,t]$ then $\theta_{t-1}=\theta_{0}$, leading $p^{\mathrm{\pi}}_{t}(\theta_{t-1})=p^{\mathrm{\pi}}_{t}(\theta_{0})=p^{*}_{t}$ according to Proposition 1.

We will do the proof separately for three possible cases of $e^{*}_{t}$: 1) $e^{*}_{t}$ never reaches upper or lower bound with all $\underline{\nu}_{t}$ and $\overline{\nu}_{t}$ equal to zero; 2) $e^{*}_{t}$ reached upper bound first; 3) $e^{*}_{t}$ reached lower bound first. These three cases are illustrated in Fig. 1.

A-1 $0<e^{*}_{t}<E$ $\forall$ $t\in[1,T]$

This cover the cases when $e^{*}_{t}$ never reached the upper or lower bound, which from (17) we know $\theta_{0}=\theta_{1}=\dotsc=\theta_{T}$, hence $p^{*}_{t}=p^{\mathrm{\pi}}_{t}(\theta_{0})$ and $e^{*}_{t}=\sigma_{t}(\theta_{0})$ for all $t\in[1,T]$, and in particular for $t=T$ we have

[TABLE]

according to (2d) and the aforementioned result. Since $O_{t}$ and $c_{T}$ are convex, $o_{t}$, $c_{T}$, and $\varphi_{t}$ (inverse of $o_{t}$) are monotonic increasing functions, it follows

[TABLE]

meaning if $x\geq-c_{T}(\sigma_{T}(x))$ then $x\geq\theta_{0}$, thus we proved condition 1-b in the Theorem. Similarly we can prove condition 2-b starting with $x\leq\theta_{0}$. Also it is trivial to see that if any $\sigma_{t}(x)$ goes above $E$ then in this case we know $\sigma_{t}(x)>\sigma_{t}(\theta_{0})$ hence $x>\theta_{0}$, and vice versa for $\sigma_{t}(x)$ goes below 0, hence we proved condition 1-a and 2-a. It is also trivial to see that if $x=c_{T}(\sigma_{T}(x))$, then the KKT condition is satisfied and $x=\theta_{0}$, which proves condition 3. Thus all conditions in this theorem are proved for this this case.

A-2 $\exists\tau\in[1,T]$ s.t. $\overline{\nu}_{\tau}>0$ and $0<e^{*}_{t}<E$ $\forall$ $t\in[1,\tau)$

This covers the cases when $e^{*}_{t}$ reached the upper bound first. An example of this case in shown in Fig. 4.

Now from (19) we can conclude if $\exists\tau\in[1,T]$ s.t. $\overline{\nu}_{\tau}>0$ and $0<e^{*}_{t}<E$ $\forall$ $t\in[1,\tau)$, then the same condition must be satisfied for all $x\geq\theta_{0}$, hence condition 1-a is proved. And from (17) we know after reaching the upper bound, all the following $\theta$ values will be greater than $\theta_{0}$ until $e^{*}_{t}$ reaches the lower bound or till the end of the operation $T$ (i.e., $e^{*}_{t}$ never reaches the lower bound). Without loss of generality, let $\gamma$ be the time that $e^{*}_{t}$ first reaches the lower bound or the end of the operation period, i.e., $e^{*}_{t}>0$ $\forall$ $t\in[1,\gamma)$, it follows

[TABLE]

which means $\sigma_{t}(x)$ either will go below [math] (condition 2-a) or $\sigma_{T}(x)\leq e^{*}_{T}$ which leads to $x\leq c_{T}(\sigma_{T}(x))$ (condition 2-b) according to (19), hence we proved this theorem for this case.

A-3 $\exists\tau\in[1,T]$ s.t. $\underline{\nu}_{\tau}>0$ and $0<e^{*}_{t}<E$ $\forall$ $t\in[1,\tau)$

This covers the case when $e^{*}_{t}$ reaches the lower bound first. This is a mirror proof to the previous case while inverting the upper and lower bound logic, hence this proof is omitted.