1

Real-Time Residential-Side Joint Energy Storage
Management and Load Scheduling with Renewable
Integration
arXiv:1511.08585v2 [cs.SY] 13 Apr 2016

Tianyi Li, Student Member, IEEE, and Min Dong, Senior Member, IEEE

Abstract—We consider joint energy storage management and
load scheduling at a residential site with integrated renewable
generation. Assuming unknown arbitrary dynamics of renewable
source, loads, and electricity price, we aim at optimizing the
load scheduling and energy storage control simultaneously in
order to minimize the overall system cost within a finite time
period. Besides incorporating battery operational constraints and
costs, we model each individual load task by its requested power
intensity and service durations, as well as the maximum and
average delay requirements. To tackle this finite time horizon
stochastic problem, we propose a real-time scheduling and storage control solution by applying a sequence of modification and
transformation to employ Lyapunov optimization that otherwise
is not directly applicable. With our proposed algorithm, we show
that the joint load scheduling and energy storage control can
in fact be separated and sequentially determined. Furthermore,
both scheduling and energy control decisions have closed-form
solutions for simple implementation. Through analysis, we show
that our proposed real-time algorithm has a bounded performance guarantee from the optimal T -slot look-ahead solution
and is asymptotically equivalent to it as the battery capacity
and time period goes to infinity. The effectiveness of joint load
scheduling and energy storage control by our proposed algorithm
is demonstrated through simulation as compared with alternative
algorithms.
Index Terms—Load scheduling, energy storage, renewable
generation, real-time algorithm, stochastic optimization, finite
time horizon

I. I NTRODUCTION
The rising global demand of energy has resulted in high
prices for electricity and also caused the growing environmental concern due to excess carbon emission from power
generation. Integrating renewable energy sources into the grid
system has become a vital green energy solution to reduce
the energy cost and build a sustainable society and economy.
Although promising, renewable energy is often intermittent
and difficult to predict, making it less reliable for both gridlevel operation and as a local energy source for consumers.
Energy storage and flexible loads are considered as two
promising management solutions to mitigate the randomness
of renewable generation, as well as to reduce electricity cost
[2], [3]. In particular, energy storage can be exploited to shift
energy across time, while flexible loads can be controlled to
This work was supported by the National Sciences and Engineering
Research Council (NSERC) of Canada under Discovery Grant RGPIN-201405181.
A preliminary version of this work [1] was presented at the 41st IEEE International Conference on Acoustic Speech and Signal Processing (ICASSP),
Shanghai, China, March 2016.
The authors are with Department of Electrical, Computer, and Software
Engineering, University of Ontario Institute of Technology, Oshawa, Ontario,
Canada L1H 7K4 (email: tianyi.li@uoit.ca, min.dong@uoit.ca).

shift demand across time. For grid operators, they can be
utilized to counter the fluctuation in renewable generation
and to increase reliability. For consumers, energy storage
and load scheduling can provide effective means for energy
management to reduce electricity cost.
As renewable penetration into the power supply increases,
the renewable generation with storage solutions at residential
homes (such as roof-top solar panels) will become increasingly
popular. Thus, developing a cost effective energy storage management system to maximally harness energy from renewable
sources is of critical importance. At the same time, many smart
appliances have been developed, creating more controllable
loads at the consumer side. They can be controlled to benefit
from the dynamic price set at the utility, and help shift the
energy demand from high-peak to low-peak periods to reduce
energy bills. Providing effective management solution that
combines both energy storage and load scheduling will be
the most promising future solution for consumers to reduce
energy costs and is the goal of this paper.
Developing an effective joint energy storage management
and load scheduling solution is important, but faces unique
challenges. For energy storage, the cost reduction by storage
comes with an additional cost from battery degradation due
to charging and discharging; finite battery capacity makes
the storage control decisions coupled over time which are
difficult to optimize. For load scheduling, while minimizing
the electricity cost, it needs to ensure the delay requirements
for each load and for the overall service are met. In particular,
load scheduling decision affects the energy usage and storage
and vise versa. Thus, storage control and load scheduling
decisions are coupled with each other and over time, making
it especially challenging for a joint design.
Energy storage management alone has been considered for
power balancing to counter the fluctuation of renewable generation and increase grid reliability [4]–[6], and for consumers
to reduce electricity costs [7]–[14]. Off-line storage control
strategies for dynamic systems have been proposed [4], [7],
[8]. In these works, renewable energy arrivals are assumed
known ahead of time and the knowledge of load statistics is
assumed. For real-time storage management design, [9] has
formulated the storage management control as a Markov Decision Process (MDP) and solved it by Dynamic Programming
(DP). Lyapunov optimization technique [15] has been recently
employed for designing real-time storage control at either
grid operator side or consumer side under different system
models and optimization goals [5], [6], [10]–[14]. Among
these works, [11], [12] have considered renewable generation
without modeling the battery operational cost. Both renewable

2

generation and battery operation cost have been modeled in
[5], [6], [13], [14]. Except for [14], which considers the storage
management design within a finite time period, all the above
works consider the long-term average system cost.
Load (demand) scheduling through demand side management has been studied by many for shaping the aggregate
load at utility through direct load control [16]–[18] or pricing
optimization [19], [20], and at consumer through load scheduling to reduce electricity bill in response to the dynamic price
[21]–[26]. With the electricity price known ahead of time,
linear programming [21], [22] and DP [23] techniques are
applied for load scheduling. Without assuming known future
prices, MDP formulation has been considered in [24], [25],
and opportunistic load scheduling based on optimal stopping
rule has been proposed in [26]. Combining both utility side
and demand side management is also considered in [27],
[28], where game theoretic approach is applied for distributed
energy management.
Few existing works consider joint optimization of energy
storage management and load scheduling. A joint design has
been developed in [29], in which the electricity price is
assumed to be known ahead of time and the storage model
is simplified without the battery operational cost. Real-time
energy storage management with flexible loads has been considered in [30] and [31]. While [30] focuses on local demand
side, [31] combines both grid operator and demand side
management using distributed storage. In these works, flexible
loads are only modeled in terms of the aggregated energy
request; There is no individual task modeling or scheduling
being conducted. Furthermore, in these works, renewable
generation, loads, and electricity price are assumed to be
independent and identically distributed. With the increasing
penetration of renewable generation and energy storage in the
grid, future energy demand and supply are expected to be quite
dynamic. Renewable generation, loads, and electricity price
may all fluctuate randomly1 with their statistics likely being
non-stationary, making them difficult to predict accurately.
However, most of existing works design solutions assuming
either the future values or statistical knowledge of them to be
known. In addition, although long-term time averaged cost is
typically considered in these existing works, the consumers
may prefer a cost saving solution in a period of time defined
by their own needs. It is important to provide a solution to
meet such need. We aim at proposing a real-time algorithm
for joint energy storage and load scheduling to address these
issues.
In this paper, we consider joint energy storage management
and load scheduling at a residential site equipped with a
renewable generator and a storage battery. For renewable
source, loads, and electricity price, we assume their dynamics
to be arbitrary which can be non-stationary and their statistics
are unknown. For the residential load, we characterize each
individual load task with its own requested power intensity
and service duration, and consider both per load maximum
delay and average delay requirements. For battery storage, we
1 Energy pricing design is considered to respond to the energy demand and
supply status and to shape the load for demand management. As a result, the
real-time price in the future grid could fluctuate much more randomly and
quickly, and likely to have complicated statistical behaviors.

actively model the battery operational constraints and cost due
to charging and discharging activities.
We aim at designing a real-time solution for joint energy
storage management and load scheduling to minimize the
overall system cost over a finite time period, subject to battery
operation and load delay constraints. The interaction of load
scheduling and energy storage, the finite battery capacity, and
finite time period for optimization complicate the scheduling
and energy control decision making over time. To tackle this
difficult stochastic problem, we develop techniques through a
sequence of problem modification and transformation which
enable us to employ Lyapunov optimization to design a realtime algorithm that otherwise is not directly applicable. Interestingly, we show that the joint load scheduling and energy
storage control can be separated and sequentially determined
in our real-time optimization algorithm. Furthermore, both
load scheduling and energy control decisions have closed-form
solutions, making the real-time algorithm simple to implement.
We further show that our proposed real-time algorithm not
only provides a bounded performance guarantee to the optimal
T -slot look-ahead solution which has full future information
available, but is also asymptotically equivalent to it as the
battery capacity and the considered time period for design
go to infinity. Simulation results demonstrate the effectiveness
of joint load scheduling and energy storage control by our
proposed algorithm as compared with alternative solutions
considering neither storage nor scheduling, or storage only.
Different from our recent work [14], in which the energy
storage problem without flexible loads has been considered,
in this work, we explore both energy storage and flexible
load scheduling to reduce system cost. Given the individual
load modeling, load delay requirements imposed, and the load
interaction with energy usage over time, it is highly non-trivial
to formulate the joint design problem, develop techniques for a
real-time solution, and provide performance analysis. Through
our developed techniques, we show that the joint optimization
of load scheduling and storage control can in fact be separated
and sequentially solved. Thus, we are able to obtain the load
scheduling solution in closed-form and apply the result in [14]
for the storage control. Furthermore, we demonstrate that a
lower system cost can be achieved with joint load scheduling
and energy storage control than with just energy storage alone.
Comparing with existing works, our proposed algorithm has
the following features and advantages: 1) The battery storage
operation and associated cost, as well as individual load and its
quality of service, are thoroughly modeled; 2) The algorithm
provides a real-time joint solution for both energy storage
control and load task scheduling; 3) The solution only relies
on the current price, renewable generation, or loads, and does
not require any statistical knowledge of them; 4) The solution
is designed for a specified period of time which may be useful
for practical needs; 5) The solution is provided in closed-form
requiring minimum complexity for practical implementation.
The rest of this paper is organized as follows. In Section II,
we describe the system model. In Section III, we formulate
the joint energy management and load scheduling problem.
In Section IV, we propose a real-time algorithm for our
joint optimization problem. In Section V, we analyze the
performance of algorithm. After presenting our simulation

3

TABLE I
L IST OF MAIN SYMBOLS

Et − Qt
Pt
Grid

Et

− Qt

Battery Dt

+

Loads

Wt

user’s load arriving at time slot t (kWh)

ρt

load intensity (kWh)

λt

load duration (number of slots)

dt

delay incurred for Wt before it is served (number of slots)

dmax
t

maximum delay allowed for Wt before it is served (number
of slots)

dw

average delay of all arrived loads within the To -slot period
(number of slots)

dmax

maximum average delay (number of slots) for the loads within
the To -slot period

Cd (·)

cost function associated with the average delay dw

Et

energy purchased from conventional grid at time slot t (kWh)

Emax

maximum amount of energy that can be bought from the grid
per slot (kWh)

Pt

unit price of buying energy at time slot t ($/kWh)

Pmax

maximum unit energy price ($/kWh)

Pmin

minimum unit energy price ($/kWh)

St

renewable energy harvested at time slot t (kWh)

A. Load Scheduling

Sw,t

