Online ADMM for Distributed Optimal Power Flow via Lagrangian Duality

Abstract

At present, the power system has the characteristics of mutual independence but interconnection, and the interconnection between the various subsystems brings certain challenges to the distributed computing of the power grid. In addition, a substantial amount of naturally uncertain renewable resources are incorporated into the power system, which will impose volatile dynamics on the grid. In this paper, an alternating direction multiplier method (ADMM) is proposed for the power system with real-time renewables to tackle the online optimal power flow (OPF) problem. Due to the adoption of the Lagrangian duality, the proposed distributed ADMM scheme utilizes consensus ADMM to solve the dual OPF problem, which only discloses boundary coupling via the Lagrangian multiplier and further reduces the amount of information communication. Given the natural uncertainty of distributed energy resources (DER), the algorithm avoids the double-loop implementation or the uncertainty of traditional distributed methods of using the boundary information as equality constraints caused by dynamic DER. It is thus capable of providing a provable performance guarantee and is inherently developed to cope with the dynamic OPF problem with renewables in an online fashion. Taking the IEEE 30-bus system as a test feeder, the simulation results verify the efficiency and robustness of the proposed algorithms in solving both the static and dynamic OPF problems; in addition, the online method can effectively avoid the violent fluctuations of the conventional generator output copying with renewables rapid variation in comparison with the offline algorithms.

1. Introduction

Since the pioneering work on optimal power flow [1] was presented by Carpentier in 1962, his calculation model containing power flow equality and inequality constraints has been an important reference paradigm for realizing the reliable and efficient operation of power systems. This work has attracted attention and obtained a vast amount of research results; please see [2] and the references therein.

OPF aims to determine the economical operating parameters for power systems, minimizing various costs, such as generation cost, line losses, and even consumer disutility. Due to the quadratic relation between the power-related quantities and the voltage phasors [3], the OPF constitutes a nonconvex quadratically constrained quadratic program (QCQP) problem. To deal with this challenge, some (non)linear programming approaches based on gradient descent, Newton-Raphson, and interior-point methods [4] have been put forward to seek solutions satisfying the Karush–Kuhn–Tucker (KKT) condition for local optimality. However, these solutions are often sensitive to initialization. By introducing semidefinite programming (SDP) [2,3,5,6], second-order cone programming (SOCP), and other convex relaxation methods [7] to perform relaxation transformation on the power flow constraints, the nonconvex OPF can be easily resolved. Specifically, the optimal solution of SDP relaxation is a lower bound on the optimal value of the original nonconvex QCQP problem [8].

1.1. Primary Motivations

In the past, solving OPF was usually based on centralized methods [4,6], which assumes that all nodes in the power systems transmit their information to a central operating unit. Actually, the distribution region in a power grid may belong to different owners, and some sensitive information intercepted by malicious attackers might cause financial losses or even security risks. With the wide application of distributed optimization, this situation will be thoroughly changed. As a distributed algorithm, ADMM is naturally suitable for distributed convex optimization, especially for area-based large-scale problems, such as OPF [9,10,11,12,13] and control of reactive power [7,14,15], as well as minimum voltage deviation. Furthermore, ADMM can effectively reduce information disclosure and further preserve local privacy [16].

Currently, the share of variable renewable energy resources, such as wind and solar, is increasing in power systems [17]. Because these DERs are dependent on weather conditions, a variety of technical challenges, such as unpredictability, uncontrollability, and fluctuations in voltage and reactive power [18,19], must be overcome to safely connect them to the grid while maintaining stability and reliability. When integrating DERs into the grid, the power flow distribution will be changed, and the boundary information of each area will also be changed along with it [20]; the original static OPF thus evolves into a dynamic OPF covering multiple continuous time sections [21]. The static OPF finds the optimal solution at a certain time section, regardless of the interrelation among time sections. In the power system with renewables, when the output power of DERs changes rapidly, the generator output in contiguous time sections changes greatly. If continuously adopting the static OPF, the result may cause the generator output to exceed its regulation rate limit between two contiguous time sections, which makes the result infeasible in practice. Therefore, it is imperative to build robustness into power system operation for the dynamic OPF problem against the uncertainties of DERs. Stochastic and robust optimization techniques have been employed to cope with this uncertainty and minimize operational costs [22]. However, these approaches typically require characterizing the uncertainty bounds or the distribution of the DER’s generation output.

1.2. Related Work and Research Gap

Following the above motivation, we survey some related work with the aim of giving a clear understanding of how our algorithm and the results relate to and, in many cases, improve upon it. Our paper is closest to the real-time OPF [23,24,25,26] or online OPF [3,27,28]. For the real-time OPF problem, ref. [23] presented an online gradient projection algorithm, in which the proposed algorithm converges to the set of local optima and an upper bound on the suboptimality gap of any local optimum is also provided as a sufficient condition to guarantee the algorithm converges to a global optimum. Subsequently, ref. [24] developed a real-time algorithm based on quasi-Newton methods, which utilizes second-order information to provide suboptimal power flow solutions. Further extensions of this topic were presented in [25,26] with a look toward the treatment of rapid and economical responses to the changes in the power system operating state by using a deep reinforcement learning method. For the online OPF problem, ref. [3] first sought to study this topic by using an online mirror descent-based algorithm. Besides this, ref. [27] also explored four strategies, including pure numerical, closed form, and numerical hybrids, as well as a Lagrange multiplier method, in which the authors further considered the tradeoff between the time-cost of the solution and guaranteeing a feasible solution. Recently, ref. [28] developed an online algorithm based on the second-order Taylor information in a distributed way by decomposing the sensitivities of the OPF targets among different areas. For more details, please see Table 1.

Despite these efforts, however, there still exist two significant limitations. First, since the OPF problem is nonconvex and may be NP-hard in the worst case, some gradient-based methods, such as [23,24], only guarantee the provision of suboptimal power flow solutions or a local optimum, and guaranteeing the global optimum requires an extra sufficient condition. Although Lavaei and Low [6] have invented a zero duality gap approach used in an SDP optimization for the static OPF problem, it still needs an extra sufficient condition. Second, although the Lagrange multiplier method has been proven in [27] to achieve the fastest solution time for a large-size problem, there are still some underlying problems when applying this method to the real-time or online OPF, such as the decomposition of the large-size OPF into several subproblems, and providing the global optimum for the nonconvex OPF model, as well as guaranteeing zero duality gap without any extra condition, and so on.

1.3. Statement of Contributions

In this paper, inspired by consensus optimization [12,29], we propose an online convex optimization scheme integrating ADMM and Lagrangian duality for resolving the online OPF problem. First, for the static OPF problem, the cost-minimization model referring to [2,11] is built, in which, power flow constraints are constructed via boundary equations comprising the voltage, the phase angle, and the coupling line parameters (admittance, etc.) on the boundary nodes among adjacent areas. To implement the distributed ADMM, we propose a simple area partitioning method based on spectral clustering to provide an optimal partition of the power system. By introducing the Lagrangian duality on the original OPF [30], a global consensus variable is employed as the Lagrangian multiplier for each sub-problem, which is simultaneously copied as a local variable for each sub-area solution [31]. After doing that, the distributed ADMM applied only discloses the boundary coupling information rather than the internal node information so as to avoid determining the coupling equation separately and including the boundary information in the objective function. For the online OPF with renewables, we propose an online convex optimization scheme to tackle the uncertainties of DERs. Specifically, at the beginning of each time slot, the power output of renewable generators must be determined (and fixed) without knowledge of loads; and at the end of the time slot, the loads are revealed so that the power output of conventional generators can be determined under the ramp constraints. The main contributions of the paper are as follows:

• The Lagrangian duality is utilized to couple the power flow equation and boundary information instead of directly establishing boundary coupling as the ADMM constraint, it can be well adapted to the real-time nature of DESs, such as solar energy, and it can also disclose less boundary information.
• For the SDP relaxation model of the static OPF, the distributed ADMM is compatibly employed to solve its convex Lagrangian duality problem instead of solving the iterative sub-problem of ADMM under the original nonconvex OPF model.
• In order to adapt the distributed ADMM to the dynamic OPF with renewables, an online scheme is proposed to cope with the uncertainties of DERs, in which the double-loop implementation is avoided, thus providing a provable performance guarantee.

