An Incentive-Based Online Optimization Framework for Distribution Grids

TL;DR

This paper introduces an online distributed optimization framework for distribution grids with DERs, enabling real-time voltage regulation and economic optimization amid dynamic conditions.

Contribution

It develops a novel online algorithm that tracks time-varying solutions for social-welfare maximization in distribution grids with real-time incentives.

Findings

Algorithm effectively tracks solutions under changing conditions.

Voltage constraints are maintained within prescribed limits.

The approach is both theoretically stable and practically effective.

Abstract

This paper formulates a time-varying social-welfare maximization problem for distribution grids with distributed energy resources (DERs) and develops online distributed algorithms to identify (and track) its solutions. In the considered setting, network operator and DER-owners pursue given operational and economic objectives, while concurrently ensuring that voltages are within prescribed limits. The proposed algorithm affords an online implementation to enable tracking of the solutions in the presence of time-varying operational conditions and changing optimization objectives. It involves a strategy where the network operator collects voltage measurements throughout the feeder to build incentive signals for the DER-owners in real time; DERs then adjust the generated/consumed powers in order to avoid the violation of the voltage constraints while maximizing given objectives. The stability of the proposed schemes is analytically established and numerically corroborated.

Tables

Table 1. TABLE I: Notation.

$𝒩$	Set of nodes, excluding node $0$; $𝒩:=\{1,\dots ,N\}$
$ℰ$	Set of distribution lines
$p_i$	Net real power injected at node $i$
$q_i$	Net reactive power injected at node $i$
$z_i$	Overall power injected at node $i$, $z_i:={[p_i,q_i]}^{⊺}$
${𝒵}_i$	Feasible set of real and reactive power at node i
$p$	$p:={[p_1,\dots ,p_N]}^{⊺}$
$q$	$q:={[q_1,\dots ,q_N]}^{⊺}$
$z$	$z:={[p^{⊺},q^{⊺}]}^{⊺}$
$V_i$	Complex voltage at node $i$
$v_i$	Voltage magnitude at node $i$
$v$	$v:={[v_1,\dots ,v_N]}^{⊺}$
${\alpha }_i$	Signal for injected real power for node $i$
${\beta }_i$	Signal for injected reactive power for node $i$
$s_i$	Overall signal $z_i:={[{\alpha }_i,{\beta }_i]}^{⊺}$
$\alpha$	$\alpha :={[{\alpha }_1,\dots ,{\alpha }_N]}^{⊺}$
$\beta$	$\alpha :={[{\beta }_1,\dots ,{\beta }_N]}^{⊺}$
$s$	Compact signal vector $s:={[{\alpha }^{⊺},{\beta }^{⊺}]}^{⊺}$
${\left[x\right]}_{+}$	Projection of $x$ onto the nonnegative orthant
${\left[x\right]}_{𝒵}$	Projection of $x$ onto the convex set $𝒵$

Equations

$${\mathcal{Z}}_i^t=\left\{(p_i,q_i):0\le p_i\le p_{i,\text{av}}^t,p_i^2+q_i^2\le {\eta }_i^2\right\}.$$

$${\mathcal{Z}}_i^t=\left\{(p_i,q_i):{\underline{p}}_i^t\le p_i\le {\overline{p}}_i^t,p_i^2+q_i^2\le {\eta }_i^2\right\},$$

$${\mathcal{Z}}_i^t=\left\{(p_i,q_i):{\underline{p}}_i^t\le p_i\le {\overline{p}}_i^t,q_i=0\right\},$$

$$P_{ij}^t$$

$$Q_{ij}^t$$

$${v_j^t}^2$$

$${\ell }_{ij}^t{v_i^t}^2$$

$${\hat{v}}^t$$

$$({\mathcal{P}}_{1,\mathbf{i}}^{\mathbf{t}})$$

$$s.t.$$

$$f_i^t(p_i^t,q_i^t|{\alpha }_i^t,{\beta }_i^t):=C_i^t(p_i^t,q_i^t)-{\alpha }_i^t{p}_i^t-{\beta }_i^t{q}_i^t$$

$$\displaystyle(p^{t*}_{i},q^{t*}_{i})=b^{t}_{i}(\alpha^{t}_{i},\beta^{t}_{i}):=\underset{(p^{t}_{i},q^{t}_{i})\in\mathcal{Z}^{t}_{i}}{\arg\min}\ \ f^{t}_{i}\big{(}p^{t}_{i},q_{i}^{t}|\alpha_{i}^{t},\beta_{i}^{t}\big{)}.$$

$$({\mathcal{P}}_2^{\mathbf{t}})$$

$$s.t.$$

$$\displaystyle p^{t}_{i}=b^{t}_{i}(\alpha^{t}_{i})=\left\{\begin{array}[]{ll}\underline{p}_{i},&\text{if}~{}\alpha^{t}_{i}/2<\underline{p}_{i}\\ \alpha^{t}_{i}/2,&\text{if}~{}\underline{p}_{i}\leq\alpha^{t}_{i}/2\leq\overline{p}_{i}\\ \overline{p}_{i},&\text{if}~{}\alpha^{t}_{i}/2>\overline{p}_{i}\end{array}\right.,$$

$$({\mathcal{P}}_3)$$

$$s.t.$$

$$\underline{v}\le R\tilde{p}+X\tilde{q}+a\le \overline{v}.$$

$${\alpha }^{\ast }$$

$${\beta }^{\ast }$$

$$\displaystyle\Big{(}\nabla_{p}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})+\gamma R\nabla_{\hat{v}}D(\hat{v}^{*})-R(\underline{\mu}^{*}-\overline{\mu}^{*})\Big{)}^{\intercal}$$

$$(p-p^{\ast })\ge 0,\forall (p,q)\in \mathcal{Z},$$

$$\displaystyle\Big{(}\nabla_{q}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})+\gamma X\nabla_{\hat{v}}D(\hat{v}^{*})-X(\underline{\mu}^{*}-\overline{\mu}^{*})\Big{)}^{\intercal}$$

$$(q-q^{\ast })\ge 0,\forall (p,q)\in \mathcal{Z},$$

$${\hat{v}}^{\ast }=R{p}^{\ast }+X{q}^{\ast }+a,$$

$$\underline{v}\le {\hat{v}}^{\ast }\le \overline{v},$$

$$({\hat{v}}^{\ast }-\underline{v}{)}^{⊺}{\underline{\mu }}^{\ast }=0,{\underline{\mu }}^{\ast }\ge 0,$$

$$(\overline{v}-{\hat{v}}^{\ast }{)}^{⊺}{\overline{\mu }}^{\ast }=0,{\overline{\mu }}^{\ast }\ge 0.$$

$$\displaystyle\Big{(}\nabla_{p}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})+\gamma R\nabla_{\hat{v}}D(\hat{v}^{*})\Big{)}^{\intercal}(p-p^{*})$$

$$\displaystyle+\Big{(}\nabla_{q}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})+\gamma X\nabla_{\hat{v}}D(\hat{v}^{*})\Big{)}^{\intercal}(q-q^{*})$$

