NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

Conference Paper

NREL/CP-5D00-76338

July 2020

Model-Free Primal-Dual Methods for 

Network Optimization with Application 

to Real-Time Optimal Power Flow

Preprint

Yue Chen,1 Andrey Bernstein,1 Adithya Devraj,2 and 

Sean Meyn2

1 National Renewable Energy Laboratory 

2 University of Florida, Gainesville 

Presented at the 2020 American Control Conference (ACC) 

July 1–3, 2020 
NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

National Renewable Energy Laboratory

15013 Denver West Parkway

Golden, CO 80401

303-275-3000 • www.nrel.gov

Conference Paper

NREL/CP-5D00-76338

July 2020

Model-Free Primal-Dual Methods for 

Network Optimization with Application 

to Real-Time Optimal Power Flow

Preprint

Yue Chen,1 Andrey Bernstein,1 Adithya Devraj,2 and 

Sean Meyn2 

1 National Renewable Energy Laboratory 

2 University of Florida, Gainesville 

Suggested Citation 

Chen, Yue, Andrey Bernstein, Adithya Devraj, and Sean Meyn. 2020. Model-Free PrimalDual Methods for Network Optimization with Application to Real-Time Optimal Power 

Flow: Preprint. Golden, CO: National Renewable Energy Laboratory. NREL/CP-5D00-

76338. https://www.nrel.gov/docs/fy20osti/76338.pdf. 
NOTICE

This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for 

Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. 

This work was supported by the Laboratory Directed Research and Development (LDRD) Program at NREL. The 

views expressed herein do not necessarily represent the views of the DOE or the U.S. Government.

This report is available at no cost from the National Renewable 

Energy Laboratory (NREL) at www.nrel.gov/publications. 

U.S. Department of Energy (DOE) reports produced after 1991 

and a growing number of pre-1991 documents are available 

free via www.OSTI.gov. 

Cover Photos by Dennis Schroeder: (clockwise, left to right) NREL 51934, NREL 45897, NREL 42160, NREL 45891, NREL 48097, 

NREL 46526.

NREL prints on paper that contains recycled content. 
Model-Free Primal-Dual Methods for Network Optimization with

Application to Real-Time Optimal Power Flow

Yue Chen1, Andrey Bernstein1, Adithya Devraj2, Sean Meyn2

Abstract— This paper examines the problem of real-time

optimization of networked systems and develops online algorithms that steer the system towards the optimal trajectory

without explicit knowledge of the system model. The problem

is modeled as a dynamic optimization problem with timevarying performance objectives and engineering constraints.

The design of the algorithms leverages the online zero-order

primal-dual projected-gradient method. In particular, the primal step that involves the gradient of the objective function

(and hence requires a networked systems model) is replaced

by its zero-order approximation with two function evaluations

using a deterministic perturbation signal. The evaluations are

performed using the measurements of the system output, hence

giving rise to a feedback interconnection, with the optimization

algorithm serving as a feedback controller. The paper provides

some insights on the stability and tracking properties of this

interconnection. Finally, the paper applies this methodology to

a real-time optimal power flow problem in power systems, and

shows its efficacy on the IEEE 37-node distribution test feeder

for reference power tracking and voltage regulation.

I. INTRODUCTION

This paper considers the problem of real-time optimization

of networked systems where the desired operation of the

network is formulated as a dynamic optimization problem

with time-varying performance objectives and engineering

constraints [1]–[3]. To illustrate the idea, consider the following time-varying model that describes the input-output

behavior of the systems in the network:

y(t) = ht(x(t)) (1)

where x(t) ∈ R

n is a vector representing controllable inputs;

y(t) ∈ R

m collects the outputs of the network; and ht(·) :

R

n → Rm is a time-varying map representing the system

model. Suppose that the desired behavior of the system is

defined via a time-varying optimization problem of the form

min

x∈X (t),y=ht(x)

ft(y) (2)

1Y. Chen and A. Bernstein are with the Power Systems Engineering

Center, National Renewable Energy Laboratory, Golden, CO 80401, USA

yue.chen@nrel.gov, andrey.bernstein@nrel.gov

2A. Devraj and S. Meyn are with the Department of Electrical and

Computer Engineering, University of Florida, Gainesville, FL 32611, USA

adithyamdevraj@ufl.edu, meyn@ece.ufl.edu

where X (t) is a convex set representing, e.g., engineering

constraints; and ft : R

m → R is a convex function

representing the time-varying performance goals.

In an ideal situation, when infinite computation and

communication capabilities are available, one could seek

solutions (or stationary points) of (2) at each time t by, e.g.,

an iterative projected-gradient method:

x

(k+1) =ProjX (t)

n

x

(k) − α(J

(k)

t

)

T∇yft(ht(x(k)

))o, (3)

where k = 0, 1, . . . is the iteration index; J

(k)

t

:=

Jht(x

(k)

) is the Jacobian matrix of the input-output map

ht; ProjX (z) := arg minx∈X kz − xk2 denotes projection;

and α > 0 is the step size. In practice, this approach might

be infeasible. First, running (3) to convergence might be

impossible due to stringent real-time constraints. Second, (3)

requires the functional knowledge of the input-output map

ht or an approximation thereof in real time. The latter is

absent in typical network optimization applications, such as

the power distribution system example considered here.

To address the first challenge, we adopt the measurementbased optimization framework of [3]–[6] wherein (3) is

replaced with a running (or online) version:

x

(k+1) =ProjX (k)

n

x

(k) − α(J

(k)

)

T∇yf(k)

(yb

(k)

)

o

, (4)