Notations: The notation “$\mathbf{i}$” is reversed for the imaginary unit. $\mathrm{Re}\left\{A\right\}$, $\mathrm{Im}\left\{A\right\}$, and $A^{*}$ represent the real part, imaginary part, and complex conjugation of a given complex scalar A, respectively. The notations ${\mathbf{A}}^{\top }$, ${\mathbf{A}}^{\mathrm{H}}$, and $\mathrm{rank}\left\{\mathbf{A}\right\}$ represent the transpose, Hermitian transpose, and rank of a given matrix $\mathbf{A}$, respectively. For given the notation $\left|x\right|$ of absolute-value for a scalar $x\in \mathbb{R}$, the ${\ell }_p$-norm ($p\ge 1$) is defined by ${\parallel \mathbf{x}\parallel }_p:=\left(\right|x_1{|}^p+\cdots +|x_N{{|}^p)}^{1/p}$ for the vector $\mathbf{x}={[x_1,\cdots ,x_N]}^{\top }\in {\mathbb{R}}^N$, which yields the ${\ell }_1$-norm when $p=1$ and the Euclidean norm when $p=2$, as well as the ${\ell }_{\infty }$-norm when $p\to \infty$. Without special statement, $\parallel ·\parallel$ commonly denotes the Euclidean norm. The notation $\langle \mathbf{x},\mathbf{y}\rangle$ denotes the inner product given by $\langle \mathbf{x},\mathbf{y}\rangle :={\mathbf{x}}^{\top }\mathbf{y}={\sum }_{i=1}^N{x}_i{y}_i$, for $\mathbf{x},\mathbf{y}\in {\mathbb{R}}^N$. Finally, we define the infimum of a set $\mathcal{S}$ as the greatest scalar that is less than or equal to all elements of $\mathcal{S}$, abbreviated as $inf\mathcal{S}$; and denote the argument of the minimum and argument of the maximum for an objective function f by $argminf$ and $argmaxf$, respectively.

Table 1. Comparisons of the proposed algorithm with the existing distributed online algorithms.

Literature	Topic	AC/DC	Renewables	Approach
[30]	Dynamic OPF	DC	N/A	Distributed dual consensus ADMM
[23]	Real-time OPF	AC	N/A	Online gradient projection
[24]	Real-time OPF	AC	N/A	Quasi-Newton methods
[25,26]	Real-time OPF	AC	N/A	Deep reinforcement learning
[27]	Online OPF	DC	N/A	Lagrange multiplier method
[3]	Online OPF	AC	Wind turbine	Online mirror descent
[28]	Online OPF	AC	Photovoltaic	Second-order Taylor-based gradient
ours.	Online OPF	AC	Photovoltaic	Online ADMM Lagrangian duality

Table 1. Comparisons of the proposed algorithm with the existing distributed online algorithms.

Literature	Topic	AC/DC	Renewables	Approach
[30]	Dynamic OPF	DC	N/A	Distributed dual consensus ADMM
[23]	Real-time OPF	AC	N/A	Online gradient projection
[24]	Real-time OPF	AC	N/A	Quasi-Newton methods
[25,26]	Real-time OPF	AC	N/A	Deep reinforcement learning
[27]	Online OPF	DC	N/A	Lagrange multiplier method
[3]	Online OPF	AC	Wind turbine	Online mirror descent
[28]	Online OPF	AC	Photovoltaic	Second-order Taylor-based gradient
ours.	Online OPF	AC	Photovoltaic	Online ADMM Lagrangian duality

2. Model and Formulation

2.1. Optimal Power Flow Formulation

Considering a power network with the set of buses $\mathcal{N}:=\{1,2,$$\cdots ,N\}$, the set of flow lines $\mathcal{E}\subseteq \mathcal{N}\times \mathcal{N}$ between the buses. Let ${\Omega }_n:=\left\{n^{\prime }\right|(n,n^{\prime })\in \mathcal{E}\}$ be the set of buses connected to bus n, where $n\in \mathcal{N}$ is the index. Denote the complex voltage amplitude and the injected current at bus n by $V_n$ and $I_n$, respectively. Each line $(n,n^{\prime })\in \mathcal{E}$ of the network is modeled as a passive device with an admittance $y_{n{n}^{\prime }}$, thus the network can be modeled as a general admittance matrix $\mathbf{Y}\in {\mathbb{C}}^{N\times N}$ whose $(n,n^{\prime })$-entry is formed as

