Distributed Dual Subgradient Methods with Averaging and Applications to Grid Optimization

Abstract

We study finite-time performance of a recently proposed distributed dual subgradient (DDSG) method for convex-constrained multi-agent optimization problems. The algorithm enjoys performance guarantees on the last primal iterate, as opposed to those derived for ergodic means for standard DDSG algorithms. Our work improves the recently published convergence rate of ${{\mathcal {O}}}(\log T/\sqrt{T})$ with decaying step-sizes to ${{\mathcal {O}}}(1/\sqrt{T})$ with constant step-size on a metric that combines sub-optimality and constraint violation. We then numerically evaluate the algorithm on three grid optimization problems. Namely, these are tie-line scheduling in multi-area power systems, coordination of distributed energy resources in radial distribution networks, and joint dispatch of transmission and distribution assets. The DDSG algorithm applies to each problem with various relaxations and linearizations of the power flow equations. The numerical experiments illustrate various properties of the DDSG algorithm–comparison with standard DDSG, impact of the number of agents, and why Nesterov-style acceleration can fail in DDSG settings.

Appendices

Appendix A: Dual Subgradient Methods Cannot Generally be Accelerated

Nesterov-type acceleration relies on smoothness of the objective function. Here, we illustrate why such acceleration is generally untenable in dual subgradient settings. We first present examples where the dual function for constrained optimization problems are nonsmooth. Then, we illustrate through examples how nonsmoothness of the objective function impairs acceleration. Consider an optimization problem of the form

(51)

with ${\varvec{\Xi }}$ being positive semidefinite (but not positive definite) and $\mathbb {X}$ being a convex polyhedral set. This is a quadratic program (QP) that simplifies to a linear program (LP) when ${\varvec{\Xi }}=0$. Associate multipliers ${\varvec{z}} \ge 0$ with ${\varvec{A}}\varvec{x} \le \varvec{b}$. Then, the dual function for the problem in (51) is given by

(52)

With parameter choices

(53)

we plot ${{\mathcal {D}}}({\varvec{z}})$ with ${\varvec{\Xi }} = 0$ and ${\varvec{\Xi }} = \mathop {\text {diag}}\left( 24, 26, 0 \right)$, respectively, in the left and right of Fig. 10. The primal and dual optimum is ${{\mathcal {P}}}^\star = {{\mathcal {P}}}^\star _{D} = 2.43$ for the QP and 2.30 for the LP.

Fig. 10

Illustrations of nonsmooth ${{\mathcal {D}}}({\varvec{z}})$ on top and a slice ${{\mathcal {D}}}(\cdot , z_2^\star )$ on the bottom for (51) with parameters in (53). The plots on the left are obtained with ${\varvec{\Xi }} = 0$ and on the right with ${\varvec{\Xi }} = \mathop {\text {diag}}\left( 24, 26, 0 \right)$

${{\mathcal {D}}}$ in Fig. 10 is nonsmooth. One might expect that for QP with ${\varvec{\Xi }} \ne 0$, ${{\mathcal {D}}}$ might be smooth. Indeed with positive definite ${\varvec{\Xi }}$, ${{\mathcal {D}}}$ becomes smooth, per [53, Lemma 2.2] or [6, Theorem 6.3.3]. Smoothness can no longer be guaranteed, when ${\varvec{\Xi }}$ is only positive semidefinite, as evidenced by Fig. 10.

To demonstrate the impact of nonsmoothness on acceleration, consider the problem

$$\begin{aligned} \underset{\varvec{x} \in \mathbb {R}^{2}}{\text {minimize}}&\quad f(\varvec{x}) := \frac{1}{2} \Vert \varvec{x} - {\varvec{x}_{\text {C}}}\Vert _2^2 + \lambda \Vert \varvec{x} \Vert _1. \end{aligned}$$

(54)

With ${\varvec{x}_{\text {C}}}= 0$, f is minimized at the origin. If $\lambda =0$, then f is smooth. For $\lambda >0$, f is nonsmooth at the origin. With ${\varvec{x}_{\text {C}}}\ne 0$, f is again smooth with $\lambda = 0$. However, with $\lambda > 0$, f is nonsmooth at the origin, but not at the optimum of f.

We compare (sub)gradient descent, i.e., $\varvec{x}(t+1) = \varvec{x}(t) - \eta \nabla f(\varvec{x}(t))$, with an accelerated variant described in Algorithm 4 on (54).

Algorithm 4

Accelerated Gradient Descent adopted from [71], Sect. 3.7.2

We start both algorithms at (0.05, 0.05) with $\eta = 0.003$. For accelerated (sub)gradient descent (AGD), we also set $\alpha (1) = 0$. The plot on the left side of Fig. 11 shows the progress of subgradient descent (SGD) and AGD on (54) with ${\varvec{x}_{\text {C}}}= 0$. Note that when $\lambda = 0$, i.e., when f is smooth, AGD outperforms SGD. However, when $\lambda = 0.01$ and f is nonsmooth at the optimum, AGD performs better initially, but both algorithms oscillate around the optimum with similar errors, eventually.