in which the iteration index k is now identified with a discrete

time instance tk at which the computation is performed;

X

(k)

:= X (tk), J

(k)

:= Jhtk(x

(k)

), f

(k)

:= ftk; and yb

(k)

represents the measurement of the network output htk(x

(k)

)

at time tk. Note that (4) represents a feedback controller

that computes the input to the system based on its output;

however, the implementation of this controller requires the

real-time model information in the form of J

(k)

.

This paper sets out to address the second challenge and

to tackle situations in which J

(k)

is not available for realtime optimization. This is motivated by networked systems

optimization problems, such as the optimal power flow in

power systems considered in Section V, wherein the topology and parameters of the network, as well as exogenous

inputs, might change rapidly due to variability of working

conditions, natural disasters, and cyber-physical attacks. To

address this challenge, we leverage a zero-order approximation of the gradient of the objective function. In particular,

we replace the gradient of1 F

(k)

(x) := f

(k)

(htk(x))

∇F

(k)

(x) = (Jhtk

(x))T∇yf

(k)

(htk(x)) (5)

1A more direct approach is to estimate the Jacobian matrix with a zeroorder approximation. This direction is a subject of ongoing research.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 1
with the two function evaluation approximation:

∇bF

(k)

(x; ξ, ) := 1

2

ξ

h

F

(k)

(x + ξ) − F

(k)

(x − ξ)

i

(6)

where ξ ∈ R

n is a perturbation (or exploration) vector, and

 > 0 is a (small) scalar. Observe that this modification now

requires measurements of the output at two (rather than one)

inputs: x

(k) ±ξ at each time step k. With this modification,

algorithm (4) can be carried out without any knowledge of

the map ht.

Although the stylized algorithm (4) and approximation (6)

are used here to illustrate the main ideas, the paper considers

a more general convex constrained optimization framework,

and develops model-free primal-dual methods to track the

optimal trajectories. Using a deterministic exploration approach reminiscent of the quasi-stochastic approximation

method [7], this paper provides design principles for the

exploration signal ξ, as well as other algorithmic parameters,

to ensure stability and tracking guarantees. In particular, we

show that under some conditions, the iterates x

(k)

converge

within a ball around the optimal solution of (2). Finally,

we apply the developed methodology to the problem of

distributed real-time voltage regulation and power setpoint

tracking in power distribution networks.

Much of the analysis is restricted to the case in which the

model map ht is linear. Extensions to nonlinear models is

a topic of current research.

Literature Review

Zero-order (or gradient-free) optimization has been a

subject of interest in the optimization, control, and machine

learning communities for decades. The seminal paper of

Kiefer and Wolfowitz [8] introduced a one-dimensional

variant of approximation (6); for a d-dimensional problem, it perturbs each dimension separately and requires 2d

function evaluations. The simultaneous perturbation stochastic approximation (SPSA) algorithm [9] uses zero-mean

independent random perturbations, requiring two function

evaluations at each step. In [10], [11], the SPSA algorithm

was extended to deterministic perturbations, to improve

convergence rates under the assumption of a vanishing stepsize and vanishing quasi-noise.

Analysis of the zero-order methods based on two function

evaluation (similar to (6)) has recently been a focus of several

papers [12]–[15]. In these papers, a stochastic exploration

signal ξ is used (typically, Gaussian iid), and convergence

analysis of (projected) gradient methods and various variants

of the method of multipliers are offered.

In the theoretical machine learning community, “bandit

optimization” refers to zero-order online optimization algorithms based on a single or multiple evaluations of the

objective function [16]–[18]. These algorithms typically have

high variance of the estimate if the evaluation is performed

using one or two function evaluations [19]. Finally, the

gradient-free technique termed “extremum-seeking control”

(ESC) [20] leverages sinusoidal perturbation signals and

single function evaluation to estimate the gradient. Stability

of the ESC feedback scheme was analyzed in, e.g., [21],

[22]. These algorithms typically require additional filters to

smooth out the noise introduced by the perturbations.

In terms of zero-order network optimization, our work is

closely related to the recent work [15], which considers a

multi-agent optimization of a class of nonconvex problems.

The paper then develops distributed algorithms based on the

zero-order approximation of the method of multipliers. In

contrast to our paper, [15] considers a stochastic exploration

signal for the gradient estimation, and typically requires

multiple function evaluations to reduce the estimation variance [15, Lemma 1]; see also [7] for the detailed analysis

of the advantage of deterministic vs. stochastic exploration.

Moreover, [15] considers only a static optimization problem.

Contrary to much previous work (with the exception of

perhaps [15]), we consider a constant step-size algorithm, necessitated by our application in a time-varying optimization

problem. Moreover, because the problem we consider here is

a constrained optimization problem, it leads to the extension

of the gradient-free methods to novel model-free primal-dual

algorithms. Due to the real-time nature of our application, we

resort here to computationally light online gradient-descent

methods rather than variants of methods of multipliers that

typically require solving an optimization problem at each

time step.

Organization: The rest of the paper is organized as

follows. Section II introduces the general network optimization problem. Section III proposes the model-free primaldual algorithm to solve the problem in Section II. The

convergence of the algorithm is analyzed in Section IV.

Section V introduces the power systems application, and the

simulation results are presented in Section VI. Section VII

concludes the paper and proposes future works.

II. NETWORK OPTIMIZATION FRAMEWORK

Consider the general network optimization problem [3],

wherein at each time step k ∈ N, the goal is to optimize the

operation of a network of N systems:

min

x∈R

n

f

(k)

0

(y

(k)

(x)) +X

N

i=1