amount of renewable energy directly supplied user’s loads to
be served at time slot t (kWh)

Sr,t

amount of renewable energy stored into battery at time slot t
(kWh)

Qt

portion of Et stored into battery at time slot t (kWh)

Rmax

maximum charging amount (kWh) per slot allowed for the
battery

Dt

amount of energy discharged from the battery at time slot t
(kWh)

Dmax

maximum discharging amount per slot allowed from the battery (kWh)

We assume the user has load tasks in various types arriving
over time slots. An example of the scheduling time line of two
loads is shown in Fig. 2. Let Wt denote the load arriving at the
beginning of time slot t. It is given by Wt = ρt λt , where ρt
and λt are the load intensity and duration for Wt , respectively.
We assume λt is an integer multiple of time slots, and the
minimum duration for any load is 1, i.e., λt ∈ {1, 2, . . .}. Let
dmax
denote the maximum delay allowed for Wt before it is
t
served (multiple of time slots), and let dt denote the actual
delay incurred for Wt before it is served. We have

Bt

battery energy level at time slot t (kWh)

Bmin

minimum energy level required in the battery (kWh)

Bmax

maximum energy level allowed in the battery (kWh)

Crc

entry cost for battery due to each charging activity ($)

Cdc

entry cost for battery due to each discharging activity ($)

xe,t

entry cost for battery at time slot t as xe,t = 1R,t Crc +
1D,t Cdc ($)

xu,t

net amount of energy change in battery at time slot t, as
xu,t = |Qt + Sr,t − Dt | (kWh)

xe

average entry cost for battery over the To -slot period ($)

xu

average net amount of energy change in battery over the To slot period (kWh)

Cu (·)

cost function associated with average usage amount xu ($)

α

weight for the cost of scheduling delay

at

energy storage control action vector at time slot t

∆u

desired change of battery energy level within To slots (kWh)

µ

weight for delay related queues in Lyapunov function

studies in Section VI, we conclude our paper in Section VII.
Notations: The main symbols used in this paper are summarized in Table I.
II. S YSTEM M ODEL
We consider a residential-side electricity consuming entity
powered by the conventional grid and a local renewable generator (RG) (e.g., wind or solar generators). An energy storage
battery is co-located with RG to store energy from both power

Sr,t
Renew. St
Gen.

Fig. 1.

−

Contrl.

t
X
τ =0

Sw,t

ρτ 1S,t (dτ )
00
11

00
11
00
11
11111
00000
00
11
0000000000
1111111111
00000
11111
0000000000
1111111111
t

The residential energy storage management system.

sources and supply power to the user. The energy storage
management (ESM) system is shown in Fig. 1. As a part of the
ESM system, a load scheduling mechanism is implemented to
schedule each load task to meet its delay requirements. We
assume the ESM system operates in discrete time slots with
t ∈ {0, 1, · · · }, and all operations are performed per time slot
t. Each component of the EMS system is described below.

dt ∈ {0, 1 . . . , dmax
}, ∀t.
t

(1)

Thus, the earliest serving time duration for Wt is [t, t + λt ],
and the latest serving time duration is [t+dmax
, t+dmax
+λt ].
t
t
We define an indicator function 1S,t (dτ ) , {1 : if t ∈ [τ +
dτ , τ + dτ + λτ ); 0 : otherwise}, for ∀τ ≤ t. It indicates
whether or not the load Wτ is being served at time slot t.
Consider a To -slot period. We define dw as the average delay
of all arrived loads within this To -slot period, given by2
dw ,

To −1
1 X
dτ .
To τ =0

(2)

Besides the per load maximum delay dmax
constraint in (1),
t
we impose a constraint on the average delay dw as
dw ∈ [0, dmax ]

(3)

where dmax is the maximum average delay for the loads
within the To -slot period. It is straightforward to see
that for constraint (3) to be effective, we have dmax ≤
maxt∈[0,To −1] {dmax
}, for ∀t. The average delay dw reflects
t
the average quality of service for the loads within the To slot period. We define a cost function Cd (dw ) associated with
dw . A longer delay reduces the quality of service and incurs
a higher cost. Thus, we assume Cd (·) to be a continuous,
convex, non-decreasing function with derivative Cd′ (·) < ∞.
2 Without

loss of generality, we start the To -period at time slot t = 0.

4

Wt1
...
0

t1

dmax
t1

Wt2
...
...
t
2
dt1

ρt
1110101011
00
00
10
00
00 . . . 11
11
001001. .0. 11
11
00
11
1000
11

ρt1

dt2

Let Bt denote the battery energy level at time slot t. With a
finite capacity, Bt is bounded by

dmax
t2

2

λt1

...
To − 1

λt2

Fig. 2. An example of load scheduling for two arrival loads Wt1 and Wt2 .

B. Energy Sources and Storage
1) Power Sources: The user can purchase energy from the
conventional grid with a real-time price. Let Et denote the
amount of energy bought in time slot t. It is bounded by
Et ∈ [0, Emax ]

(4)

where Emax is the maximum amount of energy that can be
bought from the grid per slot. This amount can be used to
directly supply the user’s loads and/or be stored into the
battery. Let Pt denote the unit price of buying energy at
time slot t. It is bounded as Pt ∈ [Pmin , Pmax ], where Pmin
and Pmax are the minimum and maximum unit energy prices,
respectively. We assume Pt is known to the user at time slot
t and is kept unchanged within the slot duration. The average
cost for the purchased energy from the grid over a To -slot
PTo −1
∆
Et Pt .
period is defined by J = T1o t=0
Renewable generator: An RG is used as an alternative
energy source in the ESM system. Let St denote the amount
of renewable energy harvested at time slot t. We assume St is
first used to supply the loads scheduled to be served at time
slot t. Denote this portion by Sw,t , we have
( t
)
X
Sw,t = min
(5)
ρτ 1S,t(dτ ), St
τ =0

where the first term in (5) represents the total energy over
those scheduled loads that need to be served at time slot t.
The remaining portion of St , if any, can be stored into the
battery. Since there is a cost associated to the battery charging
activity, we use a controller to determine whether or not to
store the remaining portion into the battery. Let Sr,t denote
the amount of renewable energy charged into the battery at
time slot t. It is bounded by
Sr,t ∈ [0, St − Sw,t ] .

(6)

2) Battery Operation: The battery can be charged from
either the grid, the renewable generator, or both at the same
time. Let Qt denote the portion of Et from the grid that is
stored into the battery at time slot t. The total charging amount
at time slot t is bounded by
Qt + Sr,t ∈ [0, Rmax ]

(7)

where Rmax is the maximum charging amount per slot for the
battery. Similarly, let Dt denote the discharging amount from
the battery at time slot t, bounded by
Dt ∈ [0, Dmax ]

(8)

where Dmax is the maximum discharging amount per slot
allowed from the battery. We assume there is no simultaneous
charging and discharging activities at the battery, i.e.,
(Qt + Sr,t ) · Dt = 0.

(9)

Bt ∈ [Bmin , Bmax ]

(10)

where Bmin and Bmax are the minimum energy required and
maximum energy allowed in the battery, respectively. The
dynamics of Bt over time due to charging and discharging
activities are given by3
Bt+1 = Bt + Qt + Sr,t − Dt .

(11)

It is known from battery technology that frequent charging/discharging activities cause a battery to degrade over time
[32]–[34].4 Both the frequency of charging or discharging and
the amount that is charged or discharged affect the battery
lifetime. Given this, we model two types of battery operational
costs associated with the charging/discharging activities: entry
cost and usage cost.
The entry cost is a fixed cost incurred due to
each charging or discharging activity. Define two indicator functions to represent charging and discharging activities as 1R,t , {1 : if Qt + Sr,t > 0; 0 : otherwise} and
1D,t , {1 : if Dt > 0; 0 : otherwise}, respectively. Let Crc
denote the entry cost for each charging activity and Cdc for
that of the discharging activity. Let xe,t denote the entry cost
at time slot t. It is given by xe,t , 1R,t Crc + 1D,t Cdc . We
define thePtime-averaged entry cost over the To -slot period as
To −1
xe,t .
xe , T1o t=0
The usage cost is defined as the cost associated with the
∆
battery charging and discharging amount. Let xu,t = |Qt +
Sr,t − Dt | denote the net amount of energy change in battery
at time slot t due to charging or discharging. From (7) and
(8), it follows that xu,t is bounded by
xu,t ∈ [0, max {Rmax , Dmax }] .

(12)

In general, the battery usage cost is associated with charge
cycles5 . Each charge cycle typically lasts for a period of
time in a day [35]. To approximate this, we consider the
average net amount of energy change inPthe battery over
To −1
the To -slot period, defined as xu , T1o t=0
xu,t . From
(12), it is straightforward to see that xu is bounded by
xu ∈ [0, max {Rmax , Dmax }]. We model the usage cost as
a function of xu , denoted by Cu (xu ). It is known that faster
charging/discharging within a fix period has a more detrimental effect on the life time of the battery. Thus, we assume
Cu (xu ) is a continuous, convex, non-decreasing function with
derivative Cu′ (xu ) < ∞.6
Based on the above, the average battery operational cost
over the To -slot period due to charging/discharging activities
is given by xe + Cu (xu ).
3 We consider an ideal battery model with no leakage of stored energy over
time and full charging/discharging efficiency.
4 The lifetime of battery can be coarsely measured by the number of full
charge cycles.
5 Charging the battery and then discharging it to the same level is considered
a charge cycle.
6 Such a convex cost function has also been adopted in literature [36], [37].

5

C. Supply and Demand Balance
For each load Wτ arrived at time slot τ , if it is scheduled
to be served at time slot t (≥ τ ), the energy supply needs
to meet the amount ρτ scheduled for Wτ . The overall energy
supply must be equal to the total demands from those loads
which need to be served at time slot t. Thus, we have supply
and demand balance relation given by
Et − Qt + Sw,t + Dt =

t
X

ρτ 1S,t (dτ ), ∀t.

(13)

τ =0

A. Problem Modification
Due to the finite battery capacity constraint, the control
actions {at } are coupled over time. To remove the time
coupling, similar to the technique used in our previous work
[14] for energy storage only problem, we remove the finite
battery capacity constraint, and instead we impose a constraint
on the change of battery energy level over the To -slot period.
Specifically, by (11), the change of P
battery energy level over
o −1
(Qt + Sr,t − Dt ).
the To -slot period is BTo − B0 = Tt=0
We now set this change to be a desired value ∆u , i.e.,
To −1
∆u
1 X
(Qt + Sr,t − Dt ) =
.
To t=0
To

