1

A Decentralized Framework for Real-Time Energy Trading in
Distribution Networks with Load and Generation Uncertainty

arXiv:1705.02575v1 [] 7 May 2017

Shahab Bahrami and M. Hadi Amini
email: bahramis@ece.ubc.ca, amini@cmu.edu
Abstract—The proliferation of small-scale renewable generators and price-responsive loads makes it a challenge for distribution network operators (DNOs) to schedule the controllable loads
of the load aggregators and the generation of the generators
in real-time. Additionally, the high computational burden and
violation of the entities’ (i.e., load aggregators’ and generators’)
privacy make a centralized framework impractical. In this paper,
we propose a decentralized energy trading algorithm that can be
executed by the entities in a real-time fashion. To address the
privacy issues, the DNO provides the entities with proper control
signals using the Lagrange relaxation technique to motivate them
towards an operating point with maximum profit for entities.
To deal with uncertainty issues, we propose a probabilistic load
model and robust framework for renewable generation. The
performance of the proposed algorithm is evaluated on an IEEE
123-node test feeder. When compared with a benchmark of not
performing load management for the aggregators, the proposed
algorithm benefits both the load aggregators and generators
by increasing their profit by 17.8% and 10.3%, respectively.
When compared with a centralized approach, our algorithm
converges to the solution of the DNO’s centralized problem with
a significantly lower running time in 50 iterations per time slot.
Keywords: price-responsive load, generation uncertainty, distributed algorithm, trading market.

I. I NTRODUCTION
One goal of the emerging smart grid is to move distribution
systems towards a smarter and more secure network through
integrating two-way communication infrastructure. The information exchange provides distribution network operators
(DNOs) with sophisticated management and monitoring systems to perform complex analyses and automated operations
in near real-time. Furthermore, drivers such as distribution
organizations have accelerated the expansion of applications
for smart grid technologies, such as smart meters, and integration of renewable energy generators. Resulting benefits include
a more efficient use of electric appliances in households to
reduce the energy bill payment for the load aggregators and
lower operation cost for the generators, as well as a higher
flexibility for the DNO to enhance the system’s technical
operation; thereby reaching a triple-win condition.
The DNO is responsible for optimal power flow (OPF)
analysis. There are challenges in solving the OPF problem
by the DNO. First, the OPF can be computationally difficult
to be solved, especially when the number of decision variables
increases by participation of the price-responsive load aggregators in the energy market. Second, the DNO may violate
the entities’ privacy, e.g., by revealing the load aggregators’
demand information and generators’ cost parameters to the
DNO. Third, the DNO is uncertain about the load demand
and renewable generation ahead of time.
There have been some efforts in the literature to tackle
the above-mentioned challenges. We divide the related works
into three main threads. The first thread is concerned with
decentralized energy management programs for a market with
multiple suppliers and multiple users. Mechnisms such as the

multi-level game methods [1], Stackelberg game [2], dual
decomposition method [3], supply bidding framework [4],
and hierarchical bidding [5] have been used. However, these
approaches did not consider the constraints imposed by the
topology and operation of the distribution network. The second
thread is concerned with including the power flow equations in
the decentralized energy management procedure. To achieve
this goal, different techniques such as convex relaxation [6]–
[9], quadratic programming [10], alternating direction method
of multipliers (ADMM) [11]–[13], and Lagrange relaxation
method [14], [15] have been used. These studies, however,
did not consider the uncertainties in the renewable generators
and load demand. Furthermore, these studies mainly focused
on off-line algorithms, which are applicable in day-ahead
markets. The third thread concerned with the online operation
of distribution systems using different mechanisms such as
real-time closed-loop control [16], differential evolution optimization [17], online gradient method [18], projected gradient
descent [19], online mirror descent [20], and graph theorybased approach [21]. These works, however, did not mentioned
how to consider the uncertainty in the load demand for users
in smart distribution networks.
In this paper, we focus on designing a distributed algorithm
for an electricity trading market with renewable energy generators and price-responsive residential load aggregators. In each
time slot (e.g., every hour), the load aggregators and generators
use the communication infrastructure in the smart grid to
exchange information with the DNO and jointly maximize
their profit, while considering the uncertainty in the future
demand and renewable generation. The privacy of each entity
is protected in the proposed framework, as the generators and
load aggregators solve their own profit maximization problem
using the locally available information. The main challenges
that we address are tackling the uncertainty in the load demand
and renewable generation, as well as determining the proper
control signals communicated between the DNO, generators,
and load aggregators that enforce the proposed distributed
algorithm to converge to the solution of the DNO’s centralized
problem with the objective of maximizing the social welfare.
The main contributions of this paper are as follows:
• Uncertainty Issues and Risk Evaluation: To address the
uncertainty in the load demand, we propose a probabilistic load estimation for the electric appliances of
the residential users served by each load aggregator. It
enables each load aggregator to schedule the electric
appliances of its users in real-time, while taking into
account the impacts of its decision on the load profile
in the upcoming time slots. We also consider an adaptive
robust decision making framework for the renewable generators to optimize the risk of power shortage based on an
adjustable confidence level. It enables the generators to
limit their cost for compensating the generation shortage.

2