f

(k)

i

(xi) (7a)

subject to : xi ∈ X (k)

i

, i = 1, . . . , N (7b)

g

(k)

j

(y

(k)

(x)) ≤ 0, j = 1, . . . , M (7c)

where X

(k)

i ⊂ R

ni

, with PN

i=1 ni = n, are convex sets

representing operational constraints on the control input xi

of system i; and y

(k)

(x) is an algebraic representation of

some observable outputs in the networked system. In this

paper, much analysis is restricted to the linear case2

y

(k)

(x) := Cx + Dw(k) ∈ R

m (8)

where C ∈ R

m×n and D ∈ Rm×w are given model

parameters, and w(k) ∈ R

w is a vector of time-varying

exogenous (uncontrollable) inputs. For example, in the power

2However, the developed algorithms, as well as our numerical results in

Section VI, are not restricted to the linear case.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 2
systems case considered in Section VI, y

(k)

(x) represents

the linearized power flow equations.

Further, the convex function f

(k)

0

(y

(k)

(x)) : R

n → R

represents the time-varying performance goals associated

with the outputs y

(k)

(x); whereas f

(k)

i

(xi) : R

ni → R

is a convex function representing the performance goals of

system i. Finally, g

(k)

j

(y

(k)

(x)) : R

m → R are convex

functions used to impose time-varying constraints on the

output y

(k)

(x).

Observe that (7) is naturally model-based—in order to

solve it, the knowledge of the (linearized) system model (8)

and the exogenous input w(k)is required. This problem will

be used next to define the desired operational trajectories for

the networked system. In Section III, we propose a modelfree algorithm to track these desired operational trajectories.

We define the following shorthand notation:

g

(k)

(y

(k)

(x)) := hg

(k)

1

(y

(k)

(x)), . . . , g

(k)

M (y

(k)

(x))iT

f

(k)

(x) := X

i

f

(k)

i

(xi)

h

(k)

(x) := f

(k)

(x) + f

(k)

0

(y

(k)

(x)).

Setting λ ∈ RM

+ as the vector of the dual variables associated

with (7c), the Lagrangian function at step k is given by:

L

(k)

(x, λ) := h

(k)

(x) + λ

T

g

(k)

(y

(k)

(x)). (9)

Finally, the regularized Lagrangian function is given by [23]:

L

(k)

p,d(x, λ) := L

(k)

(x, λ) + p

2

kxk

2

2 −

d

2

kλk

2

2

(10)

where p > 0 and d > 0 are the regularization parameters.

With this notation at hand, consider the saddle-point problem

associated with L

(k)

p,d:

max

λ∈D(k)

min

x∈X (k)

L

(k)

p,d(x, λ) k ∈ N (11)

where X

(k)

:= X

(k)

1 × . . . × X (k)

N and D(k)

is a convex

and compact set that contains the optimal dual variable of

(11); see, e.g., [3] for details. Hereafter, the optimal trajectory

z

(∗,k)

:= {x

(∗,k)

,λ

(∗,k)

}k∈N of (11) represents the desired

trajectory. It is clear that L

(k)

p,d(x, λ) is strongly convex in

x and strongly concave in λ; thus, the optimizer z

(∗,k)

of (11) is unique. Moreover, algorithms based on L

(k)

p,d(x, λ)

exhibit Q-linear (i.e., exponentially fast) convergence, which

is the major motivation for the regularization. However,

z

(∗,k) might be different from any saddle point of the

exact Lagrangian L

(k)

(x, λ), and the bound on the distance

between z

(∗,k)

and the solution of (7) can be established as,

e.g., in [23].

To proceed, we list the assumptions that will be imposed

on the different elements of (7).

Assumption 1. Slater’s condition holds at each time step k.

Assumption 2. The functions f

(k)

0

(y) and f

(k)

i

(xi) are

convex and continuously differentiable. The maps ∇f

(k)

(x),

∇f

(k)

0

(y), and ∇2f

(k)

0

(y) are Lipschitz continuous.

Assumption 3. For each j = 1, . . . , M, the function

g

(k)

j

(y) is convex and continuously differentiable. Moreover,

∇g

(k)

j

(y) and ∇2g

(k)

j

(y) are Lipschitz continuous.

Assumption 4. There exists a constant ef such that for all

k and x, y:

k∇f

(k)

0

(y) − ∇f

(k−1)

0

(y)k ≤ ef ,

k∇f

(k)

(x) − ∇f

(k−1)(x)k ≤ ef ,

k∇g

(k)

j

(y) − ∇g

(k−1)

j

(y)k ≤ ef , j = 1, . . . , M.

III. MODEL-FREE PRIMAL-DUAL METHOD

Leveraging the regularized Lagrangian L

(k)

p,d, a

measurement-based primal-dual method was introduced in

[3] to track the desired trajectory {z

(∗,k)}; we present the

method here for convenience:

x

(k+1) = ProjX (k)

n

(1 − αp)x

(k) − α



∇xf

(k)

(x

(k)

)

+ C

T∇yf

(k)

0

(yb

(k)

) +X

M

j=1

λ

(k)

j C

T∇g

(k)

j

(yb

(k)

)

o

(12a)

λ

(k+1) = ProjD(k)

n

(1 − αd)λ

(k) + αg

(k)

(yb

(k)

)

o

(12b)

where α > 0 is a constant step size, and yb

(k)

is a (possibly

noisy) measurement of y

(k)

(x

(k)

) collected at time step k

(see Assumption 5). The main advantage of (12) is that it

avoids explicit computation of (8) at each time step, and

thus it does not require knowledge or estimation of the