III. J OINT E NERGY S TORAGE M ANAGEMENT AND L OAD
S CHEDULING : P ROBLEM F ORMULATION
Our goal is to jointly optimize the load scheduling and
energy flows and storage control for the ESM system to
minimize an overall system cost over the To -slot period.
The loads, renewable generation, and price {Wt , St , Pt } have
complicated statistical behaviors which may be non-stationary
and thus are often difficult to acquire or predict in practice.
In our design, we assume arbitrary dynamics for {Wt , St , Pt }
and do not assume their statistical knowledge being known. We
intend to develop a real-time control algorithm that is capable
to handle such arbitrary and unknown system inputs.
We model the overall system cost as a weighted sum of the
cost from energy purchase and battery degradation, and the
cost of scheduling delay. Define at , [Et , Qt , Dt , Sw,t , Sr,t ]
as the control action vector for the energy flow in the ESM
system at time slot t. Our goal is to find an optimal policy
{at , dt } that minimizes the time-averaged system cost. This
optimization problem is formulated as follows
P1: min

{at ,dt }

J + xe + Cu (xu ) + αCd (dw )

s.t. (1), (3), (4), (6), (9), (13), and
0 ≤ Sr,t + Qt ≤ min {Rmax , Bmax − Bt } (14)
0 ≤ Dt ≤ min {Dmax , Bt − Bmin }
(15)
where α is the positive weight for the cost of scheduling delay;
It sets the relative weight between energy related cost and load
delay incurred by scheduling in the joint optimization.
Note that in P1, {Wt , St , Pt } are random, and their future
values are unknown at time slot t. Thus, P1 is a finite time
horizon joint stochastic optimization problem which is difficult
to solve. Joint energy storage control and load scheduling
complicates the problem, making it much more challenging
than each separate problem alone. The finite battery capacity
imposes a hard constraint on the control actions {at }, making
{at } correlated over time, due to the time-coupling dynamics
of Bt in (11). Furthermore, the finite time horizon problem
is much more difficult to tackle than the infinite time horizon
problem as considered in most existing energy storage works.
New techniques need to be developed for a real-time control
solution.
In the following, we focus on proposing a real-time algorithm to provide a suboptimal solution to P1 with a certain
performance guarantee. To do this, we first modify P1 to allow
us to design a real-time algorithm for joint energy storage
control and load scheduling at every time slot. We later discuss
how our solution can meet the constraints of P1.

(16)

Note that, ∆u is only a desired value we set, which may not be
achieved by a control algorithm at the end of To -slot period.
We will quantify the amount of mismatch with respect to ∆u
under our proposed control algorithm in Section V. By the battery capacity and (dis)charging constraints, it is easy to see that
∆
|∆u | ≤ ∆max = min{Bmax − Bmin , To max{Rmax , Dmax }}.
We now modify P1 to the follow optimization problem
by adding the new constraint (16), and removing the battery
capacity constraint (10)
P2: min

{at ,dt }

J + xe + Cu (xu ) + αCd (dw )

s.t. (1), (3) − (9), (13), (16).
Note that by removing the battery capacity constraint (10),
we remove the dependency of per-slot charging/discharging
amount on Bt in constraints (14) and (15), and replace them
by (7) and (8), respectively.
B. Problem Transformation
In P2, both battery average usage cost Cu (xu ) and scheduling delay cost Cd (dw ) are functions of time-averaged variables, which complicates the problem. Using the technique
introduced in [38], we now transform the problem into one that
only contains the time-average of the functions. Specifically,
we introduce auxiliary variables γu,t and γd,t for xu,t and dt ,
respectively, and impose the following constraints

∆

0 ≤ γu,t ≤ max{Rmax , Dmax }, ∀t
γu = xu

(17)
(18)

0 ≤ γd,t ≤ min{dmax
, dmax }, ∀t
t

(19)

γd = dw
1 PTo −1

(20)

where γi = To τ =0 γi,t , for i = u, d. The above constraints
ensure that each auxiliary variable lies in the same range as
its original variable, and its time average is the
as that
PTsame
o −1
Ci (γi,t )
of its original variable. Define Ci (γi ) , T1o t=0
as the time average of Ci (γi,t ) over To slots, for i = u, d.
Applying (18) and (20) to the objective of P2, and defining
π t , [at , dt , γu,t , γd,t ], we transform P2 into the following
optimization problem
P3: min J + xe + Cu (γu ) + αCd (γd )
{π t }

s.t. (1), (3) − (9), (13), (16) − (20)
where the terms in the objective are all To -slot time-averaged
cost functions. We can show that P2 and P3 are equivalent,

6

i.e., they have the same optimal control solution {a∗t , d∗t } (See
Appendix A).
Although P3 is still difficult to solve, it enables us to design
a dynamic control and scheduling algorithm for joint energy
storage control and load scheduling by adopting Lyapunov
optimization technique [15]. In the following, we propose
our real-time algorithm for P3, and then design parameters
to ensure the proposed solution meets the battery capacity
constraint in the original P1 which is removed in P2.
IV. J OINT E NERGY S TORAGE M ANAGEMENT AND L OAD
S CHEDULING : R EAL -T IME A LGORITHM
By Lyapunov optimization, we first introduce virtual queues
for each time-averaged inequality and equality constraints of
P3 to transform them into queue stability problems. Then, we
design a real-time algorithm based on the drift of Lyapunov
function defined on these virtual queues.
A. Virtual Queues
We introduce a virtual queue Xt to meet constraint (3),
evolving as follows
Xt+1 = max (Xt + dt − dmax , 0) .

(21)

From (2), the above results in dw ≤ dmax + (XTo − X0 )/To .
Thus, formulating the virtual queue Xt in (21) will guarantee
to meet the average delay constraint (3) with a margin (XTo −
X0 )/To . In Section V, we will further discuss this constraint
under our proposed algorithm.
For constraint (16), dividing both sides by To gives the timeaveraged net change of battery energy level per slot being
∆u /To . To meet this constraint, we introduce a virtual queue
Zt , evolving as follows
Zt+1 = Zt + Qt + Sr,t − Dt −

∆u
.
To

∆

∆u
t.
To

L(Θt ) ,


1 2
2
2
Zt + Hu,t
+ µ Xt2 + Hd,t
2

(26)

where µ is a positive weight to adjust the relative importance
of load delay related queues in the Lyapunov function. We
define a one-slot sample path Lyapunov drift as ∆(Θt ) ,
L (Θt+1 ) − L(Θt ), which only depends on the current system
inputs {Wt , St , Pt }.
Instead of directly minimizing the system cost objective in
P3, we consider the drift-plus-cost metric given by ∆(Θt ) +
V [Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]. It is a weighted sum
of the drift ∆(Θt ) and the system cost at time slot t with
V > 0 being the relative weight between the two terms.
Directly using the drift-plus-cost function to determine
control action π t is still challenging. In the following, we use
an upper bound of this drift-plus-cost function to design our
real-time algorithm. The upper bound is derived in Appendix
B as (34). Using this upper bound, we formulate a per-slot
real-time optimization problem and solve it at every time slot
t. By removing all the constant terms independent of control
action πt , we arrive at the following optimization problem
P4 : min Zt [Et + Sr,t + Sw,t − ρt 1S,t (dt )] − |Hu,t | Sw,t
πt

+ Hu,t [γu,t − (Et + Sr,t )] + |Hu,t | ρt 1S,t (dt )
+ µXt dt + µHd,t (γd,t − dt )
+ V [Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]
s.t. (1), (4) − (9), (13), (17), (19).

(23)

u
The linear time function ∆
To t in At is to ensure that constraint
(16) is satisfied. Due to this shift At , the range of Zt is
expanded to the entire real line, i.e., Zt ∈ R for Bt ∈ R+ .
Note that Ao is a design parameter. Later, we design Ao to
ensure that our control solution {at } for the energy flows in
our proposed algorithm satisfies the battery capacity constraint
(10) imposed in P1.
Finally, to meet constraints (18) and (19), we establish
virtual queues Hu,t and Hd,t , respectively, as follows

Hu,t+1 = Hu,t + γu,t − xu,t
Hd,t+1 = Hd,t + γd,t − dt .

Note that Zt and Hu,t are the virtual queues related to
the battery operation, while Xt and Hd,t are those related
to the scheduling delay. Let Θt , [Zt , Hu,t , Xt , Hd,t ] denote
the virtual queue vector. We define the quadratic Lyapunov
function L(Θt ) for Θt as follows

(22)

From Bt in (11) and Zt above, we can show that they are
different by a time-dependent shift as follows
Zt = Bt − At , where At = Ao +

B. Real-Time Algorithm

(24)
(25)

From Lyapunov optimization, it can be shown that satisfying
constraints (3), (16), (18), and (19) is equivalent to maintaining
the stability of queues Xt , Zt , Hu,t , and Hd,t , respectively
[15].

Pt
Note that the term τ =0 ρτ 1S,t (dτ ) in the upper bound (34)
is the total energy demand from the scheduled loads at time
slot t. Since delay dτ for τ ∈ {0, 1, . . . , t − 1} are determined
in previous time slot τ ≤ t − 1 by solving P4, only ρt 1S,t (dt )
is a function of π t at time slot t, and is part of the objective
of P4.
Denote the optimal solution of P4 by π ∗t
,
∗
∗
∗ ∗
, γd,t
]. After regrouping the terms in the objective
[at , dt , γu,t
of P4 with respect to different control variables, we show
that P4 can be separated into four sub-problems to be solved
sequentially and variables in π ∗t can be determined separately.
The steps are described below.
∗
S1) Determine d∗t and γd,t
by solving the following P4a1 and
P4a2 , respectively.

P4a1 : min µdt (Xt − Hd,t ) − ρt 1S,t (dt ) (Zt − |Hu,t |)
dt

s.t. (1).
P4a2 : min µHd,t γd,t + V αCd (γd,t ) s.t. (19).
γd,t

∗
S2) Determine Sw,t
in (5) using d∗t obtained in S1).
∗
∗
S3) Using Sw,t obtained in S2) in (13), determine γu,t
and

7

a∗t by solving the following P4b1 and P4b2 , respectively.
P4b1 : min Hu,t γu,t + V Cu (γu,t ) s.t. (17).
γu,t

P4b2 : min Et (Zt − Hu,t + V Pt ) + Sr,t (Zt − Hu,t )
at

+ V (1R,t Crc + 1D,t Cdc )
s.t.

(4) − (9), (13).

Remark: An important and interesting observation of the
above is that the joint optimization of load scheduling and
energy storage control can in fact be separated: The scheduling
decision d∗t is determined first in P4a1 . Based on the resulting
energy demand in the current time slot t, energy storage
control decision a∗t is then determined in P4b2 . Note that
the two sub-problems are interconnected through the current
virtual queue backlogs Zt and Hu,t for the battery energy
level Bt and battery energy net change xu,t , respectively. The
storage control decision a∗t will further change the battery
energy level and affect Hu,t at the next time slot. Thus,
although the scheduling and storage decisions are separately
determined, they are interconnected through battery energy
level and usage, and sequentially influence each other.
In the following, we solve each subproblem and obtain a
closed-form solution. As a result, the optimal solution π ∗t is
obtained in closed-form.
1) The optimal d∗t : The optimal scheduling delay d∗t for
P4a1 is given below.
Proposition 1: Let ωo , −ρt (Zt − |Hu,t |), ω1 ,
(Xt − Hd,t ).
, µdmax
µ (Xt − Hd,t ), and ωdmax
t
t
(
0 if ωo ≤ ω1
1) If Xt − Hd,t ≥ 0, then d∗t =
1 otherwise;
(
0
if ωo ≤ ωdmax
t
2) If Xt − Hd,t < 0, then d∗t =
max
dt
otherwise.