In the right subgraph of Fig. 11, we compare SGD and AGD under the same settings, but with ${\varvec{x}_{\text {C}}}= (0.1,0.1)$. Note that AGD now performs better than SGD with zero and nonzero $\lambda$. That is, when the nonsmoothness is away from the optimum, AGD can accelerate convergence to the optimum locally, as long as the iterates remain within a region around the optimum where the function is locally smooth.

Fig. 11

Comparison of SGD and AGD on (54) with ${\varvec{x}_{\text {C}}}=0$ on the left and with ${\varvec{x}_{\text {C}}}=(0.1, 0.1)$ on the right

Appendix B: Simulation Details for Sections 4, 5 and 6

Network data were obtained from MATPOWER 7.1.

1.1 Appendix B1: Data for Solving ${{\mathcal {P}}}_1$

The multi-area power system considered in Section 4 is illustrated in Figure 1. The 118-bus networks were modified as follows. Tie-line capacities were set to 100MW and their reactances were set to 0.25p.u. Capacities of transmission lines internal to each area were set to 100MW. All loads and generators at boundary buses were removed. Quadratic cost coefficients were neglected and the linear cost coefficients ${\varvec{c}}_j$ of the generators were perturbed to $\widetilde{{\varvec{c}}}_j := {\varvec{c}}_j \circ \left( 0.99 + 0.02 {\varvec{\xi }}_j \right)$, for $j = 1,\ldots ,N$, where entries of ${\varvec{\xi }}_j$ are independent ${{\mathcal {N}}}(0,1)$ (standard normal) variables. All phase angles were restricted to $[-\frac{\pi }{6}, \frac{\pi }{6}]$.

1.2 Appendix B2: Data for Solving ${{\mathcal {P}}}_2$

The 4-bus network considered in Section 5, shown in Figure 6, is modified from the IEEE 4-bus network as follows. The branch joining buses 1 and 4 was altered to connect buses 3 and 4. We enforced squared current flows as $\ell _{j,k} \in [0, 200]$ Amp$^2$, real and reactive branch power flows as $P_{j,k} \in [-1, 1]$ MW and $Q_{j,k} \in [-1, 1]$ MVAR, respectively. DER generators were added at buses 2, 3 and 4. Bus 1 defined the T &D interface. Generation capacities were fixed to [0, 1] MW and $[-1, 1]$ MVAR. Generation costs were $\alpha _{pj} (p^G_j)^2 + \beta _{pj} p^G_j + \alpha _{qj} (q^G_j)^2$ with coefficients in Table 2.

Table 2 Cost coefficients of the 4-bus network for ${{\mathcal {P}}}_2$.

For the IEEE 15-bus system shown in Figure 5, we modified the branch flow limits to mirror those for the 4-bus system. We added 7 distributed generators at buses 5, 7, 8, 10, 13, 14, 15, where bus 1 is the T &D interface, all with capacities [0, 0.2] MW and $[-0.2, 0.2]$ MVAR. Generation costs were similar to the 4-bus network with coefficients in Table 3.

Table 3 Cost coefficients of the 15-bus network for ${{\mathcal {P}}}_2$.

We randomized the real and reactive power demands at each change point by scaling each (real/reactive) load by $\left[ \omega + (\omega ' - \omega ) \xi \right]$, where $\xi \sim {{\mathcal {N}}}(0,1)$. Parameters $(\omega , \omega ')$ were varied at the change points in the sequence (0.70, 1.30), (0.80, 1.20), (0.85, 1.15), (0.75, 1.20), (0.95, 1.05). The experiment was initialized with default loads from MATPOWER.

1.3 Appendix B3: Data for Solving ${{\mathcal {P}}}_3$

In Section 6, for the 204-bus system in Figure 7, the 6-bus transmission network was modified as follows. All branch capacities are set to 5MW. All real and reactive generation capacities were set to [0, 5]MW and $[-5, 5]$ MVAR, respectively. We considered ${P}_\ell ^D + j Q_\ell ^D = (4 + j 4) \text {[MVA]}$ at each bus $\ell = 1,\ldots ,6$. Generation costs were similar to the 4-bus network with coefficients in Table 4.

Table 4 Cost coefficients of the 6-bus network for ${{\mathcal {P}}}_3$.

For all 33-bus distribution networks, all branch capacities were set to 4 MW. Four DER generators were added at buses 18, 22, 25 and 33. Bus 1 is the T &D interface. Again, we considered generation costs as for ${{\mathcal {P}}}_2$ but with coefficients ${\varvec{\alpha }}_{p\ell } = 5 \circ \left( 0.9 + 0.1 {\varvec{\xi }}_\ell \right)$, ${\varvec{\beta }}_{p\ell } = 20 \circ \left( 0.9 + 0.1 {\varvec{\xi }}'_\ell \right)$ and ${\varvec{\alpha }}_{q\ell } = 3 \circ \left( 0.9 + 0.1 {\varvec{\xi }}''_\ell \right)$ for $\ell = 1,\ldots ,n^{\text {tran}}$, where all entries of ${\varvec{\xi }}_\ell$, ${\varvec{\xi }}'_\ell$, ${\varvec{\xi }}''_\ell$ are drawn from ${{\mathcal {N}}}(0,1)$. Real and reactive power demands in the distribution networks were randomized similarly to that for ${{\mathcal {P}}}_2$ with $\omega =0.9$ and $\omega '=1.1$.