$$+({\overline{\mu }}^{\ast }-{\underline{\mu }}^{\ast }{)}^{⊺}(\hat{v}-{\hat{v}}^{\ast })\ge 0,\forall (p,q)\in \mathcal{Z},\forall \hat{v},$$

$$\displaystyle\Big{(}\nabla_{p}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})-\alpha^{*}\Big{)}^{\intercal}(p-p^{*})\geq 0,\ \forall(p,q)\in\mathcal{Z},$$

$$\displaystyle\Big{(}\nabla_{q}\sum_{i\in\mathcal{N}}C_{i}(p_{i}^{*},q_{i}^{*})-\beta^{*}\Big{)}^{\intercal}(q-q^{*})\geq 0,\ \forall(p,q)\in\mathcal{Z}.$$

$$\mathcal{L}(z,\mu )$$

$${\mathcal{L}}_{\varphi }(z,\mu )$$

$$\underset{\mu \in {\mathbb{R}}_{+}^{2N}}{\max }\underset{z\in \mathcal{Z}}{\min }$$

$$\displaystyle\big{(}\nabla_{z}f(z,\mu)-\nabla_{z}f(z^{\prime},\mu)\big{)}^{\intercal}(z-z^{\prime})\geq c\|z-z^{\prime}\|^{2}.$$

$$\displaystyle\|z_{\phi}^{*}-z^{*}\|^{2}\leq\frac{\phi}{2{c}}\big{(}\|\mu^{*}\|^{2}-\|\mu_{\phi}^{*}\|^{2}\big{)}.$$

$${\mathcal{L}}_{\varphi }(z_{\varphi }^{\ast },\mu )\le {\mathcal{L}}_{\varphi }(z_{\varphi }^{\ast },{\mu }_{\varphi }^{\ast })\le {\mathcal{L}}_{\varphi }(z,{\mu }_{\varphi }^{\ast }),\forall z,\mu .$$

Full text

I Introduction

Market-based algorithms have been recently developed to control distributed energy assets with the objective of incentivizing end-customers to provide services to the grid while maximizing economic benefits and performance objectives [1, 2, 3]. For example, end-customers may be incentivized to adjust the output powers of distributed energy resources (DERs) in real time to aid voltage regulation [4], control the aggregate network demand [5], and follow regulating signals [6].

This paper aims to design incentive-based distributed algorithms that allow network operator and end-customers to pursue given operational and economic objectives and, in doing so, ensure that voltage magnitudes are within the prescribed limits. We start with the formulation of a time-varying social welfare maximization problem that captures a variety of optimization objectives, hardware constraints, and the nonlinear power-flow equations governing the physics of distribution systems. The time-varying nature of the problem [7, 8] enables us to model optimization and operational objectives that vary in time and to capture variability of ambient conditions and non-controllable energy assets; the time-varying problem thus defines optimal trajectories for the active and reactive powers of the DERs as well as voltage levels. A linear approximation of the nonlinear power-flow equations [9, 10, 11, 12] is utilized to facilitate the development of a computationally-tractable algorithms. Even when linear power-flow models are adopted, the resultant problem is non-convex; however, we propose a convex relaxation and provide conditions under which the relaxation is exact, i.e., the optimal solutions of the relaxed problem coincide with the global optimal points of the non-convex social-welfare problem. We then design distributed algorithms to identify (and track) the solutions of the time-varying social welfare maximization problem. The algorithm enables a distributed solution of the social-welfare-maximization problem where: (i) customers do not share private information such as their cost function and the feasible set of the DERs’ output powers; and (ii) customers and network operator pursue their own economic and operational objectives, while ensuring that voltage limits are systematically satisfied throughout the network.

The first algorithm is applicable to problems that vary slowly in time, where offline iterative methods can be utilized to solve sampled instances [7] of the time-varying problem to convergence. An online algorithm is then proposed to enable tracking of the solutions in the presence of fast time-varying operational conditions and changing optimization objectives. The online algorithm involves a strategy where the network operator collects voltage measurements through the feeder to build incentive signals for the DER-owners in real time; DERs then adjust the generated/consumed powers in order to avoid the violation of the voltage constraints while maximizing given objectives. Stability of the proposed schemes is analytically established and numerically corroborated. The design of the algorithms is grounded on the decomposability of the primal-dual gradient algorithms, and convergence of the distributed algorithm to the optimal solution of the social-welfare maximization problem is shown by leveraging the contraction mapping arguments with a regularized Lagrangian.

It should be pointed out that, traditionally, voltage regulation problems arising from reverse power flows [4] have been tackled by considering local Volt/VAR and Volt/Watt controllers [13, 14, 15, 16, 17, 18, 19] or optimization-based techniques [11, 20, 21] that leverage the flexibility of power-electronics-interfaced renewable sources of energy in adjusting the output real and reactive powers. Compared to local control strategies, the proposed method allows the network operator and end-customers to pursue well-defined performance objectives; compared to existing optimization strategies, the proposed method casts the voltage regulation problem within a pricing/incentive realm, and provides insights as to how to design real-time pricing/incentive schemes.

Load control is formulated as a Stackelberg game in e.g., [22, 23], and [5], and conditions for convergence of a leader-follower strategy to an equilibrium point are derived. Based on a Stackelberg game formulation, a distributed algorithm with only local information available for both utility companies and DER-owners is developed in [2]. A two-level game setting is considered in [24]. A time-varying pricing strategy is proposed in [25] to maximize social welfare. However, the strategies outlined in [2, 22, 23, 5, 24, 25] and pertinent references therein are network agnostic, in the sense that AC power flows in the power network are ignored (and DERs are assumed to be connected to one single electrical node). It follows that network-agnostic methods do not account for voltage variations induced by the controlled DERs and must be complemented by voltage-regulation mechanisms. The frameworks proposed in, e.g., [25, 26] offer a way to account for the power flows, but their applicability is limited to a restricted class of network topologies.

As for online optimization methods for distribution systems, Centralized controllers are developed in [27, 28], based on continuous gradient steering algorithms; the framework accounts for errors in the implementable power setpoints, and convergence of the average setpoints to the minimum of the considered control objective is established. An online gradient algorithm for AC optimal power flow in single-phase radial networks is proposed in [29]; it is shown that the proposed algorithm converges to the set of local optima of a static AC OPF problem, and sufficient conditions under which the online OPF converges to a global optimum are provided. A real-time control strategy that enables DERs to maximize given performance objectives is proposed in [8]. The proposed online algorithm is close in spirit to [8]; however compared to [8], this paper casts the real-time voltage regulation problem within a time-varying game-theoretic framework and develops distributed strategies based on pricing/incentive signals; [8] does not address the design of pricing/incentive signals.

The paper is organized as follows. Section II introduces the system model and presents the problem formulation. Section III focuses on the design of iterative algorithms that afford an offline implementation, while Section IV presents the online algorithm. Section V outlines results from numerical experiments and Section VI concludes the paper. Preliminary results were presented in [39].