It also enables the DNO to prevent high voltage drop
caused by the shortage in the total renewable generation.
• Distributed Algorithm Design: To protect the privacy of
the load aggregators and generators, as well as to address
the computational complexity of the DNO’s centralized
problem, we propose a decentralized algorithm that can
be executed by the entities in real-time. Each entity
requires to share limited information to meet its local
objective, while satisfying the physical constraints of the
linearized ac power flow in the distribution network.
• Performance Evaluation: Simulations are performed on
an IEEE 123-bus test feeder with 10 generators and 113
load aggregators. When compared with the benchmark of
not performing load management, the proposed algorithm
benefits the load aggregators and generators by increasing
their profit by 17.8% and 10.3% on average, respectively.
Furthermore, it helps generators to reduce the peak-toaverage generation ratio by 13%. Our algorithm converges to the solution of the centralized problem with
a significantly lower execution time.
The rest of this paper is organized as follows. Section II
introduces the system model. Section III formulates the DNO’s
centralized and decentralized problems. A decentralized algorithm is proposed. Section IV provides the simulation results,
followed by Section V that concludes the paper. Appendices
A−F can be found in the supplementary document.
II. S YSTEM M ODEL
Consider an electricity market with a set N of N , |N |
load aggregators and a set G of G , |G| generators scattered
in a distribution network. Each load aggregator is responsible
for managing the load demand of its electricity users. Each
generator sells electricity to the market. The load aggregators
and generators use a two-way communication infrastructure
to exchange information with the DNO. The entities are also
connected to each other through the electric power distribution
network. The DNO is a neutral entity responsible for monitoring the power flow in the network. For simplicity in the
problem formulation, we assume that each bus in the network
has exactly one load aggregator or one generator. If both load
aggregator and generator are connected to the same bus, we
divide that bus into two buses connected to each other through
a line with zero impedance. If neither load aggregator nor
generator is connected to a bus, we propose to add a virtual
load aggeragtor with zero demand for that bus. It enables
us to denote the set of buses by N ∪ G and refer a load
aggregator or a generator by its bus index. We use notation
L ⊆ (N ∪ G) × (N ∪ G) to denote the set of branches.
The overall trading horizon is denoted by H , {1, . . . , H},
where H is the number of time slots with equal length (e.g.,
each time slot is one hour). Notice that the load management
decision of a load aggregator in the current time slot affects its
demand in the upcoming time slots. Meanwhile, the generators
need to match their generation level with the changes in the
load demand. Hence, generators also need to modify their
generation for the current and upcoming time slots. To avoid
an abuse of notations, hereafter, we use index h for a time slot
in general and use index t specifically for the current time slot.

The general idea of this paper for implementing a real-time
energy trading can be summarized as follows. At the beginning
of the current time slot t, the entities optimize the demand and
generation profiles over the period Ht = {t, . . . , H} ⊆ H, but
apply only the obtained decision for the current time slot t.
The scheduling procedure is performed with uncertainty about
the load demand and renewable generation in the upcoming
time slots h > t. Hence, the entities repeat the optimization
procedure at the beginning of the next time slot to update
their scheduling decision with the revealed demand/generation
information. We aim to answer two key questions:
Q.1 How do the entities interact with the DNO to determine
their optimal load an generation profiles in the current time
slot with the locally available information?
Q.2 How do the entities address the lack of information about
the demand and generation in the upcoming time slots?
A. Load Aggregator’s Model
In this subsection, we address questions Q.1 and Q.2 for residential load aggregators by modeling the electric appliances
and providing a probabilistic load estimation technique.
1) Users’ Appliances Model: Load aggregator i ∈ N is
responsible for scheduling its users’ electric appliances. An
electric appliance is either asleep or awake in the current time
slot t. Let Aasleep
(t) and Aawake
(t) denote the sets of asleep and
i
i
awake appliances in the current time slot t, respectively. An
awake appliance a ∈ Aawake
(t) is available to be scheduled
i
for operation, i.e., the load aggregator schedules the power
consumption profile ea,i (t) = (ea,i (h), h ∈ Ht ).
The awake appliance a ∈ Aawake
(t) provides the load aggrei
gator i with its scheduling horizon, utility function, and type
using the smart meter inside the household. The scheduling
horizon Ha,i ⊆ Ht defines the time interval over the upcoming
time slots, in which the appliance should be scheduled. The
utility function Ua,i (ea,i (t)) is used to model the satisfaction
of the customer in monetary units from using the appliance.
It is generally an increasing and concave function of the
consumption profile ea,i (t). The type of appliance depends
on its specifications and the customer’s preferences. Inspired
by the work in [22], we consider three types of appliances:
• The appliance a with type 1 should be operated within
the scheduling horizon Ha,i and turned off in other time slots.
Examples include the electric vehicle (EV) and dish washer.
(t) denote the set of appliances with type
Let A1i (t) ⊆ Aawake
i
1 that are awake in the current time slot t. We have
ea,i (h) = 0,

a ∈ A1i (t), i ∈ N , h 6∈ Ha,i , (1a)

max
1
emin
a,i (h) ≤ ea,i (h) ≤ ea,i (h), a ∈ Ai (t), i ∈ N , h ∈ Ha,i , (1b)
P
min
max
, a ∈ A1i (t), i ∈ N . (1c)
Ea,i
≤ h∈Ha,i ea,i (h) ≤ Ea,i

The utility obtained from using a type 1 appliance depends on the total power consumption.
P The utility can be
e (h) , e.g.,
expressed as Ua,i (ea,i (t)) = Ua,i


P
P h∈Ha,i a,i
min
Ua,i
e
(h)
=
κ
f
e
(h) − Ea,i
a,i
a,i
a,i
h∈Ha,i
h∈Ha,i
with a concave function f (·) and nonnegative constant κa,i .
• The appliances of type 2 can be operated in time slots
out of the scheduling horizon, but the customer attains a
relatively low utility, e.g., TV and personal computer. Let

3

A2i (t) ⊆ Aawake
(t) denote the set of appliances of type 2 that
i
are awake in the current time slot t. We have
a ∈ A2i (t), i ∈ N , h 6∈ Ha,i , (2a)

ea,i (h) ≥ 0,

max
2
emin
a,i (h) ≤ ea,i (h) ≤ ea,i (h), a ∈ Ai (t), i ∈ N , h ∈ Ha,i , (2b)
P
min
max
, a ∈ A2i (t), i ∈ N . (2c)
Ea,i
≤ h∈Ha,i ea,i (h) ≤ Ea,i

The utility function for type 2 appliances depends on both
the amount of power consumption and the time of consuming
the power, i.e., the customer would gain different benefits
from consuming the same amount of power at different
times, e.g., watching
the favorite TV program. We have
P
Ua,i (ea,i (t)) = h∈Ht Ua,i (ea,i (h),
P h). As a concrete example, utility function Ua,i (ea,i ) = h∈Ha,i κa,i (h)f ea,i (h) −

 P
′
emin
a,i +
k6∈Ha,i κa,i (h)f ea,i (h) with a concave function
f (·) and time dependent nonnegative coefficients κa,i (h) and
κ′a,i (h), κ′a,i (h) ≪ κa,i (h) is a viable candidate.
• The appliances of type 3 can be operated out of the
scheduling horizon without any constraint on their total power
consumption, such as lighting and refrigerator. Let A3i (t) ⊆
Aawake
(t) denote the set of appliances of type 3 that are awake
i
in the current time slot t. We have
a ∈ A3i (t), i ∈ N , h 6∈ Ha,i , (3a)