∗
Lemma 1: The optimal γi,t
, for i = d, u, is given by


if Hi,t ≥ 0

0
′
∗
γi,t = Γi 
(27)
 if Hi,t < −V βi Ci (Γi )

H

Ci′−1 − i,t
otherwise.
V βi

where βu = 1, βd = α/µ, Γu , max{Rmax , Dmax }, and
Γd , min{dmax
, dmax }.
t
Proof: See Appendix D.
3) The optimal a∗t : Once the scheduling decision d∗t for
Wt is determined,
the total energy demand from the scheduled
P
loads, i.e., tτ =0 ρτ 1S,t(d∗τ ), is determined. Given this energy
demand, P4b2 is solved to obtain the optimal control solution
∗
[Et∗ , Q∗t , Dt∗ , Sr,t
] in a∗t . This subproblem for energy storage
and control is essentially the same as in [14], where the energy
storage only problem is considered. Thus, the closed-form
solution can be readily obtained from [14]. Here we directly
state the result.
Pt
∗
Define L∗t ,
τ =0 ρτ 1S,t (dτ ) as the current energy
demand at time slot t. Define the idle state of the battery
as the state where there is no charging or discharging activity. The control solution under this idle state is denoted
id
id
by [Etid , Qid
t , Dt , Sr,t ]. By supply-demand balancing equation
∗
id
id
(13), it is given by Etid = L∗t − Sw,t
, Qid
t = Dt = Sr,t = 0.
Let ξt denote the objective value in P4b2 for the battery being
∗
in the idle state. We have ξt = (L∗t − Sw,t
)(Zt − Hu,t + V Pt ).
′
∗
′
′
′
′
Denote at = [Et , Qt , Dt , Sw,t , Sr,t ]. The optimal control
solution a∗t of P4b2 is given in three cases below.
i ) For Zt − Hu,t + V Pt ≤ 0: The battery is in either charging
or idle state. The solution a′t in charging state is given by

Dt′ = 0,



S ′ = min S − S ∗ , R
max
r,t
w,t
 t
.
∗
∗
′
Q′t = min Rmax − Sr,t
,
E
max − Lt + Sw,t



 ′
∗
′
Et = min L∗t + Rmax − Sw,t
− Sr,t
, Emax
′
If Et′ (Zt −Hu,t +V Pt )+(Zt −Hu,t )Sr,t
+V Crc 1R,t < ξt ,
then a∗t = a′t ; Otherwise, a∗t = aid
.
t

ii ) For Zt − µu Hu,t < 0 ≤ Zt − Hu,t + V Pt : The batProof: See Appendix C.
tery is either in charging, discharging, or idle state. The
are the objective
Remark: Note that ωo , ω1 , and ωdmax
t
or discharging state is given by
solution a′t in
max

 charging
values of P4a1 when dt = 0, 1, and dt , respectively.
∗
′
∗
=
min
L
−
S
,
Dmax
D

t
w,t
t

Furthermore, wo depends on the virtual queue backlogs (Zt

S ′ = min S − S ∗ , R
t
max
and Hu,t ) related to the battery energy level, while ω1 and
r,t
w,t
.
′

depend
on
the
virtual
queue
backlogs
(X
and
H
)
ωdmax
=
0,
Q
t
d,t

t
t

+
 ′  ∗
∗
related to load delay. Proposition 1 shows that the scheduling
− Dmax
Et = Lt − Sw,t
∗
′
decision for load Wt is to either immediately serve it (dt = 0)
If Et′ (Zt − Hu,t + V Pt ) + (Zt − Hu,t )Sr,t
+ V (Crc 1R,t +
∗
′
or delay its serving time (d∗t = 1 or dmax
). This decision
t
Cdc 1D,t ) < ξt , then at = at ; Otherwise, a∗t = aid
t .
depends on whether the battery energy is high enough (so Wt
will be served immediately) or the delays for the scheduled iii ) For 0 ≤ Zt − Hu,t < Zt − Hu,t + V Pt : The′ battery is in
either discharging or idle
 state. Thesolution at in dischargloads so far are low enough (so Wt will be delayed). When the
′
∗
∗

Dt = min Lt − Sw,t , Dmax
load is delayed to serve, the delay should be either minimum
′
.
ing state is given by Sr,t
= Q′t = 0,
or maximum depending on the existing scheduling delays of

 ∗
+

′
∗
the past loads.
Et = Lt − Sw,t − Dmax
′
∗
∗
If
E
(Z
−
H
+
V
P
) + V Cdc 1D,t < ξt , then a∗t = a′t ;
t
u,t
t
t
2) The optimal γd,t and γu,t : Since Cd (·) and Cu (·) are
∗
id
Otherwise, at = at .
both convex, the objectives of P4a2 and P4b1 are convex. Let
′−1
′
Ci (·) denote the first derivative of Ci (·), and Ci (·) denote
In each case above, the cost of charging or discharging is
the inverse function of Ci′ (·), for i = d, u. We obtain the compared with the cost ξt of being in an idle state, and the
∗
, i = d, u, for P4a2 and P4b1 as follows. control solution of P4b2 is the one with the minimum cost.
optimal solutions γi,t
The condition for each case depends on Zt , Hu,t and Pt , where

8

Algorithm 1 Real-Time Joint Load Scheduling and Energy
Storage Management
Set the desired value of ∆u . Set Ao and V as in (28) and (29),
respectively.
At time slot 0: Set Z0 = X0 = Hu,0 = Hd,0 = 0.
At time slot t: Obtain the current values of {Wt , St , Pt }
1) Load scheduling: Determine d∗t according to Proposi∗
tion 1 and γd,t
according to (27), respectively.
2) Energy Storage Control:
∗
a) Renewable contribution: Determine Sw,t
in (5)
using d∗t obtained above.
∗
b) Energy purchase and storage: Determine γu,t
according to (27) and a∗t according to Cases i)-iii) in
Section IV-B3.
3) Updating virtual queues: Use π ∗t to update Bt based on
(11), and Xt , Zt , Hu,t , Hd,t based on (21)−(25).
Zt and Hu,t are rated to battery energy level Bt and usage
cost xu,t , respectively. Thus, Cases i)-iii) represent the control
actions at different battery energy levels (i.e., low, moderate,
or high) and electricity prices.
4) Feasibility of a∗t to P1: Recall that we have removed
the battery capacity constraint (10) when modifying P1 to P2.
Thus, this constraint is no longer imposed in P4b2 , and our
real-time algorithm may not provide feasible control solutions
{a∗t } to P1. To ensure the solution is still feasible to the
original problem P1, we design our control parameters Ao and
V . The result readily follows [14, Proposition 2]. We omit the
details and only state the final result below.
Proposition 2: For the optimal solution a∗t of P4b2 , the
resulting Bt satisfies the battery capacity constraint (10), and
{a∗t } is feasible to P1, if Ao in (23) is given by
(
A′o
if ∆u ≥ 0
Ao =
(28)
A′o − ∆u if ∆u < 0
where A′o = Bmin + V Pmax + V Cu′ (Γu ) + Γu + Dmax +
and V ∈ [0, Vmax ] with
Vmax =

∆u
To ,

Bmax − Bmin − Rmax − Dmax − 2Γu − |∆u |
. (29)
Pmax + Cu′ (Γu )

C. Discussions
We summarize the proposed real-time joint load scheduling
and energy storage control in Algorithm 1. Due to the separation of joint optimization, the algorithm provides a clear
sequence of control decisions at each time slot t.
Recall that we modify the original joint optimization problem P1 to P3, and apply Lyapunov optimization to propose a
real-time algorithm for P3 which is to solve the per-slot optimization problem P4. We then design system parameters Ao
and Vmax to ensure our solution satisfies the battery capacity
constraint, which is removed when we modify P1 to P4. As
a result, our proposed solution by Algorithm 1 is feasible to
P1. Furthermore, we have the following discussions.
1) In modifying P1, we remove the battery capacity constraint
(10) and instead impose a new constraint (16) on the overall
change of battery energy level over To slots to be ∆u . Note

that this constraint is set as a desired outcome, i.e., ∆u is a
desired value. The actual solution a∗t in the proposed algorithm
may not satisfy this constraint at the end of the To -slot period,
and thus may not be feasible to P3. Nonetheless, setting Ao
and V as in (28) and (29) guarantees that {a∗t } satisfy the
battery capacity constraint (10) and therefore are feasible to
P1.
2) In designing the real-time algorithm by Lyapunov framework, virtual queue Xt in (21) which we introduce for the
average delay constraint (3) can only ensure the constraint
is satisfied with a margin as indicated below (21). As a
result, constraint (3) can only be approximately satisfied.
However, this relaxation is mild in practice for the average
delay performance. Note that the per-load maximum delay
constraint (1) is strictly satisfied for d∗t by Algorithm 1. We
will show in simulation that the achieved average delay dw by
Algorithm 1 in fact meets constraint (3).
3) We point out that the load scheduling and energy storage
control decisions are provided in closed-form by Algorithm 1.
Thus, the algorithm is particularly suitable for real-time implementation with a constant computational complexity O(1).
Furthermore, no statistical assumptions on the loads, renewable source, and pricing {Wt , St , Pt } are required in the algorithm. They can be non-stochastic or stochastic with arbitrary
dynamics (including non-stationary processes). This allows
the algorithm to be applied to general scenarios, especially
when these statistics are difficult to predict. Finally, despite
that Algorithm 1 is a suboptimal solution for P1, we will
show in the following that it provides a provable performance
guarantee.
V. P ERFORMANCE A NALYSIS
In this section, we analyze the performance of Algorithm 1
and discuss the mismatch involved in some constraints as a
result of the real-time algorithm design.
A. Algorithm Performance
To evaluate the proposed algorithm, we consider a T -slot
look-ahead problem. Specifically, we partition To slots into T
frames with To = M T , for M, T ∈ N+ . For each frame, we
consider the same problem as P1 but the objective is the T slot averaged cost within the frame and the constraints are all
related to time slots within the frame. In addition, we assume
{Wt , St , Pt } for the entire frame are known beforehand. Thus,
the problem becomes a non-causal static optimization problem
and we call it a T -slot lookahead problem. Let uopt
m be the
corresponding minimum objective value achieved by a T -slot
look-ahead optimal solution over the mth frame. We intend to
bound the performance of Algorithm 1 (with no knowledge of
future values of {Wt , St , Pt }) to the optimal T -slot lookahead
performance (with full knowledge of {Wt , St , Pt } in a frame).
We denote the objective value of P1 achieved by Algorithm 1 over To -slot period by u∗ (V ), where V is the weight
value used in Algorithm 1. The following theorem provides a
bound of the cost performance under our proposed real-time
algorithm to uopt
m under the T -slot lookahead optimal solution.
Theorem 1: Consider {Wt , St , Pt } being any arbitrary processes over time. For any M, T ∈ N+ satisfying To = M T ,