II Preliminaries and System Model

II-A Network Model

Consider a distribution network with $N+1$ nodes collected in the set ${\cal N}\cup\{0\}$ with ${\cal N}:=\{1,...,N\}$ and node [math] being the point of common coupling or substation, and distribution lines collected in the set ${\cal E}$. Let $V^{t}_{i}\in\mathbb{C}$ denote the line-to-ground voltage at node $i$ at time $t$, and define $v^{t}_{i}:=|V^{t}_{i}|$. Denote as $p^{t}_{i}\in\mathbb{R}$ and $q^{t}_{i}\in\mathbb{R}$ the (net) active and reactive power injections, respectively, of a distributed energy resource (DER) located at node $i\in{\cal N}$. For notational simplicity, exposition is tailored to the case where one DER is located at each node; however, the technical approach straightforwardly applies to the case where multiple DERs are connected to a node. Hereafter, $\mathcal{Z}^{t}_{i}$ denotes the feasible set of active and reactive powers $p^{t}_{i}$ and $q^{t}_{i}$ at node $i\in{\cal N}$ at time $t$. In the following, we explain how to construct this set for some types of DERs.

Photovoltaic (PV) systems: Let $p_{i,\text{av}}^{t}$ denote the available real power from a PV system at time $t$, and let $\eta_{i}$ be the rated apparent capacity. Then, the set ${\cal Z}_{i}^{t}$ is given by:

[TABLE]

Energy storage systems: The set ${\cal Z}_{i}^{t}$ for an energy storage system is given by:

[TABLE]

for given limits $\underline{p}_{i}^{t},\overline{p}_{i}^{t}$ and for a given inverter rating $\eta_{i}$. The limits $\underline{p}_{i}^{t},\overline{p}_{i}^{t}$ are updated during the operation of the battery based on the state of charge.

Variable frequency drives: For devices such as water pumps and supply fans of commercial HVAC systems, the set ${\cal Z}_{i}^{t}$ can be described as:

[TABLE]

for given limits $\underline{p}_{i}^{t},\overline{p}_{i}^{t}$. These limits can be fixed or updated by local controllers at a regular time intervals, based on the state of e.g., thermal loads.

The operating region of small-scale diesel generators can be modeled using constant box constraints. For DERs with discrete levels of output powers (e.g., electric vehicle chargers with discrete charging levels), ${\cal Z}_{i}$ represents the convex envelope of the possible operating points; see e.g., [28]. Randomization techniques can then be utilized to recover a feasible setpoint. However, the development of control strategies for DERs with discrete levels of output powers is left as a future research activity.

Voltages, currents, and powers $\{p^{t}_{i},q^{t}_{i}\}$ are related by the well-known nonlinear AC power-flow equations; assuming, for illustrative purpose, a balanced tree network, these equations read:

[TABLE]

where $\ell^{t}_{ij}$ is the squared magnitude of the current on line $(i,j)$, $P^{t}_{ij},Q^{t}_{ij}$ are real and reactive powers injected on line $(i,j)$, and $r_{ij}+\mathrm{j}x_{ij}$ is the impedance on line $(i,j)$.

To facilitate the design and analysis of computationally-tractable algorithms, the proposed approach will employ suitable linearization approaches for (1). Particularly, the following approximate linear relationship between voltage magnitudes and injected powers is utilized:

[TABLE]

where the parameters $R,X\!\!\!\in\mathbb{R}_{++}^{N\!\times\!N}$ and $a\!\!\in\mathbb{R}^{N}$ can be obtained using one of the two following approaches: i) regression-based methods, based on real-time measurements of $\{v_{i}^{t}\}$, $p^{t}$, and $q^{t}$, e.g., the recursive least-squares method [30] can be utilized to continuously update the model parameters; and, ii) suitable linearization methods for the AC power-flow equations; see e.g., [9, 13, 10, 11, 12]. In the latter case, the model parameters $R$, $X$, and $a$ can be time-varying too, by using current operating points as linearization points for the AC power-flow equations. Parameters $R$, $X$, and $a$ should be re-computed every time that the system changes topology.

The approximate model (2) is utilized to facilitate the design of computationally-affordable algorithms. Section IV will show how to leverage appropriate measurements to cope with approximation errors and systematically enforce voltage limits.

Remark 1**.**

*(Multiphase systems) For notational and exposition simplicity, the framework is outlined for a single-phase system. However, the proposed algorithmic solution is applicable to unbalanced multiphase networks. This can be obtained by substituting (2) with the linearized model recently proposed in [31] for unbalanced multiphase networks with both wye-connected and delta-connected DERs. $\Box$ *

II-B Problem Setup

The goal is to design a strategy wherein the network operator and end-customers pursue their own operational and economic objectives, while achieving a global coordination to enforce voltage regulation.

II-B1 End-customer optimization problem

Consider a cost function $C^{t}_{i}(p^{t}_{i},q^{t}_{i})$ that captures a well-defined performance objective for the customer(s) located at node $i\in{\cal N}$ at time $t$. Let $\alpha^{t}_{i}\in\mathbb{R}$ and $\beta^{t}_{i}\in\mathbb{R}$ be incentive signals produced by the network operator (e.g., distribution system operator or aggregator) for active and reactive power injections, respectively, at time $t$. Given signals $(\alpha_{i}^{t},\beta^{t}_{i})$, the following optimization problem is solved at each node $i\in{\cal N}$ at time $t$:

[TABLE]

where

[TABLE]

with $\alpha^{t}_{i}p^{t}_{i}$ and $\beta^{t}_{i}q^{t}_{i}$ representing payment to/reward from the network operator. The following standard assumption is made.

Assumption 1**.**

*Functions $C^{t}_{i}(p^{t}_{i},q^{t}_{i}),\ \forall i\in\mathcal{N}$ are continuously differentiable and strongly convex in $(p^{t}_{i},q^{t}_{i})$. Moreover, the first-order derivative of $C^{t}_{i}(p^{t}_{i},q^{t}_{i})$ is bounded in $\mathcal{Z}_{i}$. *

The assumption of bounded derivative means that an infinitesimal change in power should not lead to a jump in cost. Because (3a) is strictly convex in $(p_{i}^{t},q_{i}^{t})$ and $\mathcal{Z}^{t}_{i}$ is convex and compact, a unique solution $(p_{i}^{t*},q_{i}^{t*})$ exists for each $t$.

For future developments, consider the so-called best response strategy of node $i$, denoted as $b^{t}_{i}(\alpha^{t}_{i},\beta^{t}_{i})$, for given $\alpha^{t}_{i}$ and $\beta^{t}_{i}$:

[TABLE]

II-B2 Social-welfare problem

Consider a cost function $D^{t}(\hat{v}^{t})$ that captures network-oriented objective in voltage at time $t$. For example, to minimize the voltage deviation from the nominal value $v^{\text{nom}}$, we can set $D^{t}(\hat{v}^{t})=\frac{1}{2}\|\hat{v}^{t}-v^{\text{nom}}\|^{2}$. The following assumption is made.