$$Y_{n{n}^{\prime }}=\left\{\begin{array}{ccc}{y}_{n0}+\sum _{n^{\prime }\in {\Omega }_n}{y}_{n{n}^{\prime }},\hfill & & \mathrm{if}n=n^{\prime }\hfill \\ -y_{n{n}^{\prime }},\hfill & & \mathrm{if}n\ne n^{\prime },(n,n^{\prime })\in \mathcal{E}\hfill \\ 0,\hfill & & \mathrm{otherwise}\hfill \end{array}\right.$$

(1)

where $y_{n0}$ is the admittance-to-ground of bus n, which is normally not considered for the OPF modeling. Then, Kirchhoff’s and Ohm’s laws yield

$$I_n=\sum _{n^{\prime }\in {\Omega }_n}(V_n-V_{n^{\prime }})y_{n{n}^{\prime }},$$

(2)

which means the complex power injected to bus n can be given by

$$S_n:=P_n+\mathbf{i}{Q}_n=V_n{I}_n^{*}=V_n\sum _{n^{\prime }\in {\Omega }_n}{(V_n-V_{n^{\prime }})}^{*}{y}_{n{n}^{\prime }}^{*},$$

(3)

where $P_n$ and $Q_n$ represent the active power and reactive power injected to bus n, respectively.

We denote the set of buses with conventional generators and renewable generators by ${\mathcal{N}}_G\subseteq \mathcal{N}$ and ${\mathcal{N}}_R\subseteq \mathcal{N}$, respectively. Let $P_n^G$ and $Q_n^G$ be the active power output and reactive power output of the conventional generator at bus $n\in {\mathcal{N}}_G$, respectively. Similarly, let $P_n^R$ and $Q_n^R$ be the active power output and reactive power output of renewable power, such as photovoltaic energy, wind power, and geothermal energy, generated by bus $n\in {\mathcal{N}}_R$, respectively. Finally, let $P_n^D$ and $Q_n^D$ be the active load and reactive load at bus $n\in \mathcal{N}$, respectively. Then, the power balance condition of each bus $n\in \mathcal{N}$ can be formed as

$$\begin{array}{cc}\hfill P_n& =P_n^G+P_n^R-P_n^D,\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill Q_n& =Q_n^G+Q_n^R-Q_n^D.\hfill \end{array}$$

(4b)

It should be noted here that $P_n^G=Q_n^G=0$, if $n\in \mathcal{N}$ while $n\notin {\mathcal{N}}_G$, and all else follows.

Ref. [2] has formed four types of line capacity constraints, which are can achieve various goals in practice, such as avoiding line overheating and guaranteeing the stability of the network. To simplify the problem, researchers generally select one of them as the line capacity constraint [3]:

$$\left|V_n-V_{n^{\prime }}\right|\le \Delta V_{n{n}^{\prime }}^{max}$$

(5)

where $\Delta V_{n{n}^{\prime }}^{max}$ represents the maximum voltage difference over line $(n,n^{\prime })$, this limitation is equivalent to the line current magnitude constraint, since each line has been modeled as a simple admittance and it follows from (2) that $V_n-V_{n^{\prime }}$ is proportional to the line current.

Now, we turn to formulating the OPF problem. We consider the minimization of system operating cost as the objective function, and the power flow balance, line capacity limitation, generator capacity limitation, and voltage amplitude limitation as constraints. Specifically, the mathematical formulation of OPF is given as

$$\begin{array}{cc}\hfill \underset{\mathbf{V},{\mathbf{P}}^G,{\mathbf{Q}}^G}{minimize}& \sum _{n\in {\mathcal{N}}_{\mathcal{G}}}{f}_n\left(P_n^G\right)\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& P_n^G+P_n^R-P_n^D=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Re}\left\{V_n(V_n^{*}-V_{n^{\prime }}^{*})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(6b)

$$\begin{array}{cc}\hfill & Q_n^G+Q_n^R-Q_n^D=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Im}\left\{V_n(V_n^{*}-V_{n^{\prime }}^{*})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(6c)

$$\begin{array}{cc}\hfill & P_n^{min}\le P_n^G\le P_n^{max}\hfill \end{array}$$

(6d)

$$\begin{array}{cc}\hfill & Q_n^{min}\le Q_n^G\le Q_n^{max}\hfill \end{array}$$

(6e)

$$\begin{array}{cc}\hfill & V_n^{min}\le \left|V_n\right|\le V_n^{max}\hfill \end{array}$$

(6f)

$$\begin{array}{cc}\hfill & \left|V_n-V_{n^{\prime }}\right|\le \Delta V_{n{n}^{\prime }}^{max}\hfill \end{array}$$

(6g)

where $\mathbf{V}$, ${\mathbf{P}}^G$, ${\mathbf{P}}^R$, ${\mathbf{P}}^D$, ${\mathbf{Q}}^G$, ${\mathbf{Q}}^R$, and ${\mathbf{Q}}^D$ are the vector forms of $V_n$, $P_n^G$, $P_n^R$$P_n^D$, $Q_n^G$, $Q_n^R$, and $Q_n^D$ for $n\in \mathcal{N}$, respectively. Equations (6b) and (6c) are the power balance conditions accounting for the conservation of power at bus n, yielded by (3), (4a), and (4b). Equations (6d), (6e), and (6f) restrict the active power, reactive power, and voltage magnitude at bus n, respectively, for given limits $P_n^{min}$, $Q_n^{min}$, $V_n^{min}$, $P_n^{max}$, $Q_n^{max}$, and $V_n^{max}$. Objective (6a) can capture various costs, such as generation cost, line losses, conservation voltage reduction, etc. Generally, the generation cost is modeled by a quadratic function form as:

$$f_n\left(P_n^G\right):={\alpha }_n{\left(P_n^G\right)}^2+{\beta }_n\left(P_n^G\right)+{\gamma }_n$$

(7)

where ${\alpha }_n\ge 0$, ${\beta }_n$, and ${\gamma }_n$ represent the cost parameters of the generator at bus $n\in {\mathcal{N}}_{\mathcal{G}}$. In a nutshell, OPF minimizes the total cost ${\sum }_{n\in \mathcal{N}}{f}_n$ over the unknown parameters $\mathbf{V}$, ${\mathbf{P}}^G$, and ${\mathbf{Q}}^G$, for given the known vector ${\mathbf{P}}^R$, ${\mathbf{P}}^D$, ${\mathbf{Q}}^R$, and ${\mathbf{Q}}^D$, subject to some constraints.

2.2. Convex Relaxation Approach to OPF

It should be noted here that objective (6a) of OPF is convex while the constraints are nonconvex due to the nonlinear terms $\left|V_n\right|$ and $V_n{V}_{n^{\prime }}^{*}$. Since solving nonconvex optimization is pretty hard, ref. [6] suggested a convex relaxation of the OPF problem by defining $\mathbf{W}:=\mathbf{V}{\mathbf{V}}^{\mathrm{H}}$, then OPF constraints containing the quadratic item $V_n{V}_{n^{\prime }}^{*}$ can be expressed as a linear function of $\mathbf{W}$. In order to ensure that the nonconvex OPF is transformed into an equivalent form, it requires that $\mathbf{W}$ is Hermitian and positive semidefinite with $\mathrm{rank}\left\{\mathbf{W}\right\}=1$. Hence, the only source of nonconvexity in the equivalent OPF problem comes from the ‘rank-1’ constraint. Motivated by the technique of convex relaxation, after dropping the rank constraint $\mathrm{rank}\left\{\mathbf{W}\right\}=1$, the nonconvex OPF problem has been formulated as a rank-relaxed semidefinite programming problem:

$$\begin{array}{cc}\hfill \underset{\mathbf{W}⪰0,{\mathbf{P}}^G,{\mathbf{Q}}^G}{minimize}& \sum _{n\in {\mathcal{N}}_{\mathcal{G}}}{\alpha }_n{\left(P_n^G\right)}^2+{\beta }_n\left(P_n^G\right)+{\gamma }_n\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& P_n^G+P_n^R-P_n^D=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Re}\left\{(W_{nn}-W_{n{n}^{\prime }})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(8b)

$$\begin{array}{cc}\hfill & Q_n^G+Q_n^R-Q_n^D=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Im}\left\{(W_{nn}-W_{n{n}^{\prime }})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(8c)

$$\begin{array}{cc}\hfill & P_n^{min}\le P_n^G\le P_n^{max}\hfill \end{array}$$

(8d)

$$\begin{array}{cc}\hfill & Q_n^{min}\le Q_n^G\le Q_n^{max}\hfill \end{array}$$

(8e)

$$\begin{array}{cc}\hfill & {\left(V_n^{min}\right)}^2\le W_{nn}\le {\left(V_n^{max}\right)}^2\hfill \end{array}$$

(8f)

$$\begin{array}{cc}\hfill & W_{nn}+W_{n^{\prime }{n}^{\prime }}-W_{n{n}^{\prime }}-W_{n^{\prime }n}\le {\left(\Delta V_{n{n}^{\prime }}^{max}\right)}^2\hfill \end{array}$$

(8g)

which can be solved using the interior-point method. Although the convex relaxation approach has nearly changed the original problem, if the original problem possesses a certain nice structure and the relaxation is also carefully performed, then the relaxation gaps are almost absent, i.e., solutions to the relax problem would be optimal for the original nonconvex problem as well.

2.3. Area Partitioning Based on Spectral Clustering

The scale of actual power grids is usually large, so the overall OPF problem of the whole power grid is very complex and is hard to solve directly. Fortunately, the voltage issues are usually regional [11]. Therefore, we can try to partition the whole power grid into several sub-regions such that the OPF problem can be decomposed into several sub-problems, thus decreasing the difficulty of solving the overall OPF problem.

In this subsection, we propose a simple partitioning method based on spectral clustering [32] for providing an optimal partition of the power system. Spectral clustering is a graph partitioning method that takes into account the affinity between any two nodes in the target graph, which in our considered problem, corresponds to the computational coupling between any two buses. Specifically, such coupling relation is reflected in the admittance matrix, but the admittance matrix cannot be regarded as the Laplacian matrix directly due to the complex admittance. Hence, modular operation is essential for constructing the Laplacian matrix. In the proposed method, the k-means algorithm [33] helps to implement clustering so that the inaccuracy in clustering nonconvex regions can be avoided. First, we fix the number of sub-regions by K, then find the optimal way to partition the system into this number of fixed sub-regions. The whole procedure is given by:

(1).

Construct the Laplacian matrix $\mathcal{L}$. Denote the nondiagonal entry of $\mathcal{L}$ by the negative modulus of its corresponding complex admittance and the diagonal entry of $\mathcal{L}$ by the sum of the complex admittance modulus, i.e.,

$${\mathcal{L}}_{n{n}^{\prime }}=\left\{\begin{array}{ccc}\sum _{n^{\prime }\in {\Omega }_n}\left|y_{n{n}^{\prime }}\right|,\hfill & & \mathrm{if}n=n^{\prime }\hfill \\ -|y_{n{n}^{\prime }}|,\hfill & & \mathrm{if}n\ne n^{\prime },(n,n^{\prime })\in \mathcal{E}\hfill \\ 0,\hfill & & \mathrm{otherwise}.\hfill \end{array}\right.$$

(2).

Find the K largest eigenvalues of $\mathcal{L}$ and form the matrix $\mathbf{U}\in {\mathbb{R}}^{N\times K}$ by stacking the corresponding eigenvectors $u_1,\cdots ,u_K$ as columns.

(3).

Let $v_n\in {\mathbb{R}}^K$ be the vector corresponding to the n-th row of $\mathbf{U}$, $n=1,\cdots ,N$. Cluster the points ${\left\{v_n\right\}}_{n=1}^N$ with the k-means algorithm into clusters $C_1,\cdots ,C_K$.

(4).

Assign bus n to cluster $C_k$ if the n-th row of $\mathbf{U}$ was assigned to cluster $C_k$.

3. Distributed ADMM for Static OPF

3.1. Lagrangian Duality Based on Partition

As described above, the power system can be divided into multiple sub-regions, and a local OPF problem can be formulated for each sub-region. Assuming that there are a total of K sub-regions. Denote the set of buses in the k-th sub-region by ${\mathcal{R}}_k$ satisfying ${\mathcal{R}}_k\cap {\mathcal{R}}_{\ell }=\varnothing$ if $k\ne \ell$, and the number of buses in the k-th sub-region by $|{\mathcal{R}}_k|$. For region k, define the variable ${\mathbf{x}}_k:=\{(W_n,P_n,Q_n)\mid n\in {\mathcal{R}}_k\}$, which covers the variable $x_n$ on each node n, and represents the original variable of each region, the OPF problem can be re-formulated by a local sum form based on variable $x_k$ as:

$$\begin{array}{cc}\hfill \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{minimize}& \sum _{k=1}^K{F}_k\left({\mathbf{x}}_k\right):=\sum _{k=1}^K\sum _{n\in {\mathcal{R}}_k}{f}_n\left(x_n\right)\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& g_k\left({\mathbf{x}}_k\right)=0\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill & h_k\left({\mathbf{x}}_k\right)\le 0\hfill \end{array}$$

(9c)

where ${\mathcal{X}}_k$ represents the local feasibility constraint set defined by ${\mathcal{X}}_k:=\left\{{\mathbf{x}}_k\mid \left(8b\right)-\left(8g\right)\right\}$, the power flow constraint $g_k\left({\mathbf{x}}_k\right)$ is expressed as equation constraints in the original OPF problem, i.e., (8b) and (8c), and the line transmission capacity limit $h_k\left({\mathbf{x}}_k\right)$ is expressed as inequality constraints in the original OPF problem, i.e., (8d)–(8g).

After the OPF model of each sub-region is given, a common practice is to duplicate the voltage phase angle or the transmitted power at boundary buses of each region [11]. This approach will remove the tie lines between connected regions, which makes the OPF problem more easily solved in a decentralized manner. However, it also has one major disadvantage; it needs to determine more coupling conditions. Especially if renewable energy integrates into the power grid, which will lead to power flow fluctuations, this approach may face some difficulties in solving the dynamical OPF problem under multiple time sections.

In order to avoid the determination of the coupling equations, we define the Lagrangian $D_k({\mathbf{x}}_k;\lambda ,\mu )$ associated with the original problem (9), which takes the constraints in (9) into account by augmenting the objective function with a sum of the constraint functions. Mathematically, it is given by:

$$D_k({\mathbf{x}}_k;\lambda ,\mu ):=\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+\lambda \sum _{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+\mu \sum _{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)$$

(10)

where the global variables $\lambda \ge 0$ and $\mu$ are the Lagrange multipliers that are introduced for the equality constraint and the inequality constraint, respectively. Now we define the Lagrangian duality Function (LDF) as the minimum value of the Lagrangian over ${\mathbf{x}}_k$:

$$\mathrm{LDF}(\lambda ,\mu )=\underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{inf}\sum _{k=1}^K{D}_k({\mathbf{x}}_k;\lambda ,\mu ).$$

(11)

which aims at solving the Lagrangian (10) with respect to ${\mathbf{x}}_k$ by maximizing $\mathrm{LDF}(\lambda ,\mu )$ with respect to $\lambda$ and $\mu$ rather than ${\mathbf{x}}_k$. Because the original problem (9) is convex, according to Slater constraint qualification [34], the duality gap between (9) and its dual problem is zero. Therefore, we can obtain the solutions to problem (9) by solving its dual problem $max\mathrm{LDF}(\lambda ,\mu )$. The advantage of this is that the constraints in (8b), (8c), and (8g) (or equivalently (9b) and (9c)) are coupled, and it is difficult to decompose these constraints into several independent subproblems, but the dual function can be easily decomposed into several independent maximizations based on region partitioning.

3.2. Distributed ADMM for Lagrangian Duality

In order to distributedly solve the Lagrangian duality, this paper adopts ADMM [16], which can divide a complicated problem into multiple simple sub-problems, and each sub-problem can be solved independently and locally by the data provider. Thus it is suitable for dealing with the Lagrangian duality based on area partitioning. However, the local Lagrangian (10) contains two global variables $\lambda$ and $\mu$ that are not desired. To this end, we introduce a local variable ${\mathbf{z}}_k$ for the k-th region to copy the global variable $\mathbf{y}:={[\lambda ;\mu ]}^{\top }$ as an equation constraint, then the Lagrangian duality problem can be transformed into an equivalent distributed minimization problem with local variable ${\mathbf{z}}_k:={[{\lambda }_{z,k};{\mu }_{z,k}]}^{\top }$ and consensus variable $\mathbf{y}$:

$$\begin{array}{cc}\hfill \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{minimize}& \sum _{k=1}^K-D_k({\mathbf{x}}_k;{\mathbf{z}}_k)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \mathbf{y}-{\mathbf{z}}_k=0\hfill \end{array}$$

(12)

which has been suitable for ADMM to solve distributedly. We denote the Lagrangian multiplier by $\sigma :={[{\sigma }_1,\dots ,{\sigma }_K]}^{\top }$ with ${\sigma }_k:={[{\lambda }_{\sigma ,k};{\mu }_{\sigma ,k}]}^{\top }$ and the penalty factor by $\rho >0$, one can get the augmented Lagrangian function of problem (12):

$$L(\mathbf{y},\mathbf{z};\sigma ):=\sum _{k=1}^K\left[-D_k({\mathbf{x}}_k;{\mathbf{z}}_k)+\langle {\sigma }_k,\left(\mathbf{y}-{\mathbf{z}}_k\right)\rangle +\frac{\rho }{2}{∥\mathbf{y}-{\mathbf{z}}_k∥}^2\right].$$

where $\mathbf{z}:={[{\mathbf{z}}_1,\dots ,{\mathbf{z}}_K]}^{\top }$. For a given interval $\rho \in [0,\overline{\rho }]$, one can evolve the penalty factor $\rho$ by a similar Logarithmic barrier indicator function, i.e., $\mathbb{I}\left(\rho \right)=-(1/c)log(\rho -\overline{\rho })$ with a constant $c=1/2$, 1 or 2. Note that in the method of multipliers, updating the dual variable uses a step-size equal to $\rho$. ADMM consists of the following three iterations with index $\tau$:

$$\begin{array}{cc}\hfill & {\mathbf{y}}^{\tau +1}:=\underset{\mathbf{y}}{argmin}L(\mathbf{y},{\mathbf{z}}^{\tau };{\sigma }^{\tau })\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill & {\mathbf{z}}_k^{\tau +1}:=\underset{{\mathbf{z}}_k}{argmin}L({\mathbf{y}}^{\tau +1},\mathbf{z};{\sigma }^{\tau })\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill & {\sigma }_k^{\tau +1}:={\sigma }_k^{\tau }+\rho \left({\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k^{\tau +1}\right).\hfill \end{array}$$

(13c)

In order to resolve the problem (13), ADMM needs to be rewritten in a more intuitive form.

(1). The update of consensus variable $\mathbf{y}$:

Since the constant term cannot change the optimization solution of $\mathbf{y}$, removing the constant term yields a precise form:

$${\mathbf{y}}^{\tau +1}:=\underset{\mathbf{y}}{argmin}\sum _{k=1}^K\left({\sigma }_k^{\tau }\mathbf{y}+\frac{\rho }{2}{∥\mathbf{y}-{\mathbf{z}}_k^{\tau }∥}^2\right).$$

(14)

It is easy to see that (14) is a convex quadratic function of $\mathbf{y}$, which can be solved by taking the derivative to 0:

$${\mathbf{y}}^{\tau +1}=\frac{1}{K}\sum _{k=1}^K\left({\mathbf{z}}_k^{\tau }-\frac{1}{\rho }{\sigma }_k^{\tau }\right).$$

(15)

(2). The update of local variable ${\mathbf{z}}_k$:

Similarly, removing the constant term yields a precise form:

$${\mathbf{z}}_k^{\tau +1}:=\underset{{\mathbf{z}}_k}{argmin}\left(-D_k\left({\mathbf{x}}_k;{\mathbf{z}}_k\right)-\langle {\sigma }_k^{\tau },{\mathbf{z}}_k\rangle +\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k∥}^2\right).$$

(16)

Substituting (10) into (16) yields

$$\begin{array}{cc}\hfill {\mathbf{z}}_k^{\tau +1}=& \underset{{\mathbf{z}}_k}{argmin}\left\{-\underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+{\lambda }_{z,k}\sum _{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\mu }_{z,k}\sum _{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)\right\}\right.\hfill \\ \hfill & \left.-\langle {\sigma }_k^{\tau },{\mathbf{z}}_k\rangle +\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k∥}^2\right\}\hfill \\ \hfill =& \underset{{\mathbf{z}}_k}{argmax}\underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+〈\left[\sum _{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right);\sum _{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)\right]+{\sigma }_k^{\tau },{\mathbf{z}}_k〉\right.\hfill \\ \hfill & -\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k∥}^2\}\hfill \end{array}$$