1
M

uopt
m ≤

m=0

L(Θ0 ) − L(ΘTo )
GT
+
V
V To

C ′ (Γu )(Hu,0 − Hu,To ) + αCd′ (Γd )(Hd,0 − Hd,To )
(30)
+ u
To
where G is given in (43) and the upper bound is finite. In
particular, as To → ∞, we have
1
lim u∗ (V ) − lim
To →∞
To →∞ M

M−1
X
m=0

uopt
m

GT
.
≤
V

(31)

B. Design Approximation
1) Average scheduling delay dmax : Recall that, by Algorithm 1, using the virtual queue Xt in (21), the average delay
constraint (3) is approximately satisfied with a margin, i.e.,
dw ≤ dmax + ǫd , where ǫd , (XTo − X0 )/To is the margin.
We now bound ǫd below.
Proposition 3: Under Algorithm 1, the margin ǫd for constraint (3) is bounded as follows
s
L(Θo ) |X0 |
2G
+
+
.
(32)
|ǫd | ≤
µTo
µTo
To
Proof: See Appendix F.
Proposition 3 indicates that the margin ǫd → 0 as To → ∞.
Thus, the average delay is asymptotically satisfied. Note that,
for X0 = 0, ǫd ≥ 0. If X0 > 0, it is possible that ǫd < 0 and
constraint (3) is satisfied with a negative margin. However, this
will drive d∗t to be smaller which may cause higher system
cost. Thus, we set X0 = 0 in Algorithm 1.
2) Mismatch of ∆u : In our design, we set ∆u to be a
desired value for the change of battery energy level over
To -slot period as in new constraint (16). This value may
not bePachieved by Algorithm 1. Define the mismatch by
o −1
ǫu , τT=0
(Qτ +Sr,τ −Dτ )−∆u . The bound for ǫu follows
the result in [14, Proposition 3] and is shown below.
(33)

Note that Vmax in (29) increases as |∆u | decreases, and a
larger Vmax is preferred for better performance by Theorem 1.
PTo −1
To → ∞, constraint (16) becomes limTo →∞ T1
τ =0 (Qτ +
o
− Dτ ) = 0 as in the infinite time horizon problem [39].

7 As

Sr,τ

0.05

0

5

10

15

20

15

20

25

Hour
0.2

Sm

St

0.1

Sl

Sh

0
0

5

10

25

Hour
0.2

0

Wm

Wt

0.1

Wl
Wh

0

5

10

15

20

25

Hour

Proof: See Appendix E.
Remark: Theorem 1 shows that our proposed algorithm is
able to track the “ideal” T -slot lookahead optimal solution
with a bounded gap, for all possible M and T . Also, for
the best performance, we should always choose V = Vmax .
The bound in (31) gives the asymptotic performance as To
increase. Since Vmax in (29) increases with Bmax , it follows
that Algorithm 1 is asymptotically equivalent to the optimal
T -slot look-ahead solution as the battery capacity and To go
to infinity. Note that, as To → ∞, P1 becomes an infinite time
horizon problem with average sample path cost objective7. The
bound in (31) provides the performance gap of long-term timeaveraged sample-path system cost of our proposed algorithm
to the T -slot look-ahead policy.

|ǫu | ≤ 2Γu + Rmax + V Pmax + V Cu′ (Γu ) + Dmax .

Pt

0.1

Energy (kWh)

u∗ (V ) −

M−1
X

Load (kWh)

the To -slot average system cost under Algorithm 1 is bounded
by

Unit price (CAD)

9

Fig. 3. Mean load W t , mean renewable generation S t , and real-time

price Pt over 24 hours.

Thus, a smaller |∆u | is preferred. Note also that our simulation
study shows that the actual mismatch ǫu is much smaller than
this upper bound.
VI. S IMULATION R ESULTS
We set each slot to be 5 minutes and consider a 24-hour
duration. Thus, we have To = 288 slots for each day. We
assume Pt , St and Wt do not change within each slot. We
collect data from Ontario Energy Board [40] to set the price Pt .
As shown Fig. 3 top, it follows a three-stage price pattern as
{Ph , Pm , Pl } = {$0.118, $0.099, $0.063} and is periodic every 24 hours. We assume {St } to be solar photovoltaic energy.
It is a non-stationary process, with the mean amount S t =
E[St ] changing periodically over 24 hours, and having threestage values as {S h , S m , S l } = {1.98, 0.96, 0.005}/12 kWh
per slot8 and standard deviation as σSi = 0.4S i , for i =
h, m, l, as shown in Fig. 3 middle. We assume the load {Wt }
is a non-stationary process, having three-stage mean values
W t = E[Wt ] as {W h , W m , W l } = {2.4, 1.38, 0.6}/12 kWh
per slot9 with standard deviation as σWi = 0.2W i , for
i = h, m, l, as shown in Fig. 3 bottom. For each load Wt , we
generate λt from a uniform distribution with interval [1, 12],
and ρt = Wt /λt . We set dmax
in (1) to be identical for all t.
t
We set the battery related parameters as follows: Rmax =
Dmax = 0.165 kWh10 , Crc = Cdc = 0.001, Bmin = 0, and the
battery initial energy level B0 = 0. Unless specified, we set
Bmax = 3 kWh. We set Emax = 0.3 kWh. Also, we set the
weights α = 1 and µ = 1 as the default values. Since Vmax
increases as |∆u | decreases, to achieve best performance11, we
set ∆u = 0 and V = Vmax .
We consider an exemplary case where the battery usage
cost and the delay cost are both quadratic functions, given
8 We set the mean renewable amount to be {S , S , S }
=
m
h
l
{1.98, 0.96, 0.005} kWh per hour. Converting to per slot energy amount
with 5 minutes duration, we have {1.98, 0.96, 0.005}/12 kWh per slot. We
set these values based on that the average energy harvested by photovoltaic
is about 20 kWh within a day for a residential home.
9 Similar as S , the mean load is converted from energy amount per hour to
t
per slot. The mean values are set based on that the each residential household
on average consumes 1 ∼ 2 kWh per hour.
10 Assuming a household consumes 1 ∼ 2 kWh per hour on average, we
set the values of Rmax and Dmax to be in the same range for load supply.
Thus, per slot, we have Rmax = Dmax = 1.98 kWh/12 = 0.165 kWh.
11 A detailed study of ∆ can be found in [14].
u

10
2

by Cu (xu ) = ku xu 2 and Cd (dw ) = kd dw . The constant
ku > 0 is a battery cost coefficient depending on the battery
characteristics. We set it as ku = 0.2.12 The constant kd > 0 is
a normalization factor based on the desired maximum average
delay dmax in (3). It is set as kd = 1/(dmax )2 , such that the
cost for an average delay dw = dmax is normalized to 1.13 The
∗
′
optimal γi,t
 can be determined with Ci (Γi ) = 2ki Γi ,
 in (27)
H

and Ci′−1 − V i,t
βi

Loads

Serving
Delay
Empty

H

= − 2βi ki,ti V , for i = u, d.

} ρt

A. An Example of Load Scheduling
305

12 The value of k depends on the type of battery and current battery
u
technology. We have tested a range of values for ku and studied the effect of
ku on the system cost performance in our previous work [14].
13 The normalization by k ensures the delay cost is properly defined.
d
Intuitively, a longer delay is encouraged in order to allow more energy cost
saving, as long as the desired maximum average delay is satisfied. Thus, the
cost Cd (dw ) should not be affected by the absolute delay value dw , but rather
2
its relative value to dmax , i.e., dw /(dmax )2 .
14 We set the value of α to ensure that the weighted delay cost αC (d ) is
d w
comparable to the other two costs in the objective of P1. Thus, both delay and
energy cost are active factors in making the storage and scheduling decision
in simulation.

315

320

325

330

335

340

345

350

Fig. 4. A trace of load scheduling results over (dmax
= dmax = 18,
t
α = 0.005).
5.2

×10 -3

α=0.001
α=0.005

5
4.8

B. Effect of Scheduling Delay Constraints
1) Effect of dmax
and dmax : We study how the average
t
system cost objective of P1 under our proposed algorithm
varies with different delay requirements. We set dmax
= dmax ,
t
max
∀t, and plot the average system cost vs. dt
in Fig. 5, for
different values of weight α in the cost objective. As can
be seen, the system cost decreases as dmax (dmax
) increases.
t
This shows that relaxing the average delay constraint gives
more flexibility to load scheduling, where each load can be
scheduled at lower electricity price, resulting in lower system
cost. This demonstrates that flexible load scheduling is more
beneficial to the overall system cost. In addition, we see that
a larger value α gives more weight on minimizing the delay
in the objective, resulting in a higher system cost.
Next, we study the effect of load delay constraints on the
monetary cost, i.e., the cost of energy purchasing and battery
degradation, given by J + xe + Cu (xu ) in the objective of
P1. The monetary cost indicates how much saving a consumer
could actually have, by allowing longer service delay. In Fig. 6,
we plot the monetary cost vs. dmax for dmax
= 216. We
t
see a clear trade-off between the monetary cost and the load
delay. The trade-off curve can be used to determine the desired
operating point. For comparison, we also consider the case
in which all loads are served immediately after arrival, i.e.,
dmax
= 0. This is essentially the case with only energy storage
t

310

Time slot t

Avg. sys. cost

In Fig. 4, we show a fraction of the load scheduling results
by Algorithm 1, where we set dmax
= dmax = 18 slots and
t
14
α = 0.005. Each horizontal bar represents a scheduled load
Wt with the width representing the load intensity ρt and the
length representing the total duration from arrival to service
being completed. For a delayed load, the delay is indicated in
different color before the load is scheduled. In this example,
we see some loads are immediately scheduled, while others
are scheduled at dmax
. The total energy demand at each time
t
slot is the vertical summation over all loads that are scheduled
in this time slot.

4.6
4.4
4.2
4
3.8
3.6
3.4
60

80

100

120

140

160

180

200

220

dmax
(dmax )
t
Fig. 5.

Average system cost vs. dmax
(dmax = dmax
).
t
t

but no load scheduling. Thus, the monetary cost is independent
of dmax . We see a substantial gap between the two curves and
the gap increases with dmax . This clearly shows the benefit of
joint load scheduling and energy storage management.
2) Average delay dw : We now study the average delay dw
achieved by our proposed algorithm vs. dmax for various dmax
t
in Fig. 7. As we see, the actual averaged delay dw increases
with the average delay dmax requirement. This is because, with
a more relaxed constraint on the delay, loads can be shifted
to a later time in order to reduce the system cost, resulting in
larger average delay. However, the increase is sublinear with
respect to dmax . Similarly, we observe that increasing the per
load maximum delay dmax
increases dw . Finally, recall that
t
we study the margin for average delay constraint (3) under our
proposed algorithm in Proposition 3. To see how the resulting
dw meets constraint (3), we plot the line dmax in Fig. 7. As
we see, dw is below dmax for all values of dmax and dmax
.
t
3) Effect of µ: Weight µ is used to control the relative
importance of virtual queues related to the battery and those
to delay in Lyapunov function L(Θt ) in (26) and Lyapunov
drift ∆(Θt ). For Lyapunov drift ∆(Θt ), if µ is large, the two
queues Xt and Hd,t related to the load delay will dominate the
drift. This will affect the drift-plus-cost objective considered
in our proposed algorithm and thus the performance. To study
the effect of µ on the performance, in Fig. 8, we evaluate the
average system cost for different values of µ. We see that a
lower system cost is achieved by smaller value of µ. This is