uncontrollable input w(k)appearing in (8). Note that because

the method is gradient-based, it relies on the availability of

the network model matrix C. However, in many practical

applications, including the power system example in Section VI, C is time-varying (e.g., due to network topology

variation), and its exact value is unavailable.

We next propose a model-free variant of (12). As in the

celebrated gradient-free stochastic approximation algorithm

of Kiefer and Wolfowitz [8], the general idea is to use few

(online) evaluations of the objective function to estimate the

gradient. The gradient approximation is justified through a

Taylor series expansion as explained next.

Let F : R

n → R be a C3

function, whose first and second

derivatives are globally Lipschitz continuous (cf. Assumption

2). For given x ∈ R

n, this variable is perturbed by ±ξ, with

ξ ∈ R

n and  > 0. According to Taylor’s theorem, for any

x ∈ R

n and r ∈ R,

F(x + rξ) = F(x) + rξ

T∇F(x) + r

2

2

ξ

T∇2F(x)ξ + O(r

3

)

Taking r = , r = −, and then subtracting yields:

F(x + ξ) − F(x − ξ) = 2ξ

T∇F(x) + O(

3

) (13)

This approximation is summarized in Lemma 1 that follows, with the following shorthand notation:

∇bF(x; ξ, ) := 1

2

ξ [F(x + ξ) − F(x − ξ)] . (14)

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 3
Lemma 1. Let F : R

n → R be a C3

function, with Lipschitz

continuous ∇F and ∇2F. Then,

∇bF(x; ξ, ) = ξξT ∇F(x) + O(

2

) (15)

Note that O(

2

) = 0 if F is a quadratic function.

In applying Lemma 1 to algorithm design, the vector ξ is

updated at each iteration of the algorithm. Let {ξ

(k)

} denote

a deterministic exploration signal, with ξ

(k) ∈ R

n. It is

assumed in this paper that the exploration signal is obtained

by sampling the sinusoidal signal

ξi(t) = √2 sin(ωit), i = 1, . . . , n, (16)

with ωi 6= ωj for all i 6= j; other signals are possible

provided that they satisfy Assumption 6.

The following gradient-free variant of the primal-dual

algorithm (12) is well motivated by Lemma 1.

Model-Free Primal-Dual Algorithm

At each time step k, perform the following steps:

[S1a] (forward exploration): Apply x

(k)

+ := x

(k) +ξ

(k)

to

the system, and collect the measurement yb

(k)

+ of the output

y

(k)

(x

(k)

+ ).

[S1b] (backward exploration): Apply x

(k)

− := x

(k) − ξ

(k)

to the system, and collect the measurement yb

(k)

− of the output

y

(k)

(x

(k)

− ).

[S1c] (control application): Apply x

(k)

to the system, and

collect the measurement yb

(k)

of the output y

(k)

(x

(k)

).

[S2a] (approximate gradient): Compute the approximate

gradient for the primal step

∇Lb (k):= ∇xf

(k)

(x

(k)

)

+

1

2

ξ

(k)

h

f

(k)

0

(yb

(k)

+ ) − f

(k)

0

(yb

(k)

− )

i

+

1

2

ξ

(k)

(λ

(k)

)

T

h

g

(k)

(yb

(k)

+ ) − g

(k)

(yb

(k)

− )

i

. (17)

[S2b] (approximate primal step): Compute

x

(k+1) = ProjX (k)

n

(1 − αp)x

(k) − α∇Lb (k)

o

. (18)

[S3] (dual step): Compute

λ

(k+1) = ProjD(k)

n

(1 − αd)λ

(k) + αg

(k)

(yb

(k)

)

o

. (19)

The following remarks are in order.

Remark 1. Lemma 1 alone is not enough to obtain convergence guarantees because the matrix ξξT

is semidefinite but

not strictly definite. We formulate conditions for convergence

in the next section.

Remark 2. The algorithm allows distributed implementation. This point is illustrated in a power system optimization

example described in Section VI.

Remark 3. The algorithm requires three measurements of

the output: yb

(k)

+ , yb

(k)

− , and yb

(k)

. In fact, the third measurement can be replaced with, e.g., an average yb

(k)

:=

0.5(yb

(k)

+ + yb

(k)

− ), provided that yb

(k)

satisfies Assumption 5.

IV. CONVERGENCE ANALYSIS OUTLINE

For this analysis, we assume bounded measurement error.

Assumption 5. There exists a scalar ey < ∞ such that the

measurement error can be bounded as

sup

k≥1

kyb

(k) − y

(k)

(x

(k)

)k2 ≤ ey,

sup

k≥1

kyb

(k)

+ − y

(k)

(x

(k)

+ )k2 ≤ ey, (20)

sup

k≥1

kyb

(k)

− − y

(k)

(x

(k)

− )k2 ≤ ey.

The exploration signal plays an important role in the convergence of the proposed algorithm. Recent work on quasistochastic approximation [7], [24] and extremum seeking

control [25] has shown that certain deterministic exploration

signals can help reduce the asymptotic variance and improve

convergence, compared to random exploration (see Fig. 1 of

[24]). We impose the following assumption on the exploration signal.

Assumption 6. The exploration signal is a periodic signal

with period T, and

1

T

Z t+T

t

ξ(τ )ξ(τ )

T

dτ = I (21)

for all t, where I is the identity matrix.

It can be shown that the signal defined in (16) satisfies

Assumption 6 with T being a common integer multiple of

the sinusoidal signal periods.

This strong assumption is imposed to simplify the proofs

that follow. It is enough to have the approximation:

1