Assumption 2**.**

*Function $D^{t}(\hat{v}^{t})$ is continuously differentiable, convex, and with bounded first-order derivative at achievable voltage magnitude values. *

Because the set of power injections $(p,q)$ is compact and $\hat{v}$ is a continuous function of $(p,q)$, the achievable $\hat{v}$ values are bounded. Thus, the boundedness of the first-order derivative of $D^{t}(\hat{v}^{t})$ is a reasonable assumption.

Consider the following optimization problem, which captures both customer-oriented and network-oriented objectives of a distribution network:

[TABLE]

where $\gamma^{t}\in\mathbb{R}_{+}$ is used to trade off between the end-customer and network-oriented objectives, and $\underline{v}^{t}$ and $\overline{v}^{t}$ are vectors collecting prescribed minimum and maximum voltage magnitude limits (the inequalities are component-wise).

Note that $(\bm{\mathcal{P}^{t}_{2}})$ is usually non-convex due to the constraint (6d). This is because (6d) is usually not affine. For better illustration of the non-convexity of $(\bm{\mathcal{P}^{t}_{2}})$, consider the following simple example with real power only. Assume a quadratic cost function $C^{t}_{i}(p_{i}^{t})={p_{i}^{t}}^{2}$ and a box feasible set $p^{t}_{i}\in[\underline{p}_{i},\overline{p}_{i}]$ with upper and lower bounds for real power injections $\underline{p}_{i}$ and $\overline{p}_{i}$, and we end up with a non-convex piece-wise linear function $b^{t}_{i}$:

[TABLE]

which will become more complex if we consider more complicated $C_{i}^{t}$ and $\mathcal{Z}_{i}^{t}$.

Problem $(\bm{\mathcal{P}^{t}_{2}})$ defines optimal operational trajectory $\{p^{t*},q^{t*}\}_{t\in\mathbb{R}_{+}}$ over time for the active and reactive powers of the distribution network. One way to identify (and track) the time-varying optimal points of $(\bm{\mathcal{P}^{t}_{2}})$ consists in discretizing the temporal domain as $t_{m}:=mh,~{}k\in\mathbb{N}$, where $h$ is a given time interval, and solve $(\bm{\mathcal{P}^{t}_{2}})$ at each time $t_{m}$. Section III will focus on the case where distributed algorithms can be utilized to solve $(\bm{\mathcal{P}^{t_{m}}_{2}})$ to convergence at each time $t_{m}$. These algorithms are suitable for operational conditions where non-controllable demand/generation, ambient conditions, and cost functions are slow time-varying (i.e., the interval $h$ is “large enough” to allow convergence of the distributed algorithm). Section IV will then focus on faster time-varying operational settings, and will advocate the development of online algorithms [7] that track $\{p^{t*},q^{t*}\}_{t\in\mathbb{R}_{+}}$ over time.

III Incentive-Based Distributed Algorithm

Focusing on a particular problem instance at time $t$, $(\bm{\mathcal{P}^{t}_{2}})$ lends itself to a Stackelberg game interpretation where $\alpha^{t}$ and $\beta^{t}$ are calculated via $(\bm{\mathcal{P}^{t}_{2}})$ by the network operator (i.e., the leader) and broadcasted to all nodes $i\in{\cal N}$; subsequently, each end-consumer (i.e., the follower) computes the power setpoints $p_{i}^{t*}$ and $q_{i}^{t*}$ from $(\bm{\mathcal{P}^{t}_{1,i}})$. By design, $(p^{t*},q^{t*})$ is an optimal point of $(\bm{\mathcal{P}^{t}_{2}})$.

However, it is challenging for the network operator to solve $(\bm{\mathcal{P}^{t}_{2}})$ not only because of the non-convexity introduced by constraint (6d), but also because it requires knowledge of the end-customer’s best-response function $b_{i}^{t}$. To solve the problem, in Section III-A we first formulate a convex relaxation of $(\bm{\mathcal{P}^{t}_{2}})$ and show that its optimum gives the optimum of $(\bm{\mathcal{P}^{t}_{2}})$, and then in Section III-B we design a distributed algorithm to solve $(\bm{\mathcal{P}^{t}_{2}})$ based on the algorithm for the relaxed problem. Since the same solution procedure is used to solve $(\bm{\mathcal{P}^{t}_{2}})$ to convergence at each time $t$, the superscript $t$ will be dropped in this section. The superscript $t$ will be re-introduced in Section IV where we outline an online solution method.

III-A Convex Reformulation

We start by deriving a convex relaxation of the non-convex problem $(\bm{\mathcal{P}_{2}})$ as well as conditions under which an optimal point of $(\bm{\mathcal{P}_{2}})$ can be identified. Consider the following convex optimization problem:

[TABLE]

where we replace the non-convex constraint (6d) in $(\bm{\mathcal{P}_{2}})$ with (11d). We assume that the above problem is feasible.

Assumption 3** (Slater’s condition).**

There exists a feasible point $(\tilde{p},\tilde{q})\in\mathcal{Z}$, $\mathcal{Z}:=\mathcal{Z}_{1}\times\ldots\times\mathcal{Z}_{N}$, such that:

[TABLE]

Assumption 3 does not involve strict inequality because the constraint is linear. Given the strong convexity of the objective function (11a) in $(p_{i},q_{i})$ and the linear relation (11d), a unique optimal solution exists for problem $(\bm{\mathcal{P}_{3}})$. Notice that a solution $(p_{i}^{*},q_{i}^{*},{\hat{v}}^{*})$ of $(\bm{\mathcal{P}_{3}})$ may not be feasible for $(\bm{\mathcal{P}_{2}})$, i.e., there does not exist a $(\alpha^{*},\beta^{*})$ such that $(p_{i}^{*},q_{i}^{*})=b_{i}(\alpha_{i}^{*},\beta_{i}^{*})$. We will, however, show next that such a $(\alpha^{*},\beta^{*})$ exists, and thus the solution of $(\bm{\mathcal{P}_{3}})$ gives the solution of $(\bm{\mathcal{P}_{2}})$.

Denote by $\underline{\mu}$ and $\overline{\mu}$ the dual variables associated with the constraint (11d). Let $\hat{v}^{*}$ be the optimal voltage magnitudes produced by $(\bm{\mathcal{P}_{3}})$ and $\underline{\mu}^{*},\overline{\mu}^{*}$ the optimal dual variables. Then, we propose to design the incentive signals as follows:

[TABLE]

where $\nabla_{\hat{v}}D$ denotes the gradient of function $D$ with respect to the vector $\hat{v}$. Note that $\alpha^{*}$ and $\beta^{*}$ are composed of dual prices $\underline{\mu}^{*},\overline{\mu}^{*}$ and the marginal cost of network operator $\gamma\nabla_{\hat{v}}D(\hat{v}^{*})$, together with $R,X$ characterizing the network structure. As will be shown shortly, $\alpha^{*}$ and $\beta^{*}$ are de facto designed based on the optimality conditions of $(\bm{\mathcal{P}_{2}})$ and $(\bm{\mathcal{P}_{3}})$. The above incentive signals are bounded, which precludes the possibility of infinitely large signals.