ea,i (h) ≥ 0,

max
3
emin
a,i (h) ≤ ea,i (h) ≤ ea,i (h), a ∈ Ai (t), i ∈ N , h ∈ Ha,i . (3b)

The utility Ua,i (ea,i ) attained by the customer from using
the appliances with type 3 depends on the amount of power
consumption ea,i (t) within the scheduling horizon Ha,i , but
not the time of consumption. The customer attains a relatively low utility out of interval
Ha,i . Function Ua,i (ea,i
 )=
 P
P
′
min
h6∈Ha,i κa,i f ea,i (h) with
h∈Ha,i κa,i f ea,i (h)−ea,i +
a concave function f (·) and nonnegative constants κa,i and
κ′a,i , κ′a,i ≪ κa,i is a viable candidate.
The total utility of load aggregator i in the current time slot
t with decision vector ei (t) = (ea,i (t), a ∈ Aawake
(t)) is
i
P
(4)
Ui (ei (t)) = a∈Aawake (t) Ua,i (ea,i (t)), i ∈ N .
i

2) Load Estimation: The actual wake-up times of the
appliances are not available to the load aggregator in advance.
To address this lack of information, load aggregator i can
collect the sleep-awake historical data record of each appliance
and estimate the probability pa,i (h) that each appliance a ∈ Ai
becomes awake at each time slot h ∈ H. In appendix A, we
show the conditional probability pa,i (h | t) that the appliance
becomes awake in an upcoming time slot h > t, given that it
has not become awake until the current time slot, t, is
pa,i (h | t) =

1−

pa,i (h)
.
Pt
′
h′ =1 pa,i (h )

(5)

A load aggregator has no information about the scheduling
horizon, users’ utility, and type of the appliances ahead of
time. For decision making at the current time slot t, we consider the worst-case scenario, in which the electric appliances
that become awake in the upcoming time slots h > t should
be operated once they become awake without any control
max
nom
on power consumption, i.e., emin
a,i (h) = ea,i (h) = ea,i and
min
max
nom
Ea,i = Ea,i = Ea,i . The payment of the load aggregator

in the worst-case scenario is an upper-bound for its actual
payment. Hence, minimizing the worst-case payment implies
reducing the risk of high payment. For the current time slot
t, the worst-case expected electric demand liasleep(h) of the
currently sleeping appliances in an upcoming time slot h > t is
#
"
h
X
X
asleep
′
nom
pa,i (h | t) . (6)
li
(h) =
ea,i (h)
a∈Aasleep
(t)
i

h′ =max{t+1,h−Ta +1}

nom nom
where parameter Ta = Ea,i
/ea,i is the operation duration of
asleep
the appliance a ∈ Ai
(t) that becomes awake in upcoming
P
time slot h > t. The value of kh′ =max{t+1,h−Ta +1} pa (h′ | t)
is equal to the probability that a currently sleeping appliance
a is operating in the upcoming time slot h > t.
For the given current time slot t, we use the notation
li (t) = (li (h), h ∈ Ht ) to denote the profile of active power
consumption of the users during time interval Ht . We have
P
h ∈ Ht .
(7)
li (h) =liasleep (h) + a∈Aawake (t) ea,i (h),
i

To model the reactive power consumption for a load aggregator i, we consider a constant power factor PFi . Theqreactive
1−PF2i
.
power for load aggregator i in time h is qi (h) = li (h)
PF2i
3) Local Optimization Problem: Constraints (1a)−(7) define the feasible set Ei (t) for decision vector ei (t) of load
aggregator i in the current time slot t. Load aggregator i aims
to maximize the profit πiagg (ei (t)), which includes the total
utility in (4) minus the payment to the DNO over period Ht .
The DNO provides load aggregator i with the price ρi (h) for
a unit of active power in each time slot h. We assume that
load aggregators do not pay for the reactive power. We have
P
πiagg (ei (t)) = Ui (ei (t)) − h∈Ht li (h) ρi (h), i ∈ N , (8)
Load aggregator i solves the following optimization problem in time slot t to determine decision vector ei (t):
maximize πiagg (ei (t))

(9a)

subject to ei (t) ∈ Ei (t).

(9b)

ei (t)

B. Generator’s Model
In this subsection, we address questions Q.1 and Q.2 for the
generators by modeling the conventional and renewable units
and providing a robust optimization technique for renewables.
1) Conventional Unit: In general, the generation cost function of the conventional unit of generator j ∈ G in time
slot h is an increasing convex function of pconv
(h) [23]. The
j
2
conv
class of quadratic functions Cj (pj (h)) = aj2 pconv
j (h)
+aj1 pconv
(h) + aj0 is well-known. For the given current
j
time slot t, generator j offers the profiles of active and
conv
reactive powers pconv
(h), h ∈ Ht ) and qjconv (t) =
j (t) = (pj
conv
(qj (h), h ∈ Ht ) for the current and upcoming time slots.
2) Renewable Unit: Without loss of generality, we assume
that the renewable plants are operated at unity power factor.
Given the current time slot t, generator j with renewable unit
ren
offers an active power profile pren
j (t) = (pj (h), h ∈ Ht ) for
its renewable unit. To prevent non-credible high renewable
generation offers, the DNO charges generator j by the unit
price βj (h) ($/MW) for generation shortage in time slot h. To

4