11

−3

×10 -3

5.4
5.2
5

3.7

dmax
=0
t
dmax
=216
t

3.6

3.5

4.8
4.6

µ=1
µ=0.00001

4.4
4.2
4

3.4

3.8
3.6
60

3.3
0

50

100

150

d
Fig. 6.

x 10

3.8

Avg. sys. cost

Avg. monetary cost ($/slot)

3.9

200

250

80

100

120

140

160

180

200

220

dmax (dmax)

max

t

Monetary cost vs. dmax (α = 0.005).

Fig. 8.

Average system cost vs. dmax (dmax
= dmax ).
t
4.8

150

×10 -3

4.6

max

d

Avg. sys. cost

4.4

dw

100

max
dt =72
dmax=144
t
dmax=216
t
max
dt =288

50

0

0

50

100

150

200

250

Algorithm 1
4
3.8
3.6
3.4
3.2
1

300

dw vs. dmax (α = 0.005).

because, with smaller µ, the drifts of Xt related to delay is
less significant in the overall drift. This allows wider difference
between dw and dmax , i.e., smaller dw and lower delay cost.

C. Performance vs. Battery Capacity
We consider two other algorithms for comparison: A) No
storage or scheduling: In this case, neither energy storage
nor load scheduling is considered. Each load is served immediately using energy purchased from the conventional grid
and/or renewable generator. B) Storage only: In this method,
only battery storage is considered but every load is served
immediately without a delay. This is essentially the algorithm
provided in [14].
In Fig. 9, we compare our proposed algorithm to the above
two alternative algorithms under various battery capacity
Bmax . Since algorithm A does not use a battery, the system
cost is unchanged over Bmax and is 0.01. We do not plot
the curve as the cost is much higher than the rest of two
algorithms we considered. For algorithm B and our proposed
Algorithm 1, as can be seen, the system costs reduces as Bmax
increases. This is because a larger battery capacity allows
charging/discharging to be more flexible based on the current
demand and electricity price, resulting in a lower system cost.
Comparing the two, we see that joint load scheduling and
energy storage control provides further reduction in system
cost.

2

3

4

5

6

7

8

9

10

B max (kWh)

max

d
Fig. 7.

Storage only

4.2

Fig. 9.

Average system cost vs. Bmax (α = 0.001).

VII. C ONCLUSION

In this work, we have considered joint energy storage
management and load scheduling for the ESM system, where
renewable source, loads, and price may be non-stationary and
their statistics are unknown. For load scheduling, we have
characterized each load task by its power intensity and service
duration, and have considered the maximum per load delay
and maximum average delay requirements. For storage control,
our storage model includes details of the battery operational
constraints and cost. Aiming at minimizing the overall system
cost over a finite period of time, we have designed a real-time
algorithm for joint load scheduling and energy storage control,
where we have provided a closed-form per slot scheduling
and energy storage decisions. As a result, we have shown that
the joint load scheduling and energy storage control can in
fact be separately and sequentially determined in our real-time
algorithm. Furthermore, we have shown that our proposed realtime algorithm has a bounded performance guarantee from
an optimal T -slot look-ahead solution and is asymptotically
equivalent to the optimal T -slot look-ahead solution, as the
battery capacity and time period go to infinity. Simulation
results have demonstrated the gain of joint load scheduling
and storage control provided by our proposed algorithm over
other real-time schemes which consider neither storage nor
scheduling, or with storage only.

12

P ROOF

OF

A PPENDIX A
E QUIVALENCE OF P2 AND P3

Proof: The proof follows the same approach as in [14,
Lemma 1]. Let uo2 and uo3 denote the minimum objective values
of P2 and P3, respectively. Since the optimal solution of P2
satisfies all constraints of P3, it is a feasible solution of P3.
Thus, we have uo3 ≤ uo2 . By Jensen’s inequality and convexity
of Ci (·) for i = d, u, we have Cd (γd ) ≥ Cd (γd ) = Cd (dw )
and Cu (γu ) ≥ Cu (γu ) = Cu (xu ). This means uo3 ≥ uo2 .
Hence, we have uo2 = uo3 and P3 and P2 are equivalent.

U PPER B OUND

A PPENDIX B
D RIFT-P LUS -C OST F UNCTION

ON

The following lemma presents an upper bound on the drift
∆(Θt ).
Lemma 2: The one-slot Lyapunov drift ∆(Θt ) is upper
bounded by
!
t
X
∆u
ρτ 1S,t (dτ ) −
∆(Θt ) ≤ Zt Et + Sr,t + Sw,t −
To
τ =0

+ Hu,t γu,t − Hu,t (Et + Sr,t ) + µXt (dt − dmax ) + G
!
t
X
− |Hu,t | Sw,t −
ρτ 1S,t (dτ ) + µHd,t (γd,t − dt ) (34)
τ =0


2 
2 
∆u
u
where G = 12 max
+
,
D
+
Rmax − ∆
max
To
To

µ
1
2
2
max 2
) , (dmax
− dmax )2 +
t
2 max{Rmax , Dmax } + 2 max (d
µ max 2
) .
2 (dt
Proof: From the definition of ∆(Θt ), we have

∆(Θt ) , L(Θt+1 ) − L(Θt )
1 2
2
2
= [Zt+1
− Zt2 + (Hu,t+1
− Hu,t
)]
2

µ 2
2
2
(Xt+1 − Xt2 + Hd,t+1
− Hd,t
)
(35)
+
2
2
where from queue (22), Zt+1
− Zt2 can be presented by


2
Zt+1
− Zt2
∆u
= Zt Qt + Sr,t − Dt −
2
To
u 2
(Qt + Sr,t − Dt − ∆
To )
+
.
(36)
2
u
Note that from (16), we have ∆
To ≤ max{Rmax , Dmax }. For
u 2
a given value of ∆u , byn (7) and (8), (Qt + Sr,t − Dt − o∆
To )

∆u 2
u 2
is upper bound by max (Rmax − ∆
. By
To ) , (Dmax + To )
the supply-demand balance (13), the first term on RHS of (36)
can be replaced by


∆u
Zt Qt + Sr,t − Dt −
To
!
t
X
∆u
ρτ 1S,t (dτ ) −
= Zt Et + Sr,t + Sw,t −
. (37)
To
τ =0
2
2
From (24), Hu,t+1
− Hu,t
in (35) can be presented by

2
2
Hu,t+1
− Hu,t
(γu,t − xu,t )2
= Hu,t (γu,t − xu,t ) +
. (38)
2
2
Note that from (12) and (17), the second term of RHS in (38)
2
2
is upper bounded by (γu,t − xu,t )2 ≤ max{Rmax
, Dmax
}.

We now find the upper bound for −Ht xu,t in the first term
on RHS of (38). By the supply-demand balance (13), −Ht xu,t
can be replaced by
−Ht xu,t = −Hu,t (|Qt + Sr,t − Dt |)
!
t
X
ρτ 1S,t(dτ ) . (39)
= −Hu,t Et + Sr,t + Sw,t −
τ =0

The upper bound of −Ht xu,t in (39) is obtained as follows.
1 ) For Hu,t ≥ 0: We have
!
t
X
ρτ 1S,t (dτ )
− Hu,t Et + Sr,t + Sw,t −
τ =0
t
X

≤ −Hu,t Et + Sr,t + Sw,t −

!

ρτ 1S,t (dτ )

τ =0

= −Hu,t (Et + Sr,t ) − Hu,t Sw,t −

t
X

!

ρτ 1S,t (dτ ) .

τ =0

2 ) For Hu,t < 0: We have
− Hu,t

Et + Sr,t + Sw,t −

t
X

ρτ 1S,t (dτ )

τ =0

≤ −Hu,t |Et + Sr,t | + Sw,t −

t
X

!

ρτ 1S,t (dτ )

τ =0

≤ −Hu,t Et + Sr,t +

t
X

ρτ 1S,t (dτ ) − Sw,t

τ =0

= −Hu,t (Et + Sr,t ) + Hu,t Sw,t −

t
X

!

!
!

ρτ 1S,t (dτ ) .

τ =0

Combine the above cases for Hu,t , we have −Ht xu,t in (39)
upper bounded by
!
t
X
− Hu,t Et + Sr,t + Sw,t −
ρτ 1S,t (dτ )
τ =0

≤ −Hu,t (Et + Sr,t ) − |Hu,t | Sw,t −

t
X

!

ρτ 1S,t (dτ ) .

τ =0
2

2
From (21), we have Xt+1
≤ (Xt + dt − dmax ) . Thus,
2
2
(Xt+1 − Xt ) in (35) is bounded by
2
Xt+1
− Xt2
1
2
≤ Xt (dt − dmax ) + (dt − dmax )
(40)
2
2
where by (1), the last term
n in RHS of (40) is upperobounded
2
2
max 2
by (dt − d
) ≤ max (dmax ) , (dmax
− dmax ) .
t
2
2
From (25), Hd,t+1
− Hd,t
in (35) can be presented by
2
2
Hd,t+1
− Hd,t
(γd,t − dt )2
= Hd,t (γd,t − dt ) +
2
2
1 max 2
(41)
≤ Hd,t (γd,t − dt ) + (dt )
2
where the last inequality is derived from the bounds of γd in
(19) and dt in (1). We give the upper bond of (35) as follows

∆(Θt ) , L(Θt+1 ) − L(Θt )
!
t
X
∆u
ρτ 1S,t (dτ ) −
≤ Zt Et + Sr,t + Sw,t −
To
τ =0

13

− |Hu,t | Sw,t −

t
X

!

ρτ 1S,t (dτ )

τ =0

+ Hu,t γu,t − Hu,t (Et + Sr,t ) + µXt (dt − dmax )

(42)

where G includes all constant terms from the upper bounds
of (36), (38), (40) and (41), and is defined as
(
2 
2 )
1
∆u
∆u
G , max
Rmax −
, Dmax +
2
To
To
µ
1
2
2
max{Rmax
, Dmax
} + (dmax
)2
2
2 t

µ
− dmax )2 .
+ max (dmax )2 , (dmax
t
2

+

P ROOF

A PPENDIX E
P ROOF OF T HEOREM 1

+ µHd,t (γd,t − dt ) + G

(43)

A PPENDIX C
OF P ROPOSITION 1

Proof: To determine the optimal scheduling delay d∗t , we
need to compare the objective values of P4a1 under all serving
options. The optimal delay d∗t is the one that achieves the
minimum objective value. Because dt · 1S,t(dt ) = 0, we have
the following cases:
1) If the load is immediately served, we have 1S,t (dt ) = 1
and dt = 0. The objective value becomes ωo ;
2) If the load is delayed, we have dt > 0 and 1S,t (dt ) = 0.
The objective function is reduced to µdt (Xt −Hd,t ). For
Xt − Hd,t ≥ 0, the objective value is ω1 ; Otherwise, the
.
value is ωdmax
t
, we obtain d∗t .
Comparing ωo to ω1 or ωdmax
t