(17)

where we use the denotation ${\mathbf{z}}_k={[{\lambda }_{z,k};{\mu }_{z,k}]}^{\top }$. Since the duality gap between (9) and its dual problem is zero, there is a saddle point $S({\overline{\mathbf{x}}}_k,{\overline{\mathbf{z}}}_k)$ [31] such that

$$\underset{{\mathbf{z}}_k}{argmax}\underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}S({\mathbf{x}}_k,{\mathbf{z}}_k)=S({\overline{\mathbf{x}}}_k,{\overline{\mathbf{z}}}_k)=\underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\underset{{\mathbf{z}}_k}{argmax}S({\mathbf{x}}_k,{\mathbf{z}}_k).$$

Hence, (17) is equivalent to

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^{\tau +1}=& \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\underset{{\mathbf{z}}_k}{argmax}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+〈\left[\sum _{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right);\sum _{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)\right]+{\sigma }_k^{\tau },{\mathbf{z}}_k〉\right.\hfill \\ \hfill & -\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k∥}^2\}.\hfill \end{array}$$

(18)

For any fixed ${\mathbf{x}}_k$, the maximum on the right-hand side of (18) is uniquely attained by the following update step:

$$\begin{array}{cc}\hfill {\mathbf{z}}_k^{\tau +1}=& \left[\begin{array}{c}{\lambda }_{z,k}^{\tau +1}\\ {\mu }_{z,k}^{\tau +1}\end{array}\right]=\left[\begin{array}{c}{\lambda }_y^{\tau +1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\lambda }_{\sigma ,k}^{\tau }\right)\hfill \\ max\left\{0,{\mu }_y^{\tau +1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)+{\mu }_{\sigma ,k}^{\tau }\right)\right\}\hfill \end{array}\right]\hfill \end{array}$$