cope with the uncertainty issues, we consider a robust decision
making for generators with renewable units. Generator j can
uses the historical data record to forecast an uncertainty bound
[pmin,ren
(h), pmax,ren
(h)] for its actual renewable generation in
j
j
time slot h ∈ H. Generator j considers the cost Γj (pren
j (t))
of the worst-case scenario for generation shortage as follows:
P
min,ren
ren
(h)). (10)
Γj (pren
j (t)) ,
h∈Ht βj (h)(pj (h) − pj

The feasible set for the renewable generation profile pren
j (t)
can be defined based on all scenarios that satisfy pren
(h)
∈
j
max,ren
(h),
p
(h)].
However,
it
is
very
conservative
and
[pmin,ren
j
j
possibly inefficient to take into account all possible scenarios.
Inspired by the work in [24], we consider an adaptive robust
model. In the current time slot t, the uncertainty space for the
generation profile pren
j (t) in the time interval Ht is defined as
n
min,ren
ren
(h), pmax,ren
(h)], h ∈ Ht ,
Pjren (t) = pren
j (t) | pj (h) ∈ [pj
j
max,ren
ren
o
pj
(h) − pj (h)
P
≤
∆
(t)
, (11)
j
h∈Ht max,ren
pj
(h) − pmin,ren
(h)
j

where 0 ≤ ∆j (t) ≤ |Ht | is the confidence level parameter
for generator j in the current time slot t. The space defined
in (11) is a singleton, corresponding to the least-conservative
max,ren
(h), h ∈ Ht when ∆j (t) = 0. As
scenario pren
j (h) = pj
∆j (t) increases, the size of the uncertainty set enlarges, and
the resulting robust solution is more conservative. The space
includes all possible scenarios when ∆j (t) = |Ht |. In [24],
∆j (t) is known and fixed. Whereas, we consider parameter
∆j (t) as a variable that should be optimized by generator j.
3) Local Optimization Problem: For a given current time
slot t, generator j decides on the generation profile ψj (t) =
(pconv
(t), qjconv (t), pren
j
j (t)) and the confidence level ∆j (t). The
gen
objective of generator j is to maximize its profit πj (ψj (t)),
which is the revenue from selling active and reactive powers
minus the generation cost and financial risk in (10). That is


P
gen
(h)+pren
πj (ψj (t)) = h∈Ht pconv
j
j (h) ρj (h)

(h)) − Γj (pren
(12)
+ qjconv (h)̺j (h) − Cj (pconv
j
j (t)).

The problem for generator j ∈ G in the current time slot t is
maximize πjgen (ψj (t))

ψj (t), ∆j (t)

(13a)

subject to pmin,conv
≤ pconv
(h) ≤ pmax,conv
, h ∈ Ht , (13b)
j
j
j
qjmin,conv ≤ qjconv (h) ≤ qjmax,conv , h ∈ Ht ,

(13c)

pren
j (t)

(13d)

∈

Pjren (t).

C. DNO’s Model
We address Q.1 and Q.2 for the DNO in the following.
1) Linearized ac Power flow: We consider balanced distribution networks. The model for unbalanced networks is
a research for future work. The ac power flow equations
are nonlinear and nonconvex. An alternative is to use a
linearized ac power flow. Let p(h) = (pb (h), b ∈ N ∪ G) and
q(h) = (qb (h), b ∈ N ∪ G) denote the vectors of injected
active power pb (h) and reactive power qb (h) to all buses
b ∈ N ∪ G in time h. Let |vb (h)| and θb (h) denote the
voltage magnitude and phase angle of bus b in time h. We
define the grid-wide vectors θ(h) = (θb (h), b ∈ N ∪ G) and

|v(h)| = (|vb (h)|, b ∈ N ∪ G) in time slot h. Let Grs and
Brs denote the real and reactive parts of the entry (r, s) in
bus admittance matrix Y . Let brr and grr denote the shunt
susceptance and conductance at bus r. In Appendix B, we
show that the linearized ac power flow can be expressed as


 

θ(h)
p(h)
−B ′
G′
,
(14)
=
q(h)
−G′ −B |v(h)|
{z
}
|
Λ

where the diagonal element (r, r) of matrices B and B ′ are
Brr and Brr − brr , respectively. The non-diagonal elements
(r, s) of both B and B ′ are Brs . The diagonal element (r, r)
of matrices G and G′ are Grr and Grr − grr , respectively.
The non-diagonal elements (r, s) of both G and G′ are Grs .
In Appendix C, we show that, in time slot h, the linearized
active and reactive power flow through line (r, s) ∈ L with
resistance Rrs and reactance Xrs can be calculated as:
Rrs (|vr (h)|−|vs (h)|)+Xrs (θr (h)−θs (h))
, (15a)
2 + X2
Rrs
rs
Xrs (|vr (h)|−|vs (h)|)−Rrs (θr (h)−θs (h))
qrs (h) =
, (15b)
2 + X2
Rrs
rs
p
2 (h) is
The apparent power flow srs (h) = p2rs (h) + qrs
upper bounded by smax
,
which
implies
that
the
feasible
real
rs
and reactive powers are bounded by a circle. To linearize
the constraint, a piecewise approximation of the boundary
by a regular polygon with central angle α can be used. In
Appendix D, we obtain the following constraints
prs (h) =

prs (h) cos (mα) + qrs (h) sin (mα) ≤ smax
rs ,

(16)