T

Z t+T

t

ξ(τ )ξ(τ )

T

dτ = Σξ + O(1/T)

in which Σξ is positive definite, and the 1/T bound for the

error term is independent of t (this is the assumption used

in [7]). An example is the signal defined in (16) in which

the frequencies might not have a common multiple.

In this paper, the analysis is performed on a simplified iteration (18), wherein the projection is performed periodically,

every T /∆t time steps, with ∆t being the time discretization

step. The analysis of the original iteration is left for future

research.

Due to space constraints, the detailed analysis appears in

the full version [26]. It amounts to showing that (18) is

equivalent to the following approximate primal iteration:

x

(k+1) = Projk

n

x

(k) − α



∇xL

(k)

p,d(z

(k)

)

+ O(α + 

2 + ef + ey)

o, (22)

where z

(k)

:= {x

(k)

,λ

(k)

}, and Projkis the projection onto

X

(k)

for k a multiple of T /∆t.

The analysis then proceeds similarly to [3]. To this end,

quantify the temporal variability of the desired trajectory

z

(∗,k) via

σ

(k)

:= kz

(∗,k+1) − z(∗,k)

k2 . (23)

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 4
In addition, let

φ

(k)

: z 7→

"

∇xL

(k)

p,d(z)

−g

(k)

(y

(k)

(x)) + dλ

#

(24)

denote the primal-dual operator associated with (11). The

following holds (see, e.g., [3, Lemma 5]).

Lemma 2. There exist constants 0 < ηφ ≤ Lφ < ∞ such

that, for every k ∈ N, the map φ

(k)

(z) is strongly monotone

over X

(k)×D(k) with constant ηφ and Lipschitz over X(k)×

D(k) with coefficient Lφ.

Finally, we use [3, Theorem 4] to show that the approximate time-varying primal-dual algorithm (22)+(19) exhibits

the following tracking performance.

Theorem 1. Suppose that the step size α is chosen such

that 0 < α < 2ηφ/L2

φ

. Moreover, suppose that there exists

a scalar σ < ∞ such that supk≥1 σ

(k) ≤ σ. Then, the

sequence {z

(k)} converges Q-linearly to {z(∗,k)} up to an

asymptotic error bound given by:

lim sup

k→∞

kz

(k) − z(∗,k)

k2 ≤

αφ + σ

1 − c

(25)

where c := q

1 − 2αηφ + α2L

2

φ < 1 and φ = O(α + 

2 +

ef + ey).

V. MODEL-FREE REAL-TIME OPTIMAL POWER FLOW

The application considered in this paper is real-time

optimization of the power injections of distributed energy

resources (DERs) in a distribution system. The following

objectives are pursued: (i) the feeder head/substation power

is following a prescribed signal, and (ii) the voltages are

regulated within prescribed bounds. This is achieved via

formulating an optimal power flow (OPF) problem of the

form (7) and seeking its solutions in real time using the

model-free primal-dual framework introduced in Section III.

Consider a distribution system that is described by a

set of nodes N = {0, . . . , n}, where node 0 denotes the

feeder head (or substation) with fixed voltage; the rest

are P Q nodes. In this work, we focus on controlling the

batteries and photovoltaic (PV) systems in the distribution

system interfaced via power inverters; however, the proposed

framework can be easily extended to include other types of

DERs. Let Nbt and Npv be the number of batteries and PV,

respectively. So the total number of DERs Nder = Nbt+Npv.

A. Notation

For notation brevity, in this section, we drop the algorithm

iteration superscript (k)from variable names, unless otherwise specified.

The distribution system is described by the AC power flow

equations, defining the input-output relationship (v, P0) =