(19)

where we note that the inequality constraint $h_n\left({\mathbf{x}}_k\right)\le 0$ is equivalent to $max\{0,h_n\left({\mathbf{x}}_k\right)\}$. For the determined ${\mathbf{z}}_k^{\tau +1}$, substituting (19) into (18), one has

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^{\tau +1}=& \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+\rho 〈{\mathbf{z}}_k^{\tau +1}-{\mathbf{y}}^{\tau +1},{\mathbf{z}}_k^{\tau +1}〉-\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}-{\mathbf{z}}_k^{\tau +1}∥}^2\right\}\hfill \\ \hfill =& \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)-\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}∥}^2+\frac{\rho }{2}{∥{\mathbf{z}}_k^{\tau +1}∥}^2\right\}\hfill \end{array}$$

(20)

where the term $\frac{\rho }{2}{∥{\mathbf{y}}^{\tau +1}∥}^2$ in (20) is seen as a constant for ${\mathbf{x}}_k^{\tau +1}$ and cannot change the optimal solution. Hence, ${\mathbf{x}}_k^{\tau +1}$ is a solution to the following minimization problem:

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^{\tau +1}=& \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{argmin}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+\frac{\rho }{2}{∥\left[\begin{array}{c}{\lambda }_y^{\tau +1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\lambda }_{\sigma ,k}^{\tau }\right)\hfill \\ max\left\{0,{\mu }_y^{\tau +1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)+{\mu }_{\sigma ,k}^{\tau }\right)\right\}\hfill \end{array}\right]∥}^2\right\}.\hfill \end{array}$$

(21)

(3). The update of Lagrange multipliers ${\sigma }_k$:

The Lagrange multipliers are locally stored and updated according to (13c) for each area $k\in \{1,\dots ,K\}$, where all variables involved in this step have been calculated and communicated in the update of $\mathbf{y}$ and ${\mathbf{z}}_k$.

After analyzing the intrinsic update procedure of ADMM, it is easy to see that only the update of ${\mathbf{z}}_k$ involves the exchange of information between regions, and once the updated information is obtained from the adjacent regions, it can also be calculated locally. The update of $\mathbf{y}$ needs to separately solve a quadratic programming problem while updating the Lagrangian multiplier ${\sigma }_k$ is the simplest to carry out. Therefore, ADMM can be implemented in a completely distributed manner, requiring only local information exchange without centralized coordination. The distributed algorithm for the OPF problem eliminates the steps of separately determining the boundary coupling equation and coupling information and utilizes the power flow equation to construct the Lagrangian, which reduces the disclosure of boundary information to a certain extent.

The convergence of ADMM is judged by the original residual ${\mathbf{r}}^{\tau }$ and the dual residual ${\mathbf{s}}^{\tau }$ defined as

$$\begin{array}{cc}\hfill {\mathbf{r}}^{\tau }=& {\left[{\mathbf{y}}^{\tau }-{\mathbf{z}}_1^{\tau };\cdots ;{\mathbf{y}}^{\tau }-{\mathbf{z}}_k^{\tau };\cdots ;{\mathbf{y}}^{\tau }-{\mathbf{z}}_K^{\tau }\right]}^{\top },\hfill \\ \hfill {\mathbf{s}}^{\tau }=& -\rho {\left[{\mathbf{z}}_1^{\tau }-{\mathbf{z}}_1^{\tau -1};\cdots ;{\mathbf{z}}_k^{\tau }-{\mathbf{z}}_k^{\tau -1};\cdots ;{\mathbf{z}}_K^{\tau }-{\mathbf{z}}_K^{\tau -1}\right]}^{\top }.\hfill \end{array}$$

When $\parallel {\mathbf{r}}^{\tau }{\parallel }_2\le {\epsilon }_1$ and $\parallel {\mathbf{s}}^{\tau }{\parallel }_2\le {\epsilon }_2$ for given error bounds ${\epsilon }_1$ and ${\epsilon }_2$, the iteration stops and the optimal solutions are obtained.

4. Distributed ADMM for Online OPF

The previous sections investigate the OPF problem in a single gap and its solution; however, the whole model and solving procedure are based on the assumption that the power output of renewable sources remains unchanged such that an OPF solution in a single gap can be obtained. In the actual electric grid, the renewable power sources, especially the photovoltaic and wind turbine, are usually time-varying and fluctuate due to their uncertainty and unpredictability, which deduces a dynamic OPF problem. If we solve the dynamic OPF with a long time scale, the previous approach will lead to a double-loop implementation.

In order to avoid the double-loop implementation, we present an online ADMM to address the dynamic OPF problem with the fluctuated renewable energy and real-time loads. This is a completely distributed algorithm that makes the optimal decision in each time slot and only needs to know the instantaneous system state.

4.1. Online Convex Optimization

Recent studies showed that many online learning algorithms are designed based on online convex optimization tools, whose decision set is a convex set $\mathcal{X}\subseteq {\mathbb{R}}^N$ in N-dimensional Euclidean space, and costs are a series of bounded convex functions defined on the convex set $\mathcal{X}$. Online convex optimization can be regarded as a repeated game. At time t, the player chooses an action ${\mathbf{x}}^t\in \mathcal{X}$, then the adversary reveals a convex loss function $c^t:\mathcal{X}\to {\mathbb{R}}^N$. The player’s goal is to minimize the regret $R_c\left(T\right)$ on the number of time slots T, which is defined as:

$$R_c\left(T\right):=\sum _{t=1}^T{c}^t\left({\mathbf{x}}^t\right)-\underset{\mathbf{x}\in \mathcal{X}}{min}\sum _{t=1}^T{c}^t\left(\mathbf{x}\right).$$

(22)

Under the appropriate hypotheses, well-constructed online learning algorithms can achieve a regret bound that grows sublinearly in T, i.e., ${lim}_{T\to \infty }{R}_c\left(T\right)/T=0$.