Theorem 1**.**

*Under Assumptions 1–3, the incentive signals $(\alpha^{*},\beta^{*})$ defined by (13) are bounded. $\Box$ *

Proof.

Notice that the derivative $\nabla_{\hat{v}}D$ is bounded. To show the boundedness of $(\alpha^{*},\beta^{*})$, it is enough to show that the optimal duals $(\underline{\mu}^{*},\overline{\mu}^{*})$ are bounded.

Consider the KKT conditions for problem $(\bm{\mathcal{P}_{3}})$:

[TABLE]

Combining (14a)–(14c) results in:

[TABLE]

where the first two terms on the left of the inequality are bounded because of the bounded derivative of cost functions and the bounded set $\mathcal{Z}$. By the complementary slackness conditions (14e)-(14f), $\overline{\mu}^{*}_{i}$ and $\underline{\mu}^{*}_{i}$, $i\in\mathcal{N}$ cannot be nonzero at the same time. If $\overline{\mu}^{*}_{i}\rightarrow\infty$, then $\hat{v}^{*}_{i}=\overline{v}_{i}$ and we can choose a $(p,q)$ and thus $\hat{v}_{i}$ such that the third term on the left of (15) goes to $-\infty$ and (15) does not hold. So, $\overline{\mu}^{*}_{i}$ and thus $\overline{\mu}^{*}$ is bounded. Similarly, we can show that $\underline{\mu}^{*}$ is bounded too. The result follows.

By examining the optimality conditions of $(\bm{\mathcal{P}_{2}})$ and $(\bm{\mathcal{P}_{3}})$, we have the following result.

Theorem 2**.**

*The solutions of problem $(\bm{\mathcal{P}_{3}})$ along with the signals $(\alpha^{*},\beta^{*})$ defined in (13) are global optimal solutions of problem $(\bm{\mathcal{P}_{2}})$; i.e., problem $(\bm{\mathcal{P}_{3}})$ is an exact convex relaxation of problem $(\bm{\mathcal{P}_{2}})$. $\Box$ *

Proof.

By the signal design (13), (14a)–(14b) become

[TABLE]

Notice that the above variational inequalities imply that $(p_{i}^{*},q_{i}^{*})=b_{i}(\alpha_{i}^{*},\beta_{i}^{*}),~{}i\in\mathcal{N}$. So, the solution of problem $(\bm{\mathcal{P}_{3}})$ along with $(\alpha^{*},\beta^{*})$ defined in (13) is feasible for problem $(\bm{\mathcal{P}_{2}})$. The result follows, as $(\bm{\mathcal{P}_{3}})$ is a convex relaxation of $(\bm{\mathcal{P}_{2}})$.

From now on, we will use the optima of $(\bm{\mathcal{P}_{3}})$ and $(\bm{\mathcal{P}_{2}})$ interchangeably depending on the context. Next, based on Theorem 2, we will develop an iterative algorithm that achieves the optimum of $(\bm{\mathcal{P}_{3}})$ (and hence that of $(\bm{\mathcal{P}_{2}})$) without exposing any private information of the end-customers to the network operator.

Remark 2**.**

*Theorem 2 asserts that non-convex problem $(\bm{\mathcal{P}_{2}})$ can be solved through solving a convex problem $(\bm{\mathcal{P}_{3}})$. At first glance, it appears that the non-convexity of $(\bm{\mathcal{P}_{2}})$ comes from a non-convex representation of the feasible set that may have a convex representation as implied by $(\bm{\mathcal{P}_{3}})$. An ongoing investigation is to identify the specific problem structure to generalize the result in Theorem 2 to a larger class of problems. $\Box$ *

III-B Distributed Algorithm

For notational simplicity, let $s_{i}=[\alpha_{i},\beta_{i}]^{\intercal}$ denote the overall signals for end-customer $i$ and define $z_{i}=[p_{i},q_{i}]^{\intercal}$. Further denote by $z:=[p^{\intercal},q^{\intercal}]^{\intercal}\in\mathbb{R}^{2N}$ the vector of stacked power injections, and by $\mu:=[\underline{\mu}^{\intercal},\overline{\mu}^{\intercal}]^{\intercal}\in\mathbb{R}_{+}^{2N}$ the vector of stacked dual variables. Consider the following Lagrangian function associated with $(\bm{\mathcal{P}_{3}})$:

[TABLE]

which is obtained by keeping the constraints $z\in\mathcal{Z}$ and $\mu\in\mathbb{R}_{+}^{2N}$ implicit. Denote as $(z^{*},\mu^{*})$ a saddle-point of $\mathcal{L}(z,\mu)$.

To facilitate the development of provably convergent online algorithms (the subject of Section IV), consider the following regularized Lagrangian function:

[TABLE]

where $\phi>0$ is a predefined parameter (see e.g., [32, 7]). With the regularization term $-\frac{\phi}{2}\|\mu\|^{2}$, the resultant function $\mathcal{L}_{\phi}(z,\mu)$ is strongly concave in the dual variables. Based on (18), we proceed with the following minimax problem:

[TABLE]

In general, the unique optimizer of (19), denoted by $(z_{\phi}^{*},\mu_{\phi}^{*})$, is not a saddle-point of the Lagrangian function (18) because of the regularization term $-\frac{\phi}{2}\|\mu\|^{2}$. However, the discrepancy between the unique optimizer of (19) and the optimizers of (18) can be bounded as shown next.

Notice first that the boundedness of $\mu^{*}$ is shown in Theorem 1; $\mu_{\phi}^{*}$ can be readily shown to be bounded too. For ease of exposition, define $f(z):=\sum_{i\in{\cal N}}C_{i}(z_{i})+\gamma D(z)$ and $g(z):=\begin{bmatrix}\underline{v}-Rp-Xq-a\\ Rp+Xq+a-\overline{v}\end{bmatrix}$; this way, the Lagrangian can be re-expressed in a compact form as $\mathcal{L}(z,\mu)=f(z)+\mu^{\intercal}g(z)$ and the regularized counterpart reads $\mathcal{L}_{\phi}(z,\mu)=f(z)+\mu^{\intercal}g(z)-\frac{\phi}{2}\|\mu\|^{2}$. From Assumption 1–2, it follows that $f$ is strongly convex in $z$. Equivalently, $\nabla_{z}f(z,\mu)$ is strongly monotone in $z$. Therefore, we have the following lemma:

Lemma 1**.**

There exists a scalar $c>0$ such that $\forall z,z^{\prime}\in\mathcal{Z}$,

[TABLE]

$\Box$* *

Then, the discrepancy between $z^{*}$ and $z_{\phi}^{*}$ due to the regularization term can be bounded as follows (see also [32, Proposition 3.1]).