Proof: A T -slot sample path Lyapunov drift is defined by
∆T (Θt ) , L(Θt+T ) − L(Θt ). We upper bound it as follows

2
2
2
Zt+T
− Zt2 + Hu,t+T
− Hu,t
∆T (Θt ) =
2


2
2
2
µ Xt+T − Xt2 + Hd,t+T
− Hd,t
+
2

t+T
X−1 
∆u
≤ Zt
Qτ + Sr,τ − Dτ −
To
τ =t
"t+T −1 
#2
∆u
1 X
Qτ + Sr,τ − Dτ −
+
2 τ =t
To
"
#2
t+T
t+T −1
X−1
1 X
+ Hu,t
(γu,τ − xu,τ ) +
(γu,τ − xu,τ )
2 τ =t
τ =t
"t+T −1
#2
t+T
X−1
µ X
max
max
(dt − d
)
+ Xt
(dτ − d
)+
2 τ =t
τ =t
"t+T −1
#2
t+T
X−1
µ X
+ Hd,t
(γd,τ − xd,τ ) +
(γd,τ − xd,τ )
2 τ =t
τ =t

t+T
X−1 
∆u
≤ Zt
Qτ + Sr,τ − Dτ −
To
τ =t
+ Hu,t
+ Hd,t

A PPENDIX D
P ROOF OF L EMMA 1

t+T
X−1

(γτ − xu,τ ) + Xt

τ =t
t+T
X−1

t+T
X−1

(dτ − dmax )

τ =t

(γd,τ − xd,τ ) + GT 2

(44)

τ =t

Proof: Since µ, α and V are all positive weights, and
Ci (γt )’s are assumed to be continuous, convex and nondecreasing functions with respect to γi,t with maximum
∗
derivatives Ci′ (Γi ) < ∞, for i = d, u, the optimal γi,t
’s
are determined by examining the derivatives of the objective
functions of P4a2 and P4b1 . Note that, given γu,t in (17) and
γd,t in (19), Ci′ (γi,t ) ≥ 0 and increases with γi,t for i = d, u.
For βi defined in Lemma 1, we have
1 ) For Hi,t ≥ 0: We have µd Hd,t + V αCd′ (γd,t ) > 0 and
Hu,t +V Cu′ (γu,t ) > 0. Thus, the objectives of P4a2 and P4b1
are both monotonically increasing functions, and the minimum
∗
values are obtained with γi,t
= 0 for i = d, u.
′
2 ) For Hi,t < −V βCi (Γi ): Since V Ci′ (Γi ) ≥ V Ci′ (γi,t )
for i = d, u, we have µd Hd,t + V αCd′ (γd,t ) < 0 and
Hu,t + V Cu′ (γu,t ) < 0. The objectives of P4a2 and P4b1
are both monotonically decreasing functions. From (19), the
∗
minimum objective value of P4a2 is reached with γd,t
= Γd
max max
where Γd , min{dt , d
}; The minimum objective value
∗
of P4b1 is reached with γu,t
= Γu where by (17), we have
Γu , min{Rmax , Dmax };
∗
3 ) For −V βi Ci′ (Γi ) ≤ Hi,t ≤ 0: In this case, γd,t
and
′
∗
γu,t are the roots of µd Hd,t + V αCd (γd,t ) = 0 and Hu,t +


where G is defined in Lemma 2.
Assume To = M T . We consider a per-frame optimization
problem below, with the objective of minimizing the timeaveraged system cost within the mth frame of length T time
slots.
1
{at ,γt } T

Pf : min

(m+1)T −1

X

[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]

t=mT

s.t (1), (4), (6), (9), (13) − (15), (17) − (20).
We show that Pf is equivalent to P1 in which To is replaced
by T . Let ufm denote the minimum objective value of Pf .
The optimal solution of P1 satisfies all constraints of Pf and
opt
therefore is feasible to Pf . Thus, we have ufm ≤ um . By
Jensen’s inequality and convexity of Ci (·) for i = d, u, we
have Cd (γd ) ≥ Cd (γd ) = Cd (dw ) and Cu (γu ) ≥ Cu (γu ) =
Cu (xu ). Note that introducing the auxiliary variables γu,t with
constraints (17) and (18), and γd,t with constraints (19) and
(20) does not modify the problem. This means ufm ≥ uopt
m.
Hence, we have ufm = uopt
and
P
and
P1
are
equivalent.
m
f
From (44) and the objective of Pf , we have the T -slot driftplus-cost metric for the mth frame upper bounded by

H

∗
= Ci′−1 − V i,t
V Cu′ (γu,t ) = 0, respectively. We have γi,t
βi
for i = d, u.
∗
Thus, we have γi,t
for i = d, u as in (27).

(m+1)T −1

∆T (Θt ) + V

X

[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]

t=mT

14

(m+1)T −1

≤ Zt

X

∆u
Qt + Sr,t − Dt −
To

t=mT
t+T
X−1

+ Hu,t



X

+ Xt

Similarly, we have
max

(dτ − d

)

τ =t
t+T
X−1

(γu,τ − xu,τ ) + GT 2 + Hd,t

τ =t
(m+1)T −1

+V

t+T
X−1

(γd,τ − xd,τ )

τ =t

[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )] .

(45)

t=mT

Let {π̃t } denote a set of feasible solutions of Pf , satisfying
the following relations

(m+1)T −1 
 (m+1)T
X
X −1 
∆u
Q̃t + S̃r,t =
D̃t +
(46)
To
t=mT
(m+1)T −1

(m+1)T −1

X

γ̃i,t =

X

X

d˜t ≤

t=mT

t=mT
(m+1)T −1

x̃i,t , for i = u, d

(47)

t=mT
(m+1)T −1

t=mT

X

dmax

(48)

t=mT

with the corresponding objective value denoted as ũfm .
Note that comparing with P1, we impose per-frame constraints (46)-(48) as oppose to (16), (18), (20) and (3) for the
To -slot period, respectively. Let δ ≥ 0 denote the gap of ũfm
opt
f
to the optimal objective value uopt
m , i.e., ũm = um + δ.
Among all feasible control solutions satisfying (46)-(48),
there exists a solution which leads to δ → 0. The upper bound
in (45) can be rewritten as


(m+1)T −1
X
∆T (Θt ) + V 
[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]
t=mT


2
opt
≤ GT 2 + V T lim uopt
m + δ = GT + V T um .
δ→0

(49)

Summing both sides of (49) over m for m = 0, . . . , M − 1,
and dividing them by V M T , we have
M−1 (m+1)T
X −1
1 X
[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]
M T m=0
t=mT

+

M−1
GT
1 X opt
L(ΘTo ) − L(Θ0 )
≤
+
u .
V MT
V
M m=0 m

(50)

Since P
Ci (γi ) ≤ Ci (γi ) for the convex function Ci (·) where
To −1
γi,t for i = u, d, from (50), we have
γi , T1o t=0
!
To −1
1 X
Et Pt + xe + Cu (γu ) + αCd (γd )
To t=0
To −1
1 X
[Et Pt + xe,t + Cu (γu,t ) + αCd (γd,t )]
≤
To t=0

(51)

For a continuously differentiable convex function f (·), it
satisfies f (x) ≥ f (y) + f ′ (y)(x − y). Applying this to Cu (xu )
and Cu (γu ), we have
Cu (xu ) ≤ Cu (γu ) + Cu′ (xu )(xu − γu )
≤ Cu (γu ) + Cu′ (Γu )(xu − γu )
Hu,To − Hu,0
= Cu (γu ) − Cu′ (Γu )
To

(52)

Cd (dw ) ≤ Cd (γd ) − Cd′ (Γd )

Hd,To − Hd,0
.
To

(53)

Apply the inequalities in (52) and (53) to Cu (γu ) and Cd (γd )
respectively at LHS of (51), and combining (50) and (51), we
have the bound of the objective value u∗ (V ) of P1 in (30)
achieved by our proposed algorithm.
To show the bound in (31), as To → ∞, it is suffice to show
that both Hu,t and Hd,t in (30) are bounded. To show these
bounds, we need to show that the one-slot Lyapunov drift in
(35) is upper bounded as follows
L(Θt+1 ) − L(Θt ) ≤ G.

(54)

To show the above bound for the drift, we choose an alternative
feasible solution π̃ t satisfying the following per slot relations:
u
i) Q̃t + S̃r,t = D̃t + ∆
To ; ii) γ̃i,t = x̃i,t , for i = u, d; and
iii) d˜t ≤ dmax . With these relations, and by choosing d˜t =
dmax , the terms at RHS of (36), (38), (40) and (41) become
zeros, and we have (54). Averaging (54) over To -slot period,
we have T1o [L(ΘTo ) − L(Θ0 )] ≤ G. For any initial value of
L(Θ0 ) < +∞, by (26), we have

1  2
L(Θ0 )
2
2
. (55)
+ µ XT2o + Hd,T
Z + Hu,T
≤G+
o
o
2To To
To
p
2Tp
It follows that Hu,To ≤
o G + 2L(Θ0 ) and Hd,To ≤
p
(2To G + 2L(Θ0 ))/µ. Since 2To G + 2L(Θ0 )/To → 0 as
To → ∞, for any initial values of Hu,0 and Hd,0 , the third
term in RHS of (30) goes to zero. Thus, we have (31).

P ROOF

A PPENDIX F
OF P ROPOSITION 3

Proof: To prove ǫd is bounded, we note that
|XTo | + |X0 |
|XTo − X0 |
≤
.
(56)
T0
T0
p
From (55), it follows that |XTo | ≤ (2To G + 2L(Θ0 ))/µ.
Substituting the above upper bound of |XTo | in (56), we have
(32).
|ǫd | =

R EFERENCES
[1] T. Li and M. Dong, “Real-time joint energy storage management and
load scheduling with renewable integration,” in Proc. IEEE Int. Conf.
Acoust. Speech Signal Process. (ICASSP), Mar. 2016.
[2] A. Castillo and D. F. Gayme, “Grid-scale energy storage applications
in renewable energy integration: A survey,” Energy Convers. Manag.,
vol. 87, pp. 885–894, Nov. 2014.
[3] D. Callaway and I. Hiskens, “Achieving controllability of electric loads,”
Proc. IEEE, vol. 99, pp. 184–199, Jan. 2011.
[4] H.-I. Su and A. El Gamal, “Modeling and analysis of the role of energy
storage for renewable integration: Power balancing,” IEEE Trans. Power
Syst., vol. 28, no. 4, pp. 4109–4117, Nov. 2013.
[5] S. Sun, M. Dong, and B. Liang, “Real-time welfare-maximizing regulation allocation in dynamic aggregator-EVs system,” IEEE Trans. Smart
Grid, vol. 5, no. 3, pp. 1397–1409, May 2014.
[6] ——, “Real-time power balancing in electric grids with distributed
storage,” IEEE J. Sel. Topics Signal Process., vol. 8, no. 6, pp. 1167–
1181, Dec. 2014.
[7] Y. Wang, X. Lin, and M. Pedram, “Adaptive control for energy storage
systems in households with photovoltaic modules,” IEEE Trans. Smart
Grid, vol. 5, no. 2, pp. 992–1001, Mar. 2014.
[8] M. Lin, A. Wierman, L. L. H. Andrew, and E. Thereska, “Dynamic rightsizing for power-proportional data centers,” IEEE/ACM Trans. Netw.,
vol. 21, no. 5, pp. 1378–1391, Oct. 2013.