4.2. Online OPF Formulation

The dynamic OPF problem based on online convex optimization can be stated as follows. At each time slot $t\in \{1,2,\cdots ,T\}$, on the basis of the realizations observed so far for the renewable sources ${\{P_n^{R,\tau },Q_n^{R,\tau }\}}_{\tau =1}^t$ and the loads ${\{P_n^{D,\tau },Q_n^{D,\tau }\}}_{\tau =1}^t$, find the OPF $\{P_n^{R,t+1},Q_n^{R,t+1}\}$ in the next time slot, so as to minimize the regret. Specifically, the loss function at time slot t can be formed as

$$c^t\left({\mathbf{x}}^t\right):=-\sum _{k=1}^K{D}_k^t({\mathbf{x}}_k^t;{\lambda }^t,{\mu }^t)$$

(23)

where $D_k^t({\mathbf{x}}_k^t;{\lambda }^t,{\mu }^t)$ is an online version (at time slot t) of the Lagrangian $D_k({\mathbf{x}}_k;\lambda ,\mu )$, defined in (10), associated with the OPF problem in a single gap. Since the dual function (11) is the pointwise infimum of a family of affine functions of $(\lambda ,\mu )$, it is concave, even when the original problem (9) is not convex [34]. Hence, the loss function $c^t\left({\mathbf{x}}^t\right)$ is convex, and then the Lagrangian duality of the dynamic OPF problem can be solved by an online learning algorithm based on the online convex optimization technique.

4.3. ADMM for Online OPF

Considering the advantage of the ADMM framework and avoid a double-loop implementation, we adopt the online ADMM approach [35] to solve the following optimization problem

$$\begin{array}{cc}\hfill \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{minimize}& -\sum _{t=1}^T\sum _{k=1}^K{D}_k^t({\mathbf{x}}_k;{\mathbf{z}}_k)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \mathbf{y}-{\mathbf{z}}_k=0\hfill \end{array}$$

(24)

in an online fashion. The augmented Lagrangian is given by

$$L^t(\mathbf{y},\mathbf{z};\sigma ):=\sum _{k=1}^K\left[-D_k^t({\mathbf{x}}_k;{\mathbf{z}}_k)+\langle {\sigma }_k,\left(\mathbf{y}-{\mathbf{z}}_k\right)\rangle +\frac{\rho }{2}{∥\mathbf{y}-{\mathbf{z}}_k∥}^2\right].$$

From which we can obtain the update steps of the online ADMM approach:

$$\begin{array}{cc}\hfill {\mathbf{y}}^{t+1}:=& \underset{\mathbf{y}}{argmin}{L}^t(\mathbf{y},{\mathbf{z}}^t;{\sigma }^t)=\frac{1}{K}\sum _{k=1}^K\left({\mathbf{z}}_k^t-\frac{1}{\rho }{\sigma }_k^t\right)\hfill \end{array}$$

(25a)

$$\begin{array}{cc}\hfill {\mathbf{z}}_k^{t+1}:=& \underset{{\mathbf{z}}_k}{argmin}{L}^t({\mathbf{y}}^{t+1},\mathbf{z};{\sigma }^t)=\left[\begin{array}{c}{\lambda }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\lambda }_{\sigma ,k}^t\right)\hfill \\ max\left\{0,{\mu }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)+{\mu }_{\sigma ,k}^t\right)\right\}\hfill \end{array}\right]\hfill \end{array}$$

(25b)

$$\begin{array}{cc}\hfill {\sigma }_k^{t+1}:=& {\sigma }_k^t+\rho \left({\mathbf{y}}^{t+1}-{\mathbf{z}}_k^{t+1}\right).\hfill \end{array}$$

(25c)

At each time slot t, the online ADMM updates multiplier vectors ${\mathbf{y}}^t$ and ${\mathbf{z}}_k^t$ to solve variable ${\mathbf{x}}_k^{t+1}$ from the following minimization problem in an online fashion:

$$\begin{array}{cc}\hfill & \underset{{\mathbf{x}}_k\in {\mathcal{X}}_k}{minimize}\left\{\sum _{n\in {\mathcal{R}}_k}{f}_n\left({\mathbf{x}}_k\right)+\frac{\rho }{2}{∥\left[\begin{array}{c}{\lambda }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\lambda }_{\sigma ,k}^t\right)\hfill \\ max\left\{0,{\mu }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)+{\mu }_{\sigma ,k}^t\right)\right\}\hfill \end{array}\right]∥}^2\right\}.\hfill \end{array}$$

(26)

The above procedure can be summarized in Algorithm 1. The iteration complexity of the online ADMM has been extensively studied in the literature (see [35] and the references therein). It can be shown that the iterate $\left\{{\mathbf{x}}^t\right\}$ achieves a sublinear regret bound, i.e., $R_c\left(T\right)=\mathcal{O}\left(T\right)$, under some mild conditions on the problem.

Algorithm 1 Online ADMM for OPF.

Initialize:${\mathbf{z}}_k^0$, ${\sigma }_k^0$, $\rho >0$, and $T\in \mathbb{N}$

1: for each time slot $t=1,\cdots ,T$ do 2:  Calculate constraints $g_n\left({\mathbf{x}}_k^t\right)$ & $h_n\left({\mathbf{x}}_k^t\right)$ according to online data $P_n^{D,t}$, $Q_n^{D,t}$, $P_n^{R,t}$, $Q_n^{R,t}$ 3:  Obtain ${\mathbf{y}}^{t+1}$ according to (25a) 4:  for each sub-region $k=1,\cdots ,K$ do 5:     Calculate ${\mathbf{x}}_k^t$ via solving (26) 6:     Update ${\mathbf{z}}_k^{t+1}$ using (25b) 7:     Update ${\sigma }_k^{t+1}$ using (25c) 8:  end for 9: end for

It is worth noting that, at each time slot t, the online data stream will be fed to the online ADMM based on a single-loop iterative implementation, which strives for faster convergence. Recently, some variants of ADMM were studied that exhibit faster convergence rates while requiring only a little change in the computational effort of each iteration [36]. Here we adopt a practical, efficient method, the Peaceman–Rachford Splitting Method (PRSM) [37], to accelerate ADMM. Specifically, by adding an intermediate multiplier and introducing an underdetermined relaxation factor $\xi \in (0,1)$, the PRSM-ADMM can ensure strict convergence, and its convergence rate is also faster than ADMM. Mathematically, the update rule of the intermediate multiplier in PRSM-ADMM is given by

$${\sigma }_k^{t+0.5}:={\sigma }_k^t+\xi \rho \left({\mathbf{y}}^{t+1}-{\mathbf{z}}_k^t\right).$$

(27)

Then, we update ${\mathbf{z}}_k^t$ according to the new rule

$${\mathbf{z}}_k^{t+1}:=\left[\begin{array}{c}{\lambda }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{g}_n\left({\mathbf{x}}_k\right)+{\lambda }_{\sigma ,k}^{t+0.5}\right)\hfill \\ max\left\{0,{\mu }_y^{t+1}+\frac{1}{\rho }\left({\sum }_{n\in {\mathcal{R}}_k}{h}_n\left({\mathbf{x}}_k\right)+{\mu }_{\sigma ,k}^{t+0.5}\right)\right\}\hfill \end{array}\right].$$

(28)

The final update rule of the multiplier in PRSM-ADMM is

$${\sigma }_k^{t+1}:={\sigma }_k^{t+0.5}+\xi \rho \left({\mathbf{y}}^{t+1}-{\mathbf{z}}_k^{t+1}\right).$$

(29)

The whole procedure can be summarized in Algorithm 2.

Algorithm 2 Online PRSM-ADMM for OPF.

Initialize:${\mathbf{z}}_k^0$, ${\sigma }_k^0$, $\rho >0$, $\xi \in (0,1)$ and $T\in \mathbb{N}$

1: for each time slot $t=1,\cdots ,T$ do   2:  Calculate constraints $g_n\left({\mathbf{x}}_k^t\right)$ & $h_n\left({\mathbf{x}}_k^t\right)$ according to online data $P_n^{D,t}$, $Q_n^{D,t}$, $P_n^{R,t}$, $Q_n^{R,t}$   3:  Obtain ${\mathbf{y}}^{t+1}$ according to (25a)   4:  for each sub-region $k=1,\cdots ,K$ do   5:     Calculate ${\mathbf{x}}_k^t$ via solving (26)   6:     Update ${\sigma }_k^{t+0.5}$ using (27)   7:     Update ${\mathbf{z}}_k^{t+1}$ using (28)   8:     Update ${\sigma }_k^{t+1}$ using (29)   9:  end for 10: end for