where m = 0, . . . , 2π/α. For each bus b, we also have
vbmin ≤ |vb (h)| ≤ vbmax .
(17)
2) DNO’s Centralized Optimization Problem: The DNO
considers the impact of renewable generation shortage on the
technical operation of the network. As a concrete example, the
risk of voltage drop at different buses is a viable choice for the
DNO. The DNO considers function ΓDNO (t) of the grid-wide
renewable generation pren (t) for the voltage variations as
P
P
vb (h)|) , (18)
ΓDNO (pren (t)) , h∈Ht
b∈N ∪G (|vb (h)| − |b
where vbb (h) is voltage magnitude of bus b in the worst-case
scenario, when all renewable generators’ power are pmin,ren
(h)
j
in time slot h ∈ Ht . The uncertainty space that the DNO
considers for renewable generator j is defined as (11).
We consider the objective of maximizing the social welfare
for the DNO. Considering the grid-wide vectors e(t) and
pconv (t), and pren (t), the DNO’s objective function is
P
f DNO (e(t), pconv (t), pren (t)) , i∈N Ui (ei (t))
P
P
c DNO ren
(p (t)), (19)
− h∈Ht j∈G Cj (pconv
j (h)) −ϑ Γ

where ϑc is a positive weighting coefficient. We formulate the
DNO’s centralized problem as
maximize

e(t),pconv (t),q conv (t),
pren (t),∆(t),|v(t)|,θ(t)

f DNO (e(t), pconv (t), pren (t))

subject to constraints (9b), (13b)−(13d), (14)−(18).

(20a)
(20b)

5

Problem (20) has a concave objective function and linear
constraints due to the concavity of the load aggregators’ utility
function, the convexity of the generation cost function, and the
linearity of the ac power flow model in (14). Hence, we have:
Theorem 1 The optimal solution to the DNO’s centralized
problem in (20) is unique.
To solve the centralized problem (20), the DNO needs the
information about the load aggregators’ utilities, generators’
generation cost, and renewable units’ forecast data. However,
these information may not be available to the DNO. Instead,
we develop a decentralized algorithm by showing that the
DNO can determine ρb (h) and ̺b (h), b ∈ N ∪ M, as well as
the penalties βj (h), j ∈ G for h ∈ Ht such that when the load
aggregators solve (9) and generators solve (13), the resulting
solution coincides with the unique solution of problem (20).
III. D ECENTRALIZED A LGORITHM D ESIGN
Given the current time slot t, the decision vector of load
aggregator i is load profile ei (t) of the awake appliances.
Further, the decision vector of generator j is the generation
profile ψj (t) = (pconv
(t), qjconv (t), pren
j
j (t)) and the confidence
level ∆j (t). The DNO influences the entities by using the
nodal prices ρ(t), ̺(t), and penalties β(t).
One well-know technique to determine the proper values
of the nodal prices ρ(t), ̺(t), and penalties β(t) is to
formulate the partial Lagrangian relaxation of the DNO’s
centralized problem (20) [14], [15], [25]. Let λb (h) and
γb (h), b ∈ N ∪ G, h ∈ H denote the Lagrange multipliers
associated with the equality constraints for the injected active
power pb (h) and reactive power qb (h) in (14). We move these
constraints with their Lagrange multipliers to the objective
function of the centralized problem in (20). In equation (SDNO
12) of Appendix E, we obtain the objective function fLag
(·)
of the relaxed problem. Due to the convexity of problem
(20) and linearity of the constraints, the strong duality gap
condition (Slater’s condition) is satisfied if a feasible solution
exist [26, Ch. 5]. Thus, the optimal solution to the relaxed
problem is equal to the optimal solution to the primal problem
(20). Using the relaxed problem enables us to determine the
price signals ρ(t), ̺(t), and β(t) based on the Lagrangian
decomposition technique, such that the market equilibrium
among load aggregators and generators coincides with the
optimal solution of the centralized problem (20).
Theorem 2 The equilibrium of the energy market coincides
with the unique solution to the DNO’s centralized problem in
(20) if and only if for i ∈ N , j ∈ G, h ∈ Ht the DNO sets
p
ρi (h) = λi (h) + γi (h) (1 − PF2i )/PF2i , i ∈ N , h ∈ Ht , (21a)
ρj (h) = λj (h),
j ∈ G, h ∈ Ht , (21b)
̺j (h) = γj (h),
j ∈ G, h ∈ Ht , (21c)
P
(21d)
βj (h) = ϑc b∈N ∪G Λ−1 (|N ∪ G| + b, j),

where Λ−1 (|N ∪ G| + b, j) is the entry (|N ∪ G| + b, j) of the
inverse of matrix Λ in (14).
The proof can be found in Appendix E. It suggests a decentralized algorithm to determine the solution to problem (20).
We propose Algorithm 1 that can be executed by the load
aggregators, generators, and DNO in real-time. In Algorithm

Algorithm 1 Decentralized Energy Market Trading Algorithm.
1: Set k := 1 and ξ1 = ξ2 := 10−2 .
2: If t = 1
3: Each load aggregator i ∈ N randomly initializes its users’ appliances
load profile e1i (t).
4: Each generator j ∈ G randomly initializes its conventional generation
profiles pjconv,1 (t) and qjconv,1 (t).
5: Each generator j with renewable units initializes ∆1j (t) = 0 and set the
presumed generation levels to pren,1
(h) = pmax,ren
(h) for h ∈ Ht .
j
j
1
1
6: The DNO sets |vb (h)| = 1 pu, θb = 0, and λ1b (h) = γb1 (h) = 0, b ∈
N ∪ G, h ∈ Ht .
7: Else if t > 1
8: Load aggregators, generators, and DNO initialize their decision variables
with their values in the equilibrium at previous time slot t − 1.
9: End if
10: Repeat
11: Each load aggregator i and generator j sends its load profile lik (t) and
generation profiles pconv,k
(t), qjconv,k (t), and pren,k
(t) to the DNO.
j
j
12: DNO obtains the updated vector φk+1 (t) = (θbk+1 (t), |vbk+1 (t)|,
λk+1
(t), γbk+1 (t), b ∈ N ∪ G) using (22).
b
13: DNO uses (21a)−(21d) to compute the updated values of control signals
ρk+1 (t), ̺k+1 (t), and β k+1 (t), and sends the control signals to the
corresponding entity in each bus.
14: Each load aggregator i updates its controllable load profile ek+1
(t) by
i
solving its local problem (9)
15: Each generator j updates its generation profile ψjk (t) and decision
variable ∆kj (t), by solving its local problem in (13).
16: k := k + 1. The step size is updated.
17: Until |θbk (t) − θbk−1 (t)| ≤ ξ1 , ||vbk (t)| − |vbk−1 (t)|| ≤ ξ2 , b ∈ N ∪ G.

1, when the current time slot t begins, each load aggregator i determines the power consumption profile ea,i (t) =
(ea,i (h), h ∈ Ht ) of all awake appliances a ∈ Aawake
(t) over
i
time slots h ∈ Ht . Each generator j obtains the profiles of
active and reactive powers pconv
(t) = (pconv
(h), h ∈ Ht ) and
j
j
conv
conv
qj (t) = (qj (h), h ∈ Ht ) of the conventional unit and
ren
generation profile pren
j (t) = (pj (h), h ∈ Ht ). The entities
use the obtained scheduling decision for upcoming time slots
h ≥ t + 1 as an initial decision in the next time slot t + 1.
In each time slot t, Algorithm 1 is executed in an iterative
fashion. Let k denote the iteration index. Our algorithm
involves the initiation phase and market trading phase.
• Initiation phase: Lines 1 to 9 describe the initiation phase.
• Market trading phase: The loop involving Lines 10 to 18
describes this phase, which includes the following parts:
a) Information exchange: In Line 11, each load aggregator i
uses (7) to obtain its demand profile lik (t) = (lik (h), h ∈ Ht ),
and sends to the DNO. Each generator j sends the profiles
(t), qjconv,k (t), and pren,k
(t) to the DNO.
pconv,k
j
j
b) DNO’s update: In Line 12, the DNO receives the
information from the entities, it obtains the updated vector
φk+1(t) = (θbk+1 (t), |vbk+1 (t)|, λk+1
(t),γbk+1 (t), b ∈ N∪G) as
b