h(x, `). In particular:

• x ∈ R

2Nder collects all decision variables, where each

DER has real and reactive power xi = {xi,p, xi,q}

• ` ∈ R

2n collects the active and reactive power of

uncontrollable loads (the “noise”)

• v(x, `) ∈ R

n collects the voltage magnitudes at every

P Q node

• P0(x, `) ∈ R is the power flow at the feeder head.

The maps (v(x, `), P0(x, `)) exist under some conditions

and are defined by the power flow equations [27]. The goal

is to control the output (v(x, `), P0(x, `)) by controlling x.

To formulate a convex target OPF (7), an approximate

linearized relationship between the output y := (v, P0) and

the input (x, `) is leveraged. In particular, at each time step

k, the output is approximated as

y = Cx + D` + y0; (26)

see, e.g., [27]. Note that the relationship (26) is used only to

define the desired trajectory via (7) and (11); the actual simulation in Section VI is performed using the exact (nonlinear)

power flow equations.

B. Objective Function

The objective function (7a) consists of two parts:

f0(x) =(P0(x, `) − P

•

0

)

2

(27a)

fi(xi) =ci,px

2

i,p + ci,qx

2

i,q + ci,p• (xi,p − x

•

i,p)

2

(27b)

where ci,p, ci,q, and ci,p• are weighting parameters for DER

i. The global cost (27a) drives the substation real power P0

to follow the reference power P

•

0

, which is usually timevarying. The local DER cost (27b) penalizes the control

effort and the difference between the control and its desired

reference power. For a PV system, the reference power x

•

i,p

is the real-time maximal PV power generation; for a battery

system, the reference x

•

i,p is the power to reach the preferred

state of charge (SOC), which is usually set as the middle

point of the SOC to allow maximal operation flexibility.

C. Constraints

Constraints associated with the OPF are summarized as

follows.

1) Voltage: Each node needs to maintain a stable voltage

magnitude within a certain range. For any node i ∈ N ,

V i ≤ vi(x, `) ≤ V i (28)

2) Battery: Battery limits are imposed on battery control

variables to enable feasible operations. Each battery has its

power and energy limits:

Xi,p ≤ xi,p ≤ X i,p, SOCi ≤ SOCi ≤ SOC i

where SOCi represents the SOC of battery i.

Consider the battery charging/discharging dynamics:

SOC(k+1)

i = SOC(k)i + x

(k)

i,p ∆t

where ∆t is the time step in one iteration. The energy limit

can be rewritten as power constraints:

SOCi − SOC(k+1)

i

∆t

≤ x

(k)

i,p ≤

SOC i − SOC(k+1)

i

∆t

.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 5
Letting Xbt

i = max(Xi,p , (SOCi − SOC(k+1)

i

)/∆t) and

X bt

i = min(X bti,p , SOC i − SOC(k+1)

i

/∆t), we obtain

following constraints for battery i:

Xbt

i ≤ xi,p ≤ X bt

i

, x

2

i,p + x

2

i,q ≤ (S

bt

i

)

2

(30)

where S

bt

i

is the upper bound of battery apparent power.

3) PV: Driven by time-varying solar radiation, PV systems have time-varying limits. Let P

pv

i

and S

pv

i

denote the

available active power and apparent power, respectively. The

following constraints are imposed on PV i:

0 ≤ xi,p ≤ P

pv

i

, x

2

i,p + x

2

i,q ≤ (S

pv

i

)

2

. (31)

D. Distributed Algorithm Implementation

The OPF associated with cost function (27) and constraints

(28), (30), (31) can be reformulated as a saddle-point problem in the form of (11), which can then be solved using the

model-free primal-dual algorithm proposed in Section III.

In fact, the Lagrangian (10) can be decomposed with

respect to each DER. The cost function (27b) and constraints

(30), (31) are concerned with only local DER information. They have an explicit form of gradient with respect

to each local DER control/dual variables. In addition, the

measurement-based global information f0(yb+) and f0(yb−)

that is associated with (27a), as well as the constraint

function g(yb) that is associated with (28), can be broadcast

to each DER in real time. Therefore, each DER has necessary

information to locally update its control and dual variables.

This gather-and-broadcast control architecture was previously considered in, e.g., [4], in the model-based context.

One practical concern in the distributed implementation is

the asynchronicity in the control signals and measurements.

At each time step, it is possible that DERs do not implement

their control signals at the same time, and the measurements

might not capture the system response after all local control

implementations. It is particularly true for fast communications and control. In fact, the asynchronicity can be modeled

as an additional noise source, and it can be analyzed similarly

to, e.g., [28]. In the next section, we numerically evaluate the

sensitivity of our approach to different levels of noise.

VI. NUMERICAL STUDY

A. Simulation Setup

Numerical experiments were conducted on a single-phase

variant of the IEEE 37-node test feeder [29], with DERs

connected at different locations; see the full version [26] for

the detailed configuration.

In simulation, the load and PV data were obtained from the

real data in California [30]. They were smoothed and scaled

to an appropriate magnitude for the considered distribution

system. The total demand and available PV generation in a

typical day are given in Fig. 1.

The described distribution system was simulated using

MATPOWER [31], from where the feeder head active power

and node voltages are obtained as system outputs. The output

0 6 12 18 24

0

20

40

60

80

MW

Toal loads

0 6 12 18 24

Time / Hour

0

5

10

15

MW

Available PV power

Total

Node 705

Fig. 1. Total load and available PV generation in one day, and an example

of available PV generation at one node.

measurements were obtained by adding Gaussian noises to

the simulated outputs; in particular,

yb

(k)

i = y

(k)

i + Wy

(k)

i

(32)

where the random variable W follows the normal distribution

N (0, σ2). In simulations, unless otherwise stated, σ = 0.001,

which is consistent with the standard phasor measurement

units noise values [32], [33].

For each local controller, the sine signal was chosen as

the perturbation signal according to (16). The algorithm

parameters were configured as described in [26].

B. Simulation Results

The results described here are based on the 24-hour

simulation with the load and PV data given in Fig. 1. It

is observed from Fig. 2 that the feeder head power nicely

tracks the desired reference power and the DERs can respond

quickly to unforeseen changes of the reference signal. The

subplots on the bottom in Fig. 2 show the detailed response

when the reference signal was changed. The power tracking

is less accurate during [15, 21] hours because the power

tracking was conflicting with another objective of the voltage

regulation, which is shown next.

The performance of voltage regulation is illustrated in Fig.

3, which provides the comparison of node voltages with and

0 6 12 18 24

Time / Hour

0

20

40

60

MW

Feeder head

Reference

5 5.1

0

10

20

30

10 10.1

0

10

20

30

15 15.1

20

30

40

50

21 21.1

20

30

40

50

Fig. 2. Performance of feeder head power tracking achieved by the

gradient-free approach.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 6
Fig. 3. Improvement of voltage profile by the proposed model-free OPF algorithm.

without voltage regulation. The voltage lower limit was set

to 0.96 p.u. At around 15 hours, when some node voltages

dropped below the limit, they quickly recovered to normal

values. In contrast, without voltage regulation (no voltage

limit), some node voltages stayed below 0.96 p.u. for the

next 6 hours.

These grid services were provided through the control

of battery and PV systems, whose behaviors are presented

in Fig. 4 and Fig. 5, respectively. The feeder head power

was controlled using the active power of DERs, mainly

from batteries. The PV systems have the main objective of

minimizing curtailment, so the active PV power was merely

increased because it was already near maximal generation.

But the active PV power can be largely curtailed when the

system requests less power generation (e.g., at 15 hours).

Because the active power of DERs is determined to track the

reference feeder head power, the voltage regulation is mainly

accomplished by the reactive power of DERs. At around 15

hours, when node voltages dropped below the lower limit,

-10

0

10

MW

Real power

Battery 1 Battery 2

-10

0

10

Mvar

Reactive power

0

5

10

MVA

Apparent power

0 6 12 18 24

Time / Hour

0

10

20

30

MWh

SOC

Fig. 4. Control results of two battery systems. The dashed lines represent

hard physical constraints imposed on batteries. Corresponding to each

hard constraint, the soft constraint is slightly more conservative than the

hard constraint, shown as the dotted lines. Soft constraints enable smooth

transition of control signals to avoid hard-constraint projections.

the reactive power from both batteries and PV acted quickly

to recover these voltages to the normal range.

To study the impact of the measurement noise, we tested

the proposed algorithm with different levels of noise, distinguished by the standard deviation of the noise signal in

(32). Two metrics are used to quantify system performance.

The normalized root mean squared error (NRMSE) is used

to define power tracking performance:

NRMSE =

vuut

1

K

X

K

k=1 

P

(k)

0 − P

•(k)

0

P

•(k)

0

!2

.

The average voltage violation (AVV) is used to define voltage

regulation performance:

AVV =

1

NK

X

N

i=1

X

K

k=1



[v

(k)

i − V i

]+ + [V i − v

(k)

i

]+



where N and K are the number of nodes and simulation time

steps, respectively; [ ]+ projects a negative value to zero.

Fig. 6 shows the simulation results of NRMSE and AVV.

As the noise signal increases and crosses the point of σ =

1.6 × 10−3, the control performance degrades significantly.

It is because large measurement noise can suppress the exploration output and leads to inaccurate gradient estimation.

One idea to address this issue is to increase the exploration

signal to a reasonable magnitude. But the improvement space

is limited because of practical system considerations. In that

case, we can apply state estimation techniques to filter out

measurement noises. This is beyond the scope of this paper

and left to our future work.

0

5

10

15

MW

Real power

Aggregate

Node 705

0 6 12 18 24

Time / Hour

-2

0

2

4

Mvar

Reactive power

Fig. 5. Results of PV power control.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 7
0 0.5 1 1.5 2 2.5

Standard deviation of noise, σ ×10-3

0

0.1

0.2

0.3

NRMSE

0

0.5

1

1.5

AVV (

× 10-3 p.u.)

Power tracking error

Voltage violation

Fig. 6. The performance is impacted by the measurement noise.

VII. CONCLUSIONS AND FUTURE WORK

In this paper, we developed a model-free primal-dual

method to track the desired trajectories of networked systems. The algorithm leverages the zero-order deterministic

estimation of the gradient with two function evaluations.

We provided preliminary stability and tracking results, and

illustrated an application of the method to real-time OPF in

power systems.

Analysis of the exact iteration (18)–(19), with the projection performed at every time step, remains an interesting research topic. The promising direction is to analyze the timevarying projected primal-dual dynamics in continuous time.

Also, analysis of the asynchronous algorithm is important

from a practical perspective. Finally, on the application side,

it would be interesting to show the performance of the algorithm on other scenarios, including topology reconfigurations

(e.g., caused by natural disasters and attacks on the grid).

REFERENCES

[1] S. Rahili and W. Ren, “Distributed continuous-time convex optimization with time-varying cost functions,” IEEE Transactions on

Automatic Control, vol. 62, no. 4, pp. 1590–1605, April 2017.

[2] M. Fazlyab, C. Nowzari, G. J. Pappas, A. Ribeiro, and V. M. Preciado,

“Self-Triggered Time-Varying Convex Optimization,” in Proceedings

of the 55th IEEE Conference on Decision and Control, Las Vegas,

NV, US, December 2016, pp. 3090 – 3097.

[3] A. Bernstein, E. Dall’Anese, and A. Simonetto, “Online primaldual methods with measurement feedback for time-varying convex

optimization,” IEEE Transactions on Signal Processing, vol. 67, no. 8,

pp. 1978–1991, April 2019.

[4] E. Dall’Anese and A. Simonetto, “Optimal power flow pursuit,” IEEE

Transactions on Smart Grid, vol. 9, no. 2, pp. 942–952, March 2018.

[5] A. Hauswirth, I. Subotic, S. Bolognani, G. Hug, and F. D  ́ orfler, “Time-  ̈

varying projected dynamical systems with applications to feedback

optimization of power systems,” in 2018 IEEE Conference on Decision

and Control (CDC), Dec 2018, pp. 3258–3263.

[6] M. Colombino, E. Dall’Anese, and A. Bernstein, “Online optimization

as a feedback controller: Stability and tracking,” IEEE Transactions

on Control of Network Systems, pp. 1–1, 2019.

[7] A. Bernstein, Y. Chen, M. Colombino, E. Dall’Anese, P. Mehta,

and S. Meyn, “Optimal rate of convergence for quasi-stochastic

approximation,” To appear, IEEE CDC and arXiv:1903.07228, Mar

2019.

[8] J. Kiefer and J. Wolfowitz, “Stochastic estimation of the maximum

of a regression function,” Ann. Math. Statist., vol. 23, no. 3, pp.

462–466, 09 1952. [Online]. Available: https://doi.org/10.1214/aoms/

1177729392

[9] J. C. Spall, “Multivariate stochastic approximation using a simultaneous perturbation gradient approximation,” IEEE Transactions on

Automatic Control, vol. 37, no. 3, pp. 332–341, March 1992.

[10] S. Bhatnagar, M. C. Fu, S. I. Marcus, and I.-J. Wang, “Two-timescale

simultaneous perturbation stochastic approximation using deterministic perturbation sequences,” ACM Trans. Model. Comput. Simul.,

vol. 13, no. 2, pp. 180–209, 2003.

[11] L. Prashanth, S. Bhatnagar, N. Bhavsar, M. C. Fu, and S. Marcus,

“Random directions stochastic approximation with deterministic perturbations,” IEEE Transactions on Automatic Control, 2019.

[12] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono,

“Optimal rates for zero-order convex optimization: The power of

two function evaluations,” IEEE Transactions on Information Theory,

vol. 61, no. 5, pp. 2788–2806, May 2015.

[13] Y. Nesterov and V. Spokoiny, “Random gradient-free minimization

of convex functions,” Foundations of Computational Mathematics,

vol. 17, no. 2, pp. 527–566, Apr 2017. [Online]. Available:

https://doi.org/10.1007/s10208-015-9296-2

[14] X. Gao, B. Jiang, and S. Zhang, “On the information-adaptive

variants of the admm: An iteration complexity perspective,” Journal

of Scientific Computing, vol. 76, no. 1, pp. 327–363, Jul 2018.

[Online]. Available: https://doi.org/10.1007/s10915-017-0621-6

[15] D. Hajinezhad, M. Hong, and A. Garcia, “Zone: Zeroth order nonconvex multi-agent optimization over networks,” IEEE Transactions

on Automatic Control, pp. 1–1, 2019.

[16] A. D. Flaxman, A. T. Kalai, and H. B. McMahan, “Online convex

optimization in the bandit setting: Gradient descent without a gradient,” in Proceedings of the Sixteenth Annual ACM-SIAM Symposium

on Discrete Algorithms, ser. SODA ’05. Philadelphia, PA, USA:

Society for Industrial and Applied Mathematics, 2005, pp. 385–394.

[17] B. Awerbuch and R. Kleinberg, “Online linear optimization and

adaptive routing,” Journal of Computer and System Sciences, vol. 74,

no. 1, pp. 97 – 114, 2008, learning Theory 2004.

[18] S. Bubeck and N. Cesa-Bianchi, “Regret analysis of stochastic

and nonstochastic multi-armed bandit problems,” Foundations and

Trends
R in Machine Learning, vol. 5, no. 1, pp. 1–122, 2012.

[19] T. Chen and G. B. Giannakis, “Bandit convex optimization for scalable

and dynamic iot management,” IEEE Internet of Things Journal, vol. 6,

no. 1, pp. 1276–1286, Feb 2019.

[20] K. B. Ariyur and M. Krstic, ́ Real Time Optimization by Extremum

Seeking Control. New York, NY: John Wiley & Sons, Inc., 2003.

[21] M. Krstic and H.-H. Wang, “Stability of extremum seeking feedback  ́

for general nonlinear dynamic systems,” Automatica, vol. 36, no. 4,

pp. 595 – 601, 2000.

[22] H.-H. Wang and M. Krstic, “Extremum seeking for limit cycle mini-  ́

mization,” IEEE Transactions on Automatic Control, vol. 45, no. 12,

pp. 2432–2436, Dec 2000.

[23] J. Koshal, A. Nedic, and U. Y. Shanbhag, “Multiuser optimization:  ́

Distributed algorithms and error analysis,” SIAM J. on Optimization,

vol. 21, no. 3, pp. 1046–1081, 2011.

[24] D. Shirodkar and S. Meyn, “Quasi stochastic approximation,” in

Proceedings of the 2011 American Control Conference, June 2011,

pp. 2429–2435.

[25] K. Ariyur and M. Krstic, Real-Time Optimization by ExtremumSeeking Control, ser. Wiley-interscience publication. Wiley, 2003.

[26] Y. Chen, A. Bernstein, A. Devraj, and S. Meyn, “Model-free primaldual methods for network optimization with application to real-time

optimal power flow,” arxiv, 2019.

[27] A. Bernstein, C. Wang, E. Dall’Anese, J. Le Boudec, and C. Zhao,

“Load flow in multiphase distribution networks: Existence, uniqueness,

non-singularity and linear models,” IEEE Transactions on Power

Systems, vol. 33, no. 6, pp. 5832–5843, Nov 2018.

[28] A. Bernstein and E. Dall’Anese, “Asynchronous and distributed

tracking of time-varying fixed points,” in 2018 IEEE Conference on

Decision and Control (CDC), Dec 2018, pp. 3236–3243.

[29] W. H. Kersting, Distribution system modeling and analysis. CRC

press, 2006.

[30] J. Bank and J. Hambrick, “Development of a high resolution, real

time, distribution-level metering system and associated visualization,

modeling, and data analysis functions,” May 2013.

[31] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power

systems research and education,” IEEE Transactions on Power Systems, vol. 26, no. 1, pp. 12–19, Feb 2011.

[32] L. Xie, Y. Chen, and P. R. Kumar, “Dimensionality reduction of

synchrophasor data for early event detection: Linearized analysis,”

IEEE Transactions on Power Systems, vol. 29, no. 6, pp. 2784–2794,

Nov 2014.

[33] M. Brown, M. Biswal, S. Brahma, S. J. Ranade, and H. Cao,

“Characterizing and quantifying noise in pmu data,” in 2016 IEEE

Power and Energy Society General Meeting (PESGM), July 2016, pp.

1–5.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 8