5. Numerical Tests

In this section, we illustrate the use of the distributed ADMM on static and dynamic case studies to find the optimal power flow under equation and inequality constraints.

5.1. Simulation Setup and Area Partitioning

The distributed ADMM algorithms are tested on the IEEE 30-bus system, which contains 6 motors and 21 loads. The generator and load adopt IEEE 30-bus standard data, but the load adopts a 24-h spinning reserve, and adding random noise corresponds to the solar data. The photovoltaic solar panels used in this article as generators produce electricity in DC, which needs to transform into AC via an inverter for feeding the electric grid. In this experimental environment, the grid-tied system with the photovoltaic does not need battery storage because the grid acts as a power reserve. If the grid goes down, the inverter automatically shuts off and will not feed solar-generated electricity back into the grid. Therefore, the photovoltaic can be safely incorporated on the designated node to perform distributed power flow calculation [38], whose access points should be relatively concentrated and located as close to the reference node as possible so as to reduce the harmonics caused by grid connection and the impact on the grid. For a typical photovoltaic module [39], the max power is 300W; the voltage and current at maximum power points are $50.6$ V and $5.9$ A, respectively, and the open-circuit voltage and short-circuit current are $63.2$ V and $6.5$ A, respectively. The solar data were extracted from an open-source dataset (https://www.elia.be/, accessed on 20 October 2022). Fifteen-minute-sampled solar generation data were collected on 18 July 2020, from 00:00 to 24:00. Considering the power generation capacity and the standard unit value, the data at each time slot were standardized. For a typical DC/AC inverter [39], the input and output voltages are 354 and 230 V, respectively; the power rating is 10 kW; and the amplitude and frequency modulation index are $0.92$ and 100, respectively.

Section 2.3 introduced the area partitioning approach based on spectral clustering. The partitioning is chosen such that each load bus is assigned to the same area in line impedance as its nearest generator bus. After defining the Lagrangian matrix based on the complex admittance of the IEEE 30-bus system, we simply calculate the K-largest eigenvalues (here $K=2,3$) and cluster all of the buses into K regions. Among them, the Optimal (region) Partition by using the spectral clustering-based method proposed in this paper, is represented by $OP\left(K\right)$, and the other partitions based on the electrical distance between buses [40] or the simpler index clustering are represented by $NP\left(K\right)$. Specifically, in [40], the regions are partitioned in the way that each load is assigned to the same region as its closest generator in terms of line impedance. After that, we obtain the coupling parameter, $0.6117$ and $0.7532$ determine the optimal partitioning plan $\left\{\right(1-8,28);(9-27,29-30\left)\right\}$ and $\left\{\right(1-8,25-30);(9-20);(21-24\left)\right\}$ as the $OP\left(2\right)$ and $OP\left(3\right)$ in Figure 1a, and the specific information are shown in

Table 2. Area partitioning cases for the IEEE 30-bus system.

Partitioning Case	Subregions
$NP\left(2\right)$ [40]	$\{1-8,25-30\}$ & $\{9-24\}$
$OP\left(2\right)$ ours.	$\{1-8,28\}$ & $\{9-27,29,30\}$
$NP\left(3\right)$ [40]	$\{1-10\}$ & $\{11-20\}$ & $\{21-30\}$
$OP\left(3\right)$ ours.	$\{1-8,25-30\}$ & $\{9-20\}$ & $\{21-24\}$

.

5.2. Static OPF Simulation

In this simulation part, the maximum number of iterations $M=2000$ and the original residual and dual residual iteration stop criteria are set to ${\epsilon }_1={\epsilon }_2=1\times {10}^{-4}$.

First, we discuss the influence of the dual update step size $\rho$ in (13) on the number of iterations in different area partitioning cases. Figure 2 shows that the optimal dual update step size corresponding to three partition cases may be different. In that, we traverse $\rho \in {[1,100]}_{\mathbb{N}}$ to roughly determine the optimal value ${\rho }^{*}=15$ for the same residual accuracy, which will guarantee that the distributed ADMM can converge at the fastest speed. The figure also shows that $OP\left(3\right)$ is more sensitive to the value of $\rho$, which means the selection of $\rho$ will profoundly affect the convergence performance in this case.

The number of border branches and the number of global variables correspond to the Lagrangian multiplier variable $\mathbf{y}=(\lambda ;\mu )$ in (10), as well as the number of iterations for the given residual accuracy. Obviously, the more variables are involved, the more iterations are needed to achieve the same residual accuracy.

When the dual step size $\rho =15$, the convergence results of the distributed ADMM are discussed. Figure 3a,b show the relative error (expressed as a percentage) of the target value obtained by ADMM at each iteration over a single slot (assuming that the dynamic data of formula (4a) and (4b) are fixed) with respect to the target value obtained by the centralized method (the exact value), expressed as the gap in the target value. A set of comparisons is made between the optimal partition and the conventional partition with the same number of regions K. Since ADMM gradually converges to the optimal value during each iteration, the gap between the target value obtained by the centralized method will be narrowed in the case of convergence. As shown in Figure 3a,b, all four methods have considerable convergence effects, but ADMM with $OP\left(K\right)$ converges faster than that with $NP\left(K\right)$ and achieves a shorter oscillation period. It is worth noting here that there is a clear rise in Figure 3a,b, which is caused by changes in the penalty parameter. At the beginning of the iteration, the subproblems tend to find better local optima with little regard for connections to other subproblems, which leads to a gap oscillation.

The comparison of public information between the method in [40] and ours is shown in

Table 3. Performance comparison among different partitions.

Partition	$\mathit{OP}\left(2\right)$ Ours	$\mathit{NP}\left(2\right)$ [40]	$\mathit{OP}\left(3\right)$ Ours	$\mathit{NP}\left(3\right)$ [40]
Iterations	52	74	120	153
Time (s)	124	1018	145	1241
Public information	8	16	14	40
Consensus variable	15T	31T	25T	76T
Coupling parameter	0.6117	1.0195	0.7532	1.614
Gap	99.57%	99.89%	99.35%	99.22%

, which gives the differences in iteration number and running time for the same accuracy as well as other parameters, etc. Obviously, our method is significantly superior to [40] in convergence time. It is worth noting that $NP\left(2\right)$ may be superior to $OP\left(2\right)$ in the gap factor, but $OP\left(2\right)$ has better accuracy performance for the same iteration number, which has also been verified by Figure 4. It is clear that the amount of public information in our method is significantly less than [40]. Moreover, our method only needs the number of boundary branches between adjacent subsystems to calculate the amount of public information, while [40] additionally needs to consider the partition strategy on the basis of knowing the number of boundary branches.

Note that the convergence of the Lagrangian multiplier is also an important criterion for judging whether the algorithm is optimal. Let $\mathrm{LMC}=log\parallel {\sigma }^{\tau +1}-{\sigma }^{\tau }{\parallel }_{\infty }$, where ${log\parallel ·\parallel }_{\infty }$ represents the logarithm of the infinite norm. Figure 4 depicts the iteration of multiplier $\sigma$ and shows good convergence performance of the algorithm in three cases, especially for $OP\left(2\right)$ and $OP\left(3\right)$, which also verifies the conclusion suggested in Figure 3. In order to solve the single-slot OPF problem, the proposed SDP relaxation scheme incorporated into the distributed ADMM can achieve better convergence and robustness, as shown in Figure 4.

In Figure 5, we evaluate the convergence of the primal residuals of ADMM. We use $\mathrm{PRC}$ to record the convergence process of the primal residual before which the residual termination condition is defined. Like $\mathrm{LMC}$, $OP\left(K\right)$ has a better convergence speed than $NP\left(K\right)$ and reaches the termination condition faster with the same number of partitions.

This section focuses on the performance of the proposed distributed ADMM algorithm on SDP-OPF. After the test of the indicators, the effectiveness of the algorithm can be seen. In the four partitioning schemes, the algorithms have high accuracy and fast convergence performance. The spectral clustering method further shows its effectiveness in improving the system performance under the same number of regions. However, the modified algorithm has good adaptability to dynamic online problems because it bypasses the boundary coupling equation in essence. To better simulate the time-varying DER injected into the power grid, we introduce an online OPF model, which will be further elaborated upon in the next subsection.

5.3. Online OPF Simulation

This subsection tests the performance of the online ADMM, which will monitor the real-time data of solar-power generation and make the best decision due to its intermittency and uncontrollability. We still adopt the IEEE 30-bus system to be injected photovoltaic generators, as shown in Figure 1b, where the set of solar power buses is $\left\{2,3,14,17,24\right\}$, the set of conventional generator buses $\left\{1,2,5,8,11,13\right\}$, and buses $\left\{7,12,19,21,30\right\}$ will be used later. We still use $OP\left(2\right)$ as the area partitioning scheme due to its good performance in the static scene, which is not shown in Figure 1b for considering a more clean picture.

Figure 6 shows the active/reactive load data $P_{n,t}^D$ and $Q_{n,t}^D$, as well as the active/reactive solar power output data $P_{n,t}^R$ and $Q_{n,t}^R$, which can be found on the same website (https://www.elia.be/, accessed on 10 October 2022). Solar-power output is affected by weather and seasonal variations. The angle of the sun to the solar panel changes with the time of day and seasonal variations. Cloudy and rainy days also contribute to less effectiveness of sunlight collection. The sampling period is one day, and the sampling length is 24 h. Generally, the increase in load will naturally affect the output of generators, and various DESs may be our best way to offset the instantaneous power difference.

5.3.1. Performance Comparison

This subsection investigated the performance comparison of the proposed online algorithm with the closest literature [30]. In [30], the authors solved the static DC-OPF problem by ADMM and Lagrangian duality, while our paper addressed the static and dynamic AC-OPF problem. In order to perform the comparison, we develop a quasi-static distributed OPF problem to reduce the difference, in which we need to rewrite the cost function as the minimum economic operation index; then convex relaxation can be carried out. The resulting optimization problem can continue to perform the comparison with our algorithm in distributed OPF calculation.

For coming into use of the offline algorithm [30], at the end of time slot $t-1$, we assume that the solar power generation of time slot t is the same as that of time slot $t-1$ (persistence prediction), the active power output is determined by solving a distributed OPF problem. Specifically, it is to solve the following problems:

$$\begin{array}{cc}\hfill \underset{\mathbf{W}⪰0,{\mathbf{P}}^{G,t},{\mathbf{Q}}^{G,t}}{minimize}& \sum _{n\in {\mathcal{N}}_{\mathcal{G}}}{\alpha }_n{\left(P_n^{G,t}\right)}^2+{\beta }_n\left(P_n^{G,t}\right)+{\gamma }_n\hfill \end{array}$$

(30a)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& P_n^{G,t}+P_n^{R,t-1}-P_n^{D,t-1}=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Re}\left\{(W_{nn}-W_{n{n}^{\prime }})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(30b)

$$\begin{array}{cc}\hfill & Q_n^{G,t}+Q_n^{R,t-1}-Q_n^{D,t-1}=\sum _{n^{\prime }\in {\Omega }_n}\mathrm{Im}\left\{(W_{nn}-W_{n{n}^{\prime }})y_{n{n}^{\prime }}^{*}\right\}\hfill \end{array}$$

(30c)

$$\begin{array}{cc}\hfill & P_n^{min}\le P_n^{G,t}\le P_n^{max}\hfill \end{array}$$

(30d)

$$\begin{array}{cc}\hfill & Q_n^{min}\le Q_n^{G,t}\le Q_n^{max}\hfill \end{array}$$

(30e)

$$\begin{array}{cc}\hfill & {\left(V_n^{min}\right)}^2\le W_{nn}\le {\left(V_n^{max}\right)}^2\hfill \end{array}$$

(30f)

$$\begin{array}{cc}\hfill & W_{nn}+W_{n^{\prime }{n}^{\prime }}-W_{n{n}^{\prime }}-W_{n^{\prime }n}\le {\left(\Delta V_{n{n}^{\prime }}^{max}\right)}^2\hfill \end{array}$$

(30g)

In time slot t, after the actual solar power generation is realized, OPF is solved again. The total cost of the static algorithm is also ${\sum }_{n\in {\mathcal{N}}_{\mathcal{G}}}{\alpha }_n{\left(P_n^{G,t}\right)}^2+{\beta }_n\left(P_n^{G,t}\right)+{\gamma }_n$. Note that this offline scheme needs to solve the distributed OPF problem twice in each time slot, while the online scheme proposed in our paper only solves the OPF problem once in each time slot.

Figure 7 shows the evolution of total cost between the offline [30] and online algorithms for solving the distributed OPF, in which we still use solar power generation in Figure 6. As a whole, it can be seen that the overall online scheme is smoother than the offline scheme. In particular, when the solar power generation capacity suddenly drops or increases, the cost of the online scheme is much lower than the offline scheme. This is because the offline scheme fully trusts the prediction of renewable energy power generation and significantly reduces or increases the conventional power generation that has been invested, while the online scheme in our paper considers the impact of time-scale, which has a better performance in long-range dependence. Therefore, the online scheme can avoid the power output of conventional generators instantaneously surging when the power generation of renewables drops sharply. This conclusion has been confirmed in Figure 8, which describes the conventional power generation submitted by the offline [30] and online schemes in each time slot. It is known from Figure 6 that solar injection basically affects the power grid from the period 04:45–19:30; thus Figure 8 only intercepts the power generation in this period.

5.3.2. Online Performance Analysis

Under the time-varying load and solar power output given in Figure 6, we adopt the distributed online ADMM to solve the real-time optimal power flow. Figure 9 describes the convergence of the optimal value of OPF obtained by the algorithm over time under the background of no solar power access, 5 solar access points and 10 solar access points, as given in Figure 1b. The calculation example starts from 00:00. It is worth noting that within the initial hour, the algorithm needs to reach convergence within a certain period of time, that is, the optimal running state. In Figure 9, an obvious rising phase can be obtained. Obviously, due to the access of the high load side, the power input needs to climb to reach the balance with the load. After that, the algorithm performance reaches the optimum at this time in the configuration without solar energy access, which is the blue curve in Figure 9. Since the objective function is a linear function of active power, it will change with the load, as shown in the load curve in Figure 6. In the configurations with 5 solar access points and 10 solar access points, the power output is obviously constrained by the solar access point.

Finally, we evaluate the convergence performance of our distributed online ADMM for solving the real-time OPF problem. The algorithm makes the optimal decision for each time slot t to understand the instantaneous system state and minimizes the regret given in (22), as shown in Figure 10. In that, two versions of ADMM show very fine distinctions; PRSM-ADMM can approach convergence on a shorter time slot t and may play a greater advantage in solving online problems. It is particularly easy to see the algorithm has an obvious convergence lower bound on the time series, which verifies the good performance of the algorithm in the online case.

As a summary, we discuss the results and performance of the algorithm on online data. The algorithm has real-time calculation results for online data and can quickly reach the lower bound of regret.

6. Conclusions

This paper proposes a distributed control method for the OPF with renewables, which relies on SDP relaxation due to its nonconvexity. In summary, this work is split into two parts. First, based on the Lagrangian duality and ADMM, the method proposed is realized by solving the dual problem of the original OPF, in which the Lagrangian multipliers of the power flow equation are copied to each sub-area instead of the boundary coupling equation. Second, this distributed control method is extended to solving the dynamic OPF with renewables in an online fashion, which determines the current OPF value in each time slot, and achieves a regret bound within a limited time. Finally, the simulation results on the IEEE 30-bus system verify the high efficiency of the proposed distributed control method in solving both the static and dynamic OPF problems, and the online method can effectively avoid the violent fluctuations in conventional generator output coping with renewable’s rapid fluctuation in comparison with the offline algorithms.

Future works will mainly focus on addressing the two limitations of this study. First, since the algorithm of this paper considers the SDP based on AC power flow, it may be an interesting problem to extend to the large-scale AC-DC hybrid power model [41]. Second, for the dynamic, distributed OPF problem, although the experience proved that the SDP relaxation method has good interpretability for large-scale power grids, this method may not be the best match; it becomes an important task to seek new and more effective methods of analysis and design.