Theorem 3**.**

The difference between $z_{\phi}^{*}$ and $z^{*}$ is bounded as:

[TABLE]

$\Box$* *

Proof.

As a saddle point of (19), $(z_{\phi}^{*},\mu_{\phi}^{*})$ satisfies the following inequalities:

[TABLE]

The left inequality leads to

[TABLE]

where we set $\mu=\mu^{*}$. We next characterize the term $(\mu_{\phi}^{*}-\mu^{*})^{\intercal}g(z_{\phi}^{*})$.

(i) Leveraging the definition of convex functions, $g_{j}(z_{\phi}^{*})$ can be upper bounded as:

[TABLE]

where the second inequality is due to the fact that $g_{j}(z^{*})\leq 0$. Multiply both sides of (23) by $\mu_{\phi,j}^{*}$ (which is nonnegative) and sum up for all $j$ to have:

[TABLE]

where the second inequality is due to the first-order optimality condition $\nabla_{z}\mathcal{L}_{\phi}(z_{\phi}^{*},\mu_{\phi}^{*})^{\intercal}(z_{\phi}^{*}-z^{*})\leq 0$.

(ii) On the other hand, one has that:

[TABLE]

Multiply both sides of (25) by $-\mu_{j}^{*}$ (which is nonpositive) and sum up for all $j$ to get:

[TABLE]

where the first equality is due to the complimentary slackness condition and the second inequality is obtained from the first-order optimality condition.

Substitute (24) and (26) into (22), and use (20) to obtain (21).

The key advantage of utilizing the regularized Lagrangian is that the primal-dual gradient methods applied to (19) exhibit improved convergence properties [7] as explained next.

Hereafter, we omit the subscript $\phi$ from the optimization variables for notational simplicity, with the understanding that the updates of $z(k)$ and $\mu(k)$ are designed to solve the regularized saddle-point problem (19). Consider the following primal-dual projected gradient method, where $k$ denotes the iteration index:

[TABLE]

where $[~{}]_{\mathcal{Z}\times\mathbb{R}^{2N}_{+}}$ denotes the projection operation onto the set ${\mathcal{Z}\times\mathbb{R}^{2N}_{+}}$, and $\varepsilon_{1},\varepsilon_{2}>0$ are prescribed step sizes for the primal and the dual updates. Notice that $\nabla_{z}\mathcal{L}_{\phi}(z,\mu)$ and $\nabla_{\mu}\mathcal{L}_{\phi}(z,\mu)$ are Lipschitz continuous and strongly monotone. Therefore by virtue of [33, Sec. 3.5, Proposition 5.4], the following result holds.

Theorem 4**.**

*There exist some $\bar{\varepsilon}_{1},\bar{\varepsilon}_{2}>0$ such that for any $\varepsilon_{1}\in(0,\bar{\varepsilon}_{1}],\varepsilon_{2}\in(0,\bar{\varepsilon}_{2}]$, $\hat{T}$ is a contraction mapping. For $\varepsilon_{1}\in(0,\bar{\varepsilon}_{1}],\varepsilon_{2}\in(0,\bar{\varepsilon}_{2}]$, the sequence $\{(z(k),\mu(k))\}$ generated by (27) converges geometrically to the optimizer of (19). $\Box$ *

The proof is referred to [33] and omitted here. We can further provide analytical bound for such $\bar{\varepsilon}_{1},\bar{\varepsilon}_{2}$ for completeness. The results put as Theorem 6 are presented in the Appendix for better readability. We also refer to Section V-B1) for some numerical characterization of step sizes as regards convergence.

Given Theorems 3–4, algorithm (27) converges to within a small neighborhood of problem $(\bm{\mathcal{P}_{3}})$ (problem $(\bm{\mathcal{P}_{2}})$) whose size can be controlled by choosing a proper weight $\phi$ for the regularization term.

The decomposable structure of (27) naturally enables the following iterative distributed algorithm:

[TABLE]

where the power setpoints of each device are computed locally through (28a) and (28b)–(28f) are performed at the network operator. The resultant scheme is tabulated as Algorithm 1. Notice that each end-customer $i$ does not share its cost function $C_{i}$ or its feasible set $\mathcal{Z}_{i}$ with the network operator; rather, the end-customer transmits to the network operator only the resultant power injections $z_{i}(k)$. Indeed, the results of Theorem 3–4 apply to (28) too.

Remark 3**.**

*In (28), $\alpha$ and $\beta$ are utilized by the end-customers to construct the primal gradient $\nabla_{z}\mathcal{L}_{\phi}$. This strategy enables a distributed implementation of the primal-dual projected gradient method (27) without the end-customers knowing information of the network. $\Box$ *

IV Online Algorithm

Algorithm 1 involves an iterative procedure where (28) is repeated until convergence to (within a small neighborhood of) an optimum of problem $(\bm{\mathcal{P}_{2}^{t}})$ at each time $t$. This requires considering a discrete-time model where time is divided into slots of equal duration, indexed by $t_{m}\in\mathbb{N}_{+}=\{0,1,2\cdots\}$, where the timeslot duration is expected to be much longer than the convergence time of Algorithm 1. However, in the case of fast changing operational conditions and cost functions, it is desirable to use a small timeslot duration, i.e., to sample $(\bm{\mathcal{P}_{2}^{t}})$ at a small sampling interval to track $\{{p^{t}}^{*},{q^{t}}^{*}\}_{t\in\mathbb{R}_{+}}$ [7]. In this case, it may not be possible to solve $(\bm{\mathcal{P}_{2}^{t_{m}}})$ to convergence within the timeslot, and instead only $K>0$ iterations may be performed. Also, problem $(\bm{\mathcal{P}_{2}^{t_{m}}})$ uses linear approximation (2) of the power-flow model. In this section, we will develop an online algorithm that continuously pursues the optima of $(\bm{\mathcal{P}_{2}^{t_{m}}})_{t_{m}\in\mathbb{N}_{+}}$, and characterize the “loss” of optimality caused by the finite iterations as well as the approximation error.

IV-A Online Algorithm

As in Section III-B, the voltage magnitudes $v_{i},~{}i\in\mathcal{N}$ will be measured, but they will follow the power-flow equations (1) instead of its linear approximation (2). We assume that the approximation error is bounded (see also [34, 35]).

Assumption 4**.**

*There exists a constant $e>0$ such that $|v_{i}^{t_{m}}(z)-\hat{v}_{i}^{t_{m}}(z)|\leq e,\ i\in\mathcal{N}$ for all $z\in\mathcal{Z}$ at any time $t_{m}$. *

As the linearized power-flow model (2) is a very accurate approximation under normal operating condition[9, 11, 12], the bound $e$ expects to be small.

With the measurement of the voltage magnitudes, the proposed algorithm, formally described as Algorithm 2, executes the following steps at time $t_{m}$ ($k$ denotes the iteration index):

[TABLE]