+
DNO
φk+1 (t) = φk (t) + ǫk ∇φk (t) fLag
(·) ,
(22)

where ∇ is the gradient operator, and [·]+ is the projection
onto the feasible set defined by constraints in (20b). Recall
DNO
that fLag
(·) is the objective function of the DNO’s relaxed
problem, which is given in equation (S-12) in Appendix E.
The DNO uses (21a)−(21d) to compute the updated prices
ρk+1 (t) and ̺k+1 (t) and penalties β k+1 (t) for all buses.
c) Load aggregator’s update: When load aggregator i
receives the control signal ρk+1
(h), h ∈ Ht from the DNO,
i
in Line 9, it updates its controllable load profile ek+1
(t)
i

6

Bus 17

1.1
0.95

0.8
7 am

12 pm

6 pm

12 am

0.95

1.6
Bus 23
1.5
1.4
1.3
1.2
6 am 7 am 12 pm
2.1

Bus 90

0.8

1.9

6 pm

12 am

6 am

6 pm

12 am

6 am

Bus 110

1.7

0.65
0.5
7 am

12 pm

6 pm

12 am

6 am

1.5
7 am

12 pm

Profit ($)

Figure 1. Load demand profiles over 24 hours in buses 17, 23, 90, and 110
with and without appliances scheduling.
2000
1750
1500
1250
1000
750
500

Without load scheduling
With load scheduling

17

23

90

110

Conventional
Generation (MW)

Figure 2. The profit for load aggregators with and without load scheduling.
25
20

Without Demand Response
With Demand Response

15
10
7 am

12 pm

6 pm

12 am

6 am

Renewable
Generation (MW)

Time (hour)
5
4
3
2
1
0

4.5
3
1.5
0

96

192

288

384

0
7 am 12 pm

478

6 pm

12 am

6 am

Time (hour)

Time (hour)

Figure 3. (a) The generation of the conventional unit of generator 15. (b)
The PV panel historical data samples. (c) The offered output power of the
PV panel of generator 15 in the market over the day.
Conventional
Generation (MW)

IV. P ERFORMANCE E VALUATION
In this section, we evaluate the performance of our proposed
decentralized algorithm on an IEEE 123-node test feeder. The
original test system is unbalanced. For all unbalanced case
studies, we construct a balanced test system by ignoring the
phase-to-phase admittance and replace all multi-phase lines
with a one-phase line with average inductance and resistance
of the phases. The data for the test system can be found in [27].
We consider the configuration, where all switched are open.
The slack bus is the substation bus. The trading horizon is one
day with H = 24 time slots. We add 10 generators at different
buses and assume that each load aggregator serves between
100 to 500. In Appendix F, we provide the simulation setup
[28], [29]. For the benchmark scenario, we consider a system
without demand response program for all load aggregators;
thus, users operate their appliances right after they become
awake. We perform simulations using Matlab R2016b in a PC
with processor Intel(R) Core(TM) i7-3770K CPU@3.5 GHz.
1) Load aggregators’ strategy: Each load aggregator executes Algorithm 1 to schedule the appliances of its users. Fig.
1 shows the load profile of the load aggregators in buses 17,
23, 90, and 110 in the benchmark scenario and the scenario
with load scheduling. Peak shaving can be observed in the
load profiles. Since Algorithm 1 is executed in real-time, the
load aggregators can only modify the demand for upcoming
time slots using the revealed information about the awake
appliances in the current time slot and the estimated load
demand for future time slots. Results for all load aggregators
verify that by executing Algorithm 1, the peak load demand is
reduced by 14.5% in average. Load scheduling is performed
by each load aggregator with the goal of increasing the profit
in (8). Fig. 2 shows that the profit of the load aggregators 17,
23, 90, and 110 is increased. Specifically, results show that
the profit for all load aggregators is increased by 17.8% on
average, since they can benefit from the price fluctuations by
modifying the operation of their users’ appliances.
2) Generators’ strategy: On the other hand, generators can
benefit from the users’ load scheduling to reduce their peak
generation, and thus their generation cost during peak hours.
For example, Fig. 3 (a) shows the active output power profile
from the conventional unit of the generator in bus 15. The
peak generation level is reduced from 25 MW to 20 MW (i.e.
20% reduction). Generator 15 also has a PV panel with the

1.25

Renewable
Generation (MW)

by solving its local problem (9), which is convex and can
be efficiently solved at each iteration. Note that the utility
function in (8) is a concave function.
d) Generator’s update: When generator j receives signals
ρk+1
(t), ̺k+1
(t) and βjk+1 (t), it updates its generation profile
j
j
conv,k
k
ψj (t) = (pj
(t), qjconv,k (t), pren,k
(t)) and decision variable
j
k
∆j (t), by solving its local problem in (13). This problem
is a linear problem and can be solved efficiently by the
generator using its local information about its conventional
and renewable units.
e) Step size update: We use a nonsummable diminishing
step size. In Line 17, the step size is updated.
We emphasize that in Algorithm 1, the DNO needs to
consider on bus (e.g., the substation bus) as the slack bus.

22

Without Demand Response
With Demand Response

19
16
13
10
7 am

12 pm

6 pm

12 am

6 am

Time (hour)
7.5
6
4.5
3
1.5

6
4.5
3
0

96

192

288

Time (hour)

384

478

1.5
7 am 12 pm

6 pm

12 am

6 am

Time (hour)

Figure 4. (a) The generation of the conventional unit of generator 60. (b)
The wind turbine historical data samples. (c) The presumed output power of
the wind turbine of generator 60.

historical generation record shown in Fig. 3 (b). The generator
executes Algorithm 1 and responds to the penalties βj (t) from
the DNO in each time slot to set the least-risk generation level
for its PV unit. Fig. 3 (c) shows the offers for generation 15
over the day. Note that the offers may not be equal to the actual
PV panel’s generation in real-time, but the generation profile
in Fig. 3 (c) results in the optimal risk of energy shortage for
PV panel in bus 15. We also show the generation profile of the
conventional unit of generator 60 in Fig. 4 (a). The reduction
in peak generation can be observed. This generator has wind
turbine with the historical generation record shown in Fig. 4

7

Wind Turbine
Generation (MW)

PV Panel
Generation (MW)

c

5
4
3
2
1
0
7 am

c

=5

12 pm

c

7.5

c

=10

6 pm

c

=10

c

=20

=50

12 am