15

[9] Y. Zhang and M. van der Schaar, “Structure-aware stochastic storage
management in smart grids,” IEEE J. Sel. Topics Signal Process., vol. 8,
no. 6, pp. 1098–1110, Dec. 2014.
[10] R. Urgaonkar, B. Urgaonkar, M. J. Neely, and A. Sivasubramaniam,
“Optimal power cost management using stored energy in data centers,”
in Proc. ACM SIGMETRICS, Jun. 2011, pp. 115–120.
[11] L. Huang, J. Walrand, and K. Ramchandran, “Optimal demand response
with energy storage management,” in Proc. 3rd IEEE Int. Smart Grid
Commun. (SmartGridComm), Nov. 2012, pp. 61–66.
[12] S. Salinas, M. Li, P. Li, and Y. Fu, “Dynamic energy management for
the smart grid with distributed energy resources,” IEEE Trans. Smart
Grid, vol. 4, no. 4, pp. 2139–2151, Dec. 2013.
[13] T. Li and M. Dong, “Online control for energy storage management with
renewable energy integration,” in Proc. IEEE Int. Conf. Acoust. Speech
Signal Process. (ICASSP), May 2013, pp. 5248–5252.
[14] ——, “Real-time energy storage management with renewable integration: Finite-time horizon approach,” IEEE J. Select. Areas Commun.,
vol. 33, no. 12, pp. 2524–2539, Dec. 2015.
[15] M. J. Neely, Stochastic Network Optimization with Application to
Communication and Queueing Systems. USA: Morgan & Claypool,
2010.
[16] C. Chen, J. Wang, and S. Kishore, “A distributed direct load control
approach for large-scale residential demand response,” IEEE Trans.
Power Syst., vol. 29, no. 5, pp. 2219–2228, Sep. 2014.
[17] M. He, S. Murugesan, and J. Zhang, “A multi-timescale scheduling
approach for stochastic reliability in smart grids with wind generation
and opportunistic demand,” IEEE Trans. Smart Grid, vol. 4, no. 1, pp.
521–529, Mar. 2013.
[18] I. Koutsopoulos and L. Tassiulas, “Optimal control policies for power
demand scheduling in the smart grid,” IEEE J. Sel. Areas Commum.,
vol. 30, no. 6, pp. 1049–1060, Jul. 2012.
[19] C. Joe-Wong, S. Sen, S. Ha, and M. Chiang, “Optimized day-ahead
pricing for smart grids with device-specific scheduling flexibility,” IEEE
J. Sel. Areas Commum., vol. 30, no. 6, pp. 1075–1085, Jul. 2012.
[20] P. Samadi, A. H. Mohsenian-Rad, R. Schober, V. W. S. Wong, and
J. Jatskevich, “Optimal real-time pricing algorithm based on utility
maximization for smart grid,” in Proc. 1st IEEE Int. Conf. Smart Grid
Commun. (SmartGridComm), Oct. 2010, pp. 415–420.
[21] A. H. Mohsenian-Rad and A. Leon-Garcia, “Optimal residential load
control with price prediction in real-time electricity pricing environments,” IEEE Trans. Smart Grid, vol. 1, no. 2, pp. 120–133, Sep. 2010.
[22] P. Du and N. Lu, “Appliance commitment for household load scheduling,” IEEE Trans. Smart Grid, vol. 2, no. 2, pp. 411–419, Jun. 2011.
[23] S. Kishore and L. V. Snyder, “Control mechanisms for residential
electricity demand in smartgrids,” in Proc. 1st IEEE Int. Conf. Smart
Grid Commun. (SmartGridComm), Oct. 2010, pp. 443–448.
[24] T. T. Kim and H. V. Poor, “Scheduling power consumption with price
uncertainty,” IEEE Trans. Smart Grid, vol. 2, no. 3, pp. 519–527, Sep.
2011.
[25] T. H. Chang, M. Alizadeh, and A. Scaglione, “Real-time power balancing via decentralized coordinated home energy scheduling,” IEEE Trans.
Smart Grid, vol. 4, no. 3, pp. 1490–1504, Sep. 2013.
[26] P. Yi, X. Dong, A. Iwayemi, C. Zhou, and S. Li, “Real-time opportunistic
scheduling for residential demand response,” IEEE Trans. Smart Grid,
vol. 4, no. 1, pp. 227–234, Mar. 2013.
[27] A. H. Mohsenian-Rad, V. W. S. Wong, J. Jatskevich, R. Schober,
and A. Leon-Garcia, “Autonomous demand-side management based
on game-theoretic energy consumption scheduling for the future smart
grid,” IEEE Trans. Smart Grid, vol. 1, no. 3, pp. 320–331, Dec. 2010.
[28] B. G. Kim, S. Ren, M. van der Schaar, and J. W. Lee, “Bidirectional
energy trading and residential load scheduling with electric vehicles in
the smart grid,” IEEE J. Select. Areas Commun., vol. 31, no. 7, pp.
1219–1234, Jul. 2013.
[29] S. Chen, N. B. Shroff, and P. Sinha, “Heterogeneous delay tolerant task
scheduling and energy management in the smart grid with renewable
energy,” IEEE J. Sel. Areas Commum., vol. 31, no. 7, pp. 1258–1267,
Jul. 2013.
[30] Y. Huang, S. Mao, and R. M. Nelms, “Adaptive electricity scheduling
in microgrids,” IEEE Trans. Smart Grid, vol. 5, no. 1, pp. 270–281, Jan.
2014.
[31] S. Sun, M. Dong, and B. Liang, “Distributed real-time power balancing
in renewable-integrated power grids with storage and flexible loads,”
IEEE Trans. Smart Grid, vol. PP, no. 99, pp. 1–1, 2015.
[32] P. Ramadass, B. Haran, R. White, and B. N. Popov, “Performance study
of commercial licoo2 and spinel-based li-ion cells,” J. Power Sources,
vol. 111, pp. 210–220, Sep. 2002.
[33] Z. Ma, S. Zou, and X. Liu, “A distributed charging coordination for
large-scale plug-in electric vehicles considering battery degradation

[34]

[35]
[36]
[37]

[38]
[39]
[40]

cost,” IEEE Trans. Control Syst. Technol., vol. 23, no. 5, pp. 2044–2052,
Sep. 2015.
N. Michelusi, L. Badia, R. Carli, L. Corradini, and M. Zorzi, “Energy
management policies for harvesting-based wireless sensor devices with
battery degradation,” IEEE Trans. Commun., vol. 61, no. 12, pp. 4934–
4947, Dec. 2013.
I. Duggal and B. Venkatesh, “Short-term scheduling of thermal generators and battery storage with depth of discharge-based cost model,”
IEEE Trans. Power Syst., vol. 30, no. 4, pp. 2110–2118, Jul. 2015.
Y. Guo, M. Pan, Y. Fang, and P. P. Khargonekar, “Decentralized
coordination of energy utilization for residential households in the smart
grid,” IEEE Trans. Smart Grid, vol. 4, no. 3, pp. 1341–1350, Sep. 2013.
Z. Ma, S. Zou, and X. Liu, “A distributed charging coordination for
large-scale plug-in electric vehicles considering battery degradation
cost,” IEEE Trans. Control Syst. Technol., vol. 23, no. 5, pp. 2044–2052,
Sep. 2015.
M. J. Neely, “Universal scheduling for networks with arbitrary traffic,
channels, and mobility,” Tech. Rep. arXiv: 1001.0960v1, Jan. 2010.
T. Li and M. Dong, “Real-time energy storage management with
renewable energy of arbitrary generation dynamics,” in Proc. of Asilomar
Conf. on Signals, Systems and Computers, Nov. 2013, pp. 1714–1718.
Ontario Energy Board (2012). Electricity Price. [Online]. Available:
http://www.ontarioenergyboard.ca/OEB/Consumers/Electricity

Tianyi Li (S’11) received the B.S degrees in communication engineering from North China University of Technology, Beijing, China, and in telecommunication engineering technology from Southern
Polytechnic State University, Marietta, GA, USA
both in 2009, and the M.Sc degree in electrical engineering from Northwestern University, Evanston,
IL, USA in 2011. He obtained his Ph.D degree in
2015 from the Department of Electrical, Computer
and Software Engineering, University of Ontario
Institute of Technology, Oshawa, Ontario, Canada,
and currently is a research associate in the same university. He received
the best student paper in the Technical Track of Signal Processing for
Communications and Networking (SP-COM) at IEEE ICASSP in 2016. His
research interests include stochastic optimization and control in energy storage
management and demand side management of smart grid.

Min Dong (S’00-M’05-SM’09) received the B.Eng.
degree from Tsinghua University, Beijing, China, in
1998, and the Ph.D. degree in electrical and computer engineering with minor in applied mathematics
from Cornell University, Ithaca, NY, in 2004. From
2004 to 2008, she was with Corporate Research
and Development, Qualcomm Inc., San Diego, CA.
In 2008, she joined the Department of Electrical,
Computer and Software Engineering at University
of Ontario Institute of Technology, Ontario, Canada,
where she is currently an Associate Professor. She
also holds a status-only Associate Professor appointment with the Department
of Electrical and Computer Engineering at University of Toronto. Her research
interests are in the areas of statistical signal processing for communication networks, cooperative communications and networking techniques, and stochastic
network optimization in dynamic networks and systems.
Dr. Dong received the Early Researcher Award from Ontario Ministry of
Research and Innovation in 2012, the Best Paper Award at IEEE ICCC in
2012, and the 2004 IEEE Signal Processing Society Best Paper Award. She
is a co-author of the best student paper in the Technical Track of Signal
Processing for Communications and Networking (SP-COM) at IEEE ICASSP
in 2016. She served as an Associate Editor for the IEEE TRANSACTIONS
ON SIGNAL PROCESSING (2010-2014), and as an Associate Editor for the
IEEE SIGNAL PROCESSING LETTERS (2009-2013). She was a symposium
lead co-chair of the Communications and Networks to Enable the Smart Grid
Symposium at the IEEE International Conference on Smart Grid Communications (SmartGridComm) in 2014. She has been an elected member of
IEEE Signal Processing Society Signal Processing for Communications and
Networking (SP-COM) Technical Committee since 2013.