Iterations (29) are performed $K>0$ times during each timeslot $t_{m}$. When $K=1$, only one iteration (29) is computed per timeslot [7, 8, 29, 28, 27]. It is worth noticing that a centralized implementation of (29) requires collecting the time-varying $C_{i}^{t}$ and $\mathcal{Z}_{i}^{t}$ at the network operator at each iteration; hence, a a centralized implementation would incur a higher communication overhead.

In the next subsection, we will analyze the convergence and tracking capability of the above online algorithm.

Remark 4**.**

*(local controller) The proposed algorithms produce setpoints $(p_{i}^{t},q_{i}^{t})\in{\cal Z}_{i}^{t}$ for the output powers of the DERs. It is assumed that the DERs are endowed with local controllers that are designed so that, upon receiving the setpoint, the output powers are driven to the commanded setpoints. Relevant dynamical models for the output powers of inverters operating in a grid-connected mode are discussed in e.g., [36, 37] and can be found in datasheets of commercially available DERs. Assumption 4 accounts for both measurement errors and bounds the discrepancy between the commanded setpoint and the actual output powers; Assumption 4 is valid, for example, when the DER’s response to a step-change in the setpoint follows a first-order model [36, 37]. $\Box$ *

Remark 5**.**

*(implementation) The proposed real-time algorithms update the setpoints of the DERs on a second or subsecond timescale to maximize the operational objectives while coping with the variability of available renewable-based generation and non-controllable energy assets. The algorithm does not control any anti-islanding and ride-through parameters. Considerations regarding the recloser-fuse coordination problem (which affects anti-islanding and ride-through configurations) pertain to the deployment of the DERs and given interconnection agreements. In case of event where the DERs are required to shut off, the algorithm will not produce any setpoint; the algorithm will re-start producing setpoints once the DERs are allowed to reconnect to the system. $\Box$ *

IV-B Performance Analysis

In this subsection, the hatted symbols (e.g., $\hat{\alpha}^{t}(k)$, $\hat{\beta}^{t}(k)$) refer to the iterates produced by the algorithm (28) (i.e., under the linear approximation (2)), while the non-hatted symbols (e.g., $\alpha^{t}(k)$, $\beta^{t}(k)$) refer to those produced by (29) (i.e., under the nonlinear model (1)).

By comparing (28) and (29), Assumption 4 leads to the following bounds:

[TABLE]

for some $\tilde{v}^{t_{m}}$ and, therefore:

[TABLE]

Let $\delta:=[\delta_{1,1},\ldots,\delta_{1,N},\delta_{2,1},\ldots,\delta_{2,N}]\in\mathbb{R}_{+}^{2N}$, and collect the primal and dual variables in the vector $y:=(z,\mu)$ for notational simplicity. Consequently, the following holds:

[TABLE]

where $\rho:=[\varepsilon_{2}e\cdot\mathbf{1}_{1\times 2N},\delta^{\intercal}]^{\intercal}$ and $T^{t_{m}}(\cdot)$ is the counterpart of $\hat{T}^{t_{m}}(\cdot)$ for the iterates (29) at time $t_{m}$. Let $\Delta\leq\bar{\Delta}<1$ be the contraction modulus for $\hat{T}^{t_{m}}(\cdot)$; with appropriate step sizes $\varepsilon^{t_{m}}_{1}$ and $\varepsilon^{t_{m}}_{2}$ chosen according to Theorem 4, by definition, we have that:

[TABLE]

Recall that $\hat{y}^{t_{m}*}$ denotes an optimizer of $\mathcal{L}_{\phi}^{t_{m}}$. Since $\mathcal{L}_{\phi}^{t_{m}}$ is a time-varying problem, consider capturing the variation of an optimizer over two consecutive time instants as:

[TABLE]

where $0<\sigma<+\infty$ [7].

Then, the following result characterizes the discrepancy between the powers produced by (29) and an optimizer of $\mathcal{L}_{\phi}^{t_{m}}$.

Theorem 5**.**

Under Assumptions 1–4 and step sizes chosen according to Theorem 4, the sequence $\{{y}^{t_{m}}\}$ generated by Algorithm 2 converges as

[TABLE]

$\Box$* *

Proof.

We can characterize the distance between the operating point achieved by (29) in $K$ iterations and the optimizer of $\mathcal{L}^{t_{m}}_{\phi}$ as follows:

[TABLE]

where: (36) follows from (30) and (31); (40) can be obtained by repeating (35)–(36) for $K$ times; and, (40) follows from (32). We repeat steps (35)–(40) recursively over time instants $t_{m},\ldots,t_{0}$ to obtain:

[TABLE]

When $m\rightarrow\infty$, (33) follows.

The result (33) bounds the maximum discrepancy between the setpoints ${y}^{t_{m}}(K)$ generated by Algorithm 2 and a time-varying optimizer of $(\bm{\mathcal{P}^{t}_{2}})$. The bound (33) depends on:

i) The underlying dynamics of the distribution system; in fact, we recall that the dynamics of non-controllable power assets, constraints, and operational conditions translate into temporal variations of the optimizers of $(\bm{\mathcal{P}^{t}_{2}})$ [7], which is characterized by the parameter $\sigma$. When the variation of $\hat{y}^{t_{m}*}$ is smooth in time, bound (33) becomes tighter.

ii) The approximation error introduced by the linearized power-flow equation, which is implicitly captured by $\rho$.

The result (33) can also be interpreted as an input-to-state stability, when one adopts the trajectory $\{\hat{y}^{t_{m}*}\}$ as a reference frame. Future research efforts will aim at characterizing the discrepancy between ${y}^{t_{m}}(K)$ and the optimal point of $(\bm{\mathcal{P}^{t}_{2}})$ when the nonlinear AC power-flow equations are utilized.

Finally, when $K=1$, (33) boils down to:

[TABLE]

The bound in (41) is close in spirit to [8], and it can be further simplified as $\frac{\|\rho\|}{1-\Delta}$ for the case of time-invariant (i.e. static) settings [7].

V Numeric Examples

V-A Simulation Setup

Consider a modified version of the IEEE 37-node test feeder shown in Figure 1. The modified network is obtained by considering the phase “c” of the original system and by replacing the loads specified in the original dataset with real load data measured from feeders in Anatolia, California, during a week of August 2012 [38]. Particularly, the data have a granularity of 1 second, and represent the loading of secondary transformers. Line impedances, shunt admittances, as well as active and reactive loads are adopted from the respective data set. It is assumed that 18 PV systems are located at nodes $4$, $7$, $10$, $13$, $17$, $20$, $22$, $23$, $26$, $28$, $29$, $30$, $31$, $32$, $33$, $34$, $35$, and $36$, and their generation profiles are simulated based on the real solar irradiance data available in [38]. The ratings of these inverters are $300$ kVA for $i=3$, $350$ kVA for $i=15,16$, and $200$ kVA for the remaining inverters. Loads and the power available from a PV system with capacity of $50$ kW are reported in Fig. 2 for illustrative purposes.