c

=20

6 am

c

=50

=5

6
4.5
3
1.5
7 am

12 pm

6 pm

12 am

6 am

Time (hour)

Figure 5. (a) The offered output power of the PV panel of generator 15 with
different values of coefficient ϑc . (b) The offered output power of the wind
turbine of generator 60 with different values of coefficient ϑc .

PAR

1.5
1.4
1.3
1.2
1.1

1

15

18

36

51

60

79

86

100

108

Generator bus number

Figure 6. The PAR in the generation of the generators with and without
demand response.

Profit ($)

30

103

25
20
15
10

1

15

18

36

51

60

79

86

100

108

Figure 7. The profit of the generators with and without demand response.

(b). Fig. 3 (c) shows the offers for wind turbine in bus 60.
The offers for renewable units’ generation mainly depends
on the conservativeness of the DNO. In particular, when the
weight coefficient ϑc in (20a) is large, the DNO is risk-averse
and forces the generators to prevent generation shortage in
their renewable units. On the other hand, small coefficient
ϑc means the DNO encourages the generators to offer higher
amount of renewable generations. As an example, Figs. 5 (a)
and (b) show the renewable generation profiles of generators
15 and 60 for different values of ϑc . When ϑc increases from 5
to 50, the generation levels decrease, since the penalties βj (t)
in (21d) increase. Hence, the generators offer lower renewable
generations to reduce their cost of generation shortage.
To quantify the peak shaving, we consider the PAR of the
generation. Fig. 6 shows that the PAR is reduced for the
generators by 13% on average. A lower PAR means a lower
generation cost, and thus a higher profit. Fig. 7 confirms that
the generators’ profit is increased by 10.3% on average.
3) Algorithm convergence: We study the required number
of iterations for convergence, which can be interpreted as an
indicator of the number of message exchanges among the load
aggregators, generators and DNO. The angle and magnitude
of the voltage of the buses depend on all generators’ and
load aggregators’ decision variables. Thus, the convergence
of the these variables is a viable indicator of the convergence

of all decision variables in the system. Since the values of
the voltage angles of all buses can be added by a constant,
we illustrate the convergence of the phase angle difference
between the voltages of the buses at the end nodes of the
lines. As an example, we consider time slot 1 and we provide the convergence of δ14 (1) − δ11 (1), δ70 (1) − δ71 (1),
δ34 (1) − δ13 (1), and δ72 (1) − δ76 (1) in Fig. 8 (a). We also
show the voltage magnitude of the buses 19, 38, 76, and 30
in Fig. 8 (b). We can observe that 50 iterations are enough
for convergence. The average running time of the algorithm
for different initial conditions is 5 seconds for 100 random
initial conditions. The low convergence rate and running time
make the proposed algorithm implementable for real-time
interactions among entities in an energy market.
We also evaluate the running time of Algorithm 1 for larger
test systems to show its scalability. Meanwhile, we compare
its running time with a centralized algorithm, where the DNO
solves problem problem (20). We use MOSEK solver to
solve the DNO’s centralized problem (20). The control signals
in (21a)−(21d) are determined in order to obtain the same
solution for both the centralized and decentralized algorithms,
and simulation results confirm this. We provide the average
running time of Algorithm 1 and the centralized approach for
six test systems (all can be found in [27] except the system
with 1500 buses, which is a part of 8500-bus test system) in
Fig. 9. We can observe the the centralized algorithm suffer
from a high running time due to a large number of decision
variables and constraints. On the other hand, Algorithm 1
is executed by each entity to solve its own optimization
problem with its locally available information in a distributed
fashion. Hence, the number of decision variables for each
entity becomes independent of the size of the test system. The
overall running time of Algorithm 1 increases almost linearly
with the number of buses due to the increase in the required
number of iterations for convergence in larger test systems.
We also compare the performance of Algorithm 1 for the
scenario with uncertainty in the load demand and renewable
generations and the scenario with complete information. As an
example, we consider the load profile of the load aggregator
110 and the generation profile of generator 15 with a PV panel
in Figs. 10 (a) and (b). The lack of information makes the load
aggregators more conservative, since it considers the worstcase for the electric appliances in the upcoming time slots.
Whereas, when the load aggregator has complete information,
it can better manage the electric appliances especially during
the peak hours. A lower peak demand for the load aggregators
results in a lower peak in the generation level of the conventional units. The conventional unit may provide more power
during the times when the PV panel has high generation (e.g.
around 6 pm), as the PV generation is known and the generator
does not take risk to offer a high renewable generation.
V. C ONCLUSION
In this paper, we proposed a real-time decentralized algorithm for energy trading among load aggregators and generators. Our proposed approach considers the uncertainty at
both generation and demand sides. In our model, the DNO
sends control signals to the entities to encourage them towards

8

Running Time (s)

Figure 8. The convergence of phase differences and voltage magnitudes.
1000
750

Centeralized Algorithm
Decenteralized Algorithm

500
250
0

13-bus

34-bus

37-bus

123-bus

1300-bus

8500-bus

Test System

2.1
1.9
1.7
1.5
7 am

Conventional
Generation (MW)

Load Demand

Figure 9. The running time of the centralized and decentralized algorithms.

25

12 pm

6 pm

Without Demand Response

With Uncertainty

12 am

6 am

Without Uncertainty

20
15
10
7 am

12 pm

6 pm

12 am

6 am

Time (hour)

Figure 10. (a) Load scheduling with uncertainty and complete information.
(b) Generation profile with uncertainty and complete information.