The voltage limits $\overline{v}_{i}$ and $\underline{v}_{i}$ are set to $1.05$ p.u. and $0.95$ p.u. respectively, for $\forall i\in\mathcal{N}$. Various step sizes $\varepsilon_{1}$ and $\varepsilon_{2}$ are tested to provide examples of cases where the algorithm converges as well as cases where it is not convergent. The customers’ objective functions are set uniformly to $C_{i}^{t}(p_{i}^{t},q_{i}^{t})=c_{p}(p^{t}_{i,\textrm{av}}-p^{t}_{i})^{2}+c_{q}q_{i}^{t2}$, in an effort to minimize the amount of real power curtailed from the available power $p^{t}_{i,\textrm{pv}}$ based on irradiance conditions at time $t$, and the amount of reactive power injected or absorbed. The coefficients are set to $c_{p}=3$ and $c_{q}=1$. The network-oriented objective is set to $D(v^{t})=\frac{1}{2}\|v^{t}-v^{\text{nom}}\|_{2}^{2}$ to penalize voltage deviation from the nominal value $v^{\text{nom}}=1$ p.u. Without loss of generality, we demonstrate our results with the trade-off parameter $\gamma$ set to either [math] or $1$. For $\gamma=1$, it is possible to trade off the customer-oriented objectives for flatness of the voltage profile. The regularization parameter $\phi$ is set to $10^{-4}$.

V-B Iterative Algorithm

We first test Algorithm 1 and show how the algorithm can address overvoltages in distribution systems [4]. To this end, we focus on a single timeslot at 12 pm.

V-B1 Convergence

Let $\gamma=0$ for simplicity. Recall from Theorem 4 that step sizes $\varepsilon_{1}$ and $\varepsilon_{2}$ both affect the convergence properties. For simplicity, set $\varepsilon_{2}=0.01$, and consider tuning $\varepsilon_{1}$ to achieve convergence. Similar results can be observed by fixing $\varepsilon_{1}$ and tuning $\varepsilon_{2}$, or tuning both $\varepsilon_{1}$ and $\varepsilon_{2}$. As shown in Figure 3, when $\varepsilon_{1}$ is increased from $0.01$ to $0.3$, we observe faster convergence. However, when we further increase $\varepsilon_{1}$ beyond $0.4$, an oscillatory behavior is observed.

V-B2 Voltage regulation

The results are plotted in Figure 4 corresponding to the case where $\varepsilon_{1}=\varepsilon_{2}=0.01$. We show voltage profiles in three scenarios: (i) uncontrolled setting, where the PV systems operate at unity power factor and inject the maximum available power without any curtailment (blue dots), (ii) controlled voltages with $\gamma=0$ (red dots), and (iii) controlled voltages with $\gamma=1$ (yellow dots). It is clear that in the uncontrolled case (i) the voltage values exceed the limit of 1.05 p.u. (black dashed line) due to large reverse power flows, while the controlled scenarios (ii) and (iii) show voltage within limits. Furthermore, voltage values achieved by (iii) are closer to the nominal value than those by (ii), because (iii) also penalizes voltage deviation from 1 p.u.

V-C Online Algorithm

Next, Algorithm 2 is tested based on the irradiance and load profiles shown in Figure 2. One iteration (i.e., $K=1$) is performed every second (i.e., $h=1$ second). In the following, the performance of the proposed online algorithm are illustrated for both cases of $\gamma=0$ and $\gamma=1$. We will provide the voltage profiles as well as the profiles of the incentive signals. In addition, incentive signal profiles under $K=5$ are provided for comparative purpose.

V-C1 $\gamma=0$

In this case, the function $D(v)$ is disregarded. The voltage profiles obtained when the PV inverters operate according to business-as-usual practices and when they implement the proposed Algorithm 2 are provided for nodes 2, 28, and 35 in Fig. 7. In the uncontrolled case, voltage values exceed the upper limit during the mid-day hours because of the reverse power flows; in contrast, the proposed algorithm enforces voltage regulation, even when only one iteration is performed every second. Fig. 7 illustrates that the incentive signals become nonzero when voltages would violate the limits. The negative signals incentivize the customers to curtail active power and produce negative reactive power.

V-C2 $\gamma=1$

The voltage profiles obtained when $\gamma=1$ are plotted in Fig. 10. In this case, the voltage magnitudes are driven closer to the nominal value, at the cost of curtailing more real power and absorbing more reactive power. Voltages are clearly within limits.

V-C3 $K=5$

We repeat the simulations with five iterations per second (i.e., $K=5$), and plot the signal profiles based on the last iteration of each second. The results are presented in Fig. 7 and Fig. 10. As expected, the incentive signals generated with more iterations provide more accurate (see Fig. 7 vs Fig. 7) and more steady (see Fig. 10 vs Fig. 10) tracking of voltage changes. The resultant controlled voltage profiles with $K=5$ are omitted because they are not largely different from Fig. 7 and Fig. 10 with $K=1$; nevertheless, by examining the results in details we have found that voltage profiles with $K=5$ commit less voltage violation when $\gamma=0$, and enjoy smaller (time) variance in general.

VI Conclusion

This paper considers a time-varying social welfare maximization problem modeling network operator and DER-owners operational objectives as well as voltage constraints. The formulated problem is non-convex; however, we propose a convex relaxation and we provide conditions under which the optimal solutions of the relaxed problem coincide with the optimal points of the non-convex social-welfare problem. We then design distributed algorithms to identify the solutions of the time-varying social welfare maximization problem. An online algorithm is proposed to enable tracking of the solutions in the presence of fast time-varying operational conditions and changing optimization objectives. Stability of the proposed schemes is analytically established and numerically corroborated. Future research directions include the extension of the proposed framework to control DERs with discrete power levels and devices involving discrete decision variables.

Complementing the convergence results of Theorem 4, in the following we provide sufficient conditions on the stepsizes $\varepsilon_{1}$ and $\varepsilon_{2}$ that guarantee the operator $\hat{T}$ in (27) to be a contraction.

Theorem 6**.**

If the stepsizes $\varepsilon_{1}$ and $\varepsilon_{2}$ satisfy the following conditions for any $i\in\mathcal{N}$:

[TABLE]

*then $\hat{T}$ is a contraction. $\Box$ *

Proof.

Let $\nabla\hat{T}\in\mathbb{R}^{4N\times 4N}$ denote Jacobian matrix of $\hat{T}$, and let $\nabla\hat{T}_{ij}$ denote the element on row $i$ and column $j$ of matrix $\nabla\hat{T}$. To prove that $\hat{T}$ is a contraction, it is sufficient to have the following condition:

[TABLE]

which is satisfied if the following three inequalities hold:

[TABLE]

Conditions (42) and (43) are in fact equivalent. Therefore, (42) are sufficient for $\hat{T}$ to be a contraction.

Remark 6**.**

*Conditions (42) together with assumptions in this paper guarantee the existence of small enough step sizes $\varepsilon_{1}$ and $\varepsilon_{2}$ to achieve convergence. This result is consistent with Theorem 4. *