optimizing their objectives independently, while meeting the
physical constraints of the power network. This study uses
linearized ac optimal power flow formulation to increase the
accuracy of the obtained solution. Further, we solve the problem in a real-time fashion to obtain the most recent optimal
solution for each time-step. To evaluate the performance of the
proposed decentralized algorithm, we used an IEEE 123-bus
test feeder connected to some renewable generators. Although
we considered the stochastic nature of renewables as well as
ac power flow formulation, our algorithm converges in 50
iterations. We evaluated the price responsive load profiles and
generation values and showed that our method can benefit
the load aggregators by increasing their profit by 17.8%, and
the generators by reducing the PAR by 13% and increasing
their profit by 10.3%. Our algorithm benefit the DNO by
maintaining the privacy issues and a lower computational time
compared to the centralized approach.
R EFERENCES
[1] B. Chai, J. Chen, Z. Yang, and Y. Zhang, “Demand response management with multiple utility companies: A two-level game approach,” IEEE
Trans. on Smart Grid, vol. 5, no. 2, pp. 722–731, Mar. 2014.
[2] S. Maharjan, Q. Zhu, Y. Zhang, S. Gjessing, and T. Basar, “Dependable
demand response management in the smart grid: A stackelberg game
approach,” IEEE Trans. on Smart Grid, vol. 4, no. 1, pp. 120–132, Mar.
2013.
[3] R. Deng, Z. Yang, F. Hou, M.-Y. Chow, and J. Chen, “Distributed realtime demand response in multiseller-multibuyer smart distribution grid,”
IEEE Trans. on Power Systems, vol. PP, no. 99, pp. 1–11, Oct. 2014.
[4] F. Kamyab et al., “Demand response program in smart grid using supply
function bidding mechanism,” IEEE Trans. on Smart Grid, vol. 7, no. 3,
pp. 1277 – 1284, 2016.

[5] M. Parvania, M. Fotuhi-Firuzabad, and M. Shahidehpour, “ISO’s optimal
strategies for scheduling the hourly demand response in day-ahead
markets,” IEEE Trans. on Power Systems, vol. 29, no. 6, pp. 2636–2645,
2014.
[6] L. Gan, N. Li, U. Topcu, and S. H. Low, “Exact convex relaxation
of optimal power flow in radial networks,” IEEE Trans. on Automatic
Control, vol. 60, no. 1, pp. 72–87, 2015.
[7] W. Shi, N. Li, X. Xie, C. Chu, and R. Gadh, “Optimal residential demand
response in distribution networks,” IEEE Journal on Selected Areas in
Comm., vol. 32, no. 7, pp. 1441–1450, Jun. 2014.
[8] N. Li, L. Gan, L. Chen, and S. Low, “An optimization-based demand
response in radial distribution networks,” in Proc. of IEEE Globecom,
Anaheim, CA, Anaheim, CA 2012.
[9] E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal power
flow for smart microgrids,” IEEE Trans. on Smart Grid, vol. 4, no. 3,
pp. 1464–1475, Sept. 2013.
[10] A. G. Bakirtzis and P. N. Biskas, “A decentralized solution to the DCOPF of interconnected power systems,” IEEE Trans. on Power Systems,
vol. 18, no. 3, pp. 1007–1013, Aug. 2003.
[11] T. Erseghe, “Distributed optimal power flow using ADMM,” IEEE
Trans. on Power Systems, vol. 29, no. 5, pp. 2370–2380, Sept. 2014.
[12] Q. Peng and S. H. Low, “Distributed algorithm for optimal power flow
on a radial network,” in Proc. of IEEE Conf. on Decision and Control,
Dec. 2014, pp. 167–172.
[13] S. Magnusson, P. C. Weeraddana, and C. Fischione, “A distributed
approach for the optimal power-flow problem based on ADMM and
sequential convex approximations,” IEEE Trans. on Control of Network
Systems, vol. 2, no. 3, pp. 238–253, Sept. 2015.
[14] J. M. Arroyo and F. D. Galiana, “Energy and reserve pricing in security
and network-constrained electricity markets,” IEEE Trans. on Power
Systems, vol. 20, no. 2, pp. 634–643, May 2005.
[15] S. Mhanna, A. C. Chapman, and G. Verbič, “A fast distributed algorithm
for large-scale demand response aggregation,” IEEE Trans. on Smart
Grid, vol. 7, no. 4, pp. 2094–2107, 2016.
[16] M. J. Dolan, E. M. Davidson, I. Kockar, G. W. Ault, and S. D. McArthur,
“Distribution power flow management utilizing an online optimal power
flow technique,” IEEE Trans. on Power Systems, vol. 27, no. 2, pp. 790–
799, May 2012.
[17] E. Belic, N. Lukac, K. Dezelak, B. Zalik, and G. Stumberger, “Gpubased online optimization of low voltage distribution network operation,”
accepted for publication in IEEE Trans. on Smart Grid, 2017.
[18] L. Gan and S. H. Low, “An online gradient algorithm for optimal power
flow on radial networks,” IEEE J. on Selected Areas in Comm., vol. 34,
no. 3, pp. 625–638, Mar. 2016.
[19] A. Hauswirth, S. Bolognani, G. Hug, and F. Dörfler, “Projected gradient
descent on Riemannian manifolds with applications to online power
system optimization,” in Proc. of Allerton Conf. on Communications,
Control and Computing, Sept. 2016, pp. 225–232.
[20] S.-J. Kim, G. B. Giannakis, and K. Y. Lee, “Online optimal power flow
with renewables,” in Proc. of Asilomar Conf. on Signals, Systems and
Computers, Nov. 2014, pp. 355–360.
[21] D. Mehta, A. Ravindran, B. Joshi, and S. Kamalasadan, “Graph theory
based online optimal power flow control of power grid with distributed
flexible ac transmission systems (d-facts) devices,” in Proc. of North
American Power Symposium (NAPS), Oct. 2015, pp. 1–6.
[22] N. Li, L. Chen, and S. H. Low, “Optimal demand response based on
utility maximization in power networks,” in IEEE Power and Energy
Society General Meeting, Jul. 2011, pp. 1–8.
[23] T. Li and M. Shahidehpour, “Price-based unit commitment: a case of
lagrangian relaxation versus mixed integer programming,” IEEE Trans.
on Power Systems, vol. 20, no. 4, pp. 2015–2025, 2005.
[24] D. Bertsimas, E. Litvinov, X. A. Sun, J. Zhao, and T. Zheng, “Adaptive
robust optimization for the security constrained unit commitment problem,” IEEE Trans.s on Power Systems, vol. 28, pp. 52–63, Feb. 2013.
[25] G. Hug-Glanzmann and G. Andersson, “Decentralized optimal power
flow control for overlapping areas in power systems,” IEEE Trans. on
Power Systems, vol. 24, no. 1, pp. 327–336, Feb. 2009.
[26] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge
University Press, 2004.
[27] [Online]. Available: https://ewh.ieee.org/soc/pes/dsacom/testfeeders/
[28] Independent Electricty System Operator (IESO). [Online]. Available:
http://www.ieso.ca
[29] [Online]. Available: http://www.torontohydro.com/sites/electricsystem/\residential/yourb

