Correction: Distributed optimal power flow

[This corrects the article DOI: 10.1371/journal.pone.0251948.].


Introduction
In modern societies, demand for electricity is expected to be satisfied continuously via controllable generation technologies. An event is a situation in which the demand is not fulfilled. Ten-in-one, a widely used reliability criterion for events, means that an event should occur just once in a 10-year span. To meet this standard, system operators schedule the generation portfolio and the grid systems in advance. For example, a day-ahead unit commitment determines the 24-hourly dispatches, along with unit commitment decisions, to meet varying hourly demands. For each hour, the demand profiles are assumed to be constant, which defines the process's steady-state operation. In the absence of an unexpected disturbance, stochastic hourly demand is the unique source of uncertainty in traditional power system operation. Over the last decade or so, renewable energy resources and smart grid technologies have been integrated into systems to improve energy efficiency and reduce greenhouse gas emissions. This integration has introduced uncertainty into the operation of power systems, presenting a new challenge. If high-precision forecasting could be introduced to estimate future energy resources and control demand, existing operations' tools would remain useful, assuming they could be integrated into the expected effective demand (� demand-expected demand reduction-expected renewable energy resources). Unfortunately, even though the precision of forecasting tools has improved, the errors in their long-term forecasts, for a day ahead, for example, are not yet sufficiently small for reliable operation. Frequent decisions are a potential way to accommodate the uncertainty. For example, a day-ahead 24-hour unit commitment (UC) decision is made once in a daily cycle. If the errors in 2-hour ahead forecasts are small enough, then the UC decision with a forecast every 2 hours would still be a reliable tool for the power system's operation. The computational capability to support such decisions plays a key role in this process.
Optimal power flow (OPF) is a backbone in the steady-state operation of power systems. The characteristics of OPF are highly nonlinear and nonconvex. The computational complexity associated with these characteristics makes power flow (PF) analysis a non-deterministic polynomial-time (NP)-hard problem [1]. In most operational practices, a linear approximation of OPF, namely, direct current (DC) OPF, is pursued. Although easy to solve, DC OPF does not address voltage problems, losses, and the dispatch of reactive power generation. Due to these issues, DC OPF may not be feasible. To address the problem correctly, it is ideal to aim for a nonlinear and nonconvex OPF. In addition to the nonconvex nature of the full OPF, uncertainties increase the number of variables in traditional, central decision-making processes. Therefore, frequent but short-term decisions concerning large-scale power systems can be challenging. With the recent advancements in hardware in multi-core machines, distributed computation becomes an attractive approach for enhancing computational efficiency. An exemplary area in power system analysis concerns the use of distributed computation for OPF. Motivated readers can find information related to distributed approaches to solving OPF problems [2,3]. Within distributed computation, the alternating direction method of multipliers The dotted line indicates that the best linear fit for the relationship is N sub / Nb 0.88 . If the optimization problem is solved by SDP, the computation time is found in W N 3 sub À � [21] that yields a computation time that is proportional to ϑ(Nb 2.64 ). In addition, the approaches require significant communication costs as well. Therefore, the overall computation cost is much higher than that of the central optimization. In addition to the computational efficiency, it is not guaranteed that the voltages at the boundary of each subsystem represent the voltages of other nodes inside the same subsystem correctly. If they do not, convergences may not be observed because the information exchange is limited to the boundary buses. The slow convergences and/or non-monotonous convergences reported in previous studies [14][15][16][17] indicate the insufficient representativeness of the boundary buses. A relatively fast convergence is reported by Engelmann et al. [18] in exchange of per-step communication costs by sharing the sensitivities in addition to the local primal variables. The increase in the communication cost is found where n i is the dimension of the local gradient at the i th group. Although the progress at each iteration is faster than those of other ADMM approaches, the communication cost itself is much higher than the total computation costs of the nonconvex heuristic solvers or of the SDP solvers. As a result, although these studies are worth exploring, we conclude that the aggregation approach is not practical in terms of computational efficiency and the inferior quality of the solutions for a scalable algorithm due to the tradeoff issue, nonconvexity, and the modeling problem. A new network model for OPF is necessary for the distributed computation. The contributions of this paper are 1) a new network model that yields direct communication among nodes regardless of the system size, 2) a distributed, fast, and efficient algorithm to solve highly nonlinear and nonconvex OPF problems, and 3) a scalable algorithm that guarantees convergence to a local minimizer. The paper is organized as follows: the theory section proposes a new network model designed for distributed computation and presents an algorithm to solve a full OPF; the next section describes the details of the implementation of the proposed algorithm; the results and discussion present the results for the OPF and comparisons with those from other studies; the following section provides conclusions and research directions for further improving the computational efficiency; and the appendices sketch the proofs of the ranks of the matrices associated with OPF problems and of the convergence.

Proposed network model and algorithm
We propose a new network model for a nodal OPF to keep communication costs manageable regardless of the system size. For this purpose, the desired properties of the model are as follows: 1. The model must be compatible with PF studies for which Kirchhoff's laws and voltage magnitudes are well defined.
2. Each node is a short distance away from the rest of the nodes to minimize the communication costs.
3. The voltages at a node and those at the rest of the nodes are linearly related.
In the power flow studies, such as PF, OPF, state estimation, and probabilistic PF, voltages are the variables. The constraints in the studies consist of the power flows and injections, as well as the voltage magnitudes in terms of voltages. In the Cartesian coordinate system, the power flows and injections are quadratic in voltage. For example, the power flow over the line connecting Nodes i and j at i is f  See S1 Appendix for the proof of the claims. For a real-valued symmetric rank-4 matrix M j , a real-valued eigen pair λ and u exist such h i |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl fflffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } Note that the eigenpairs are dependent solely on the system parameters (not on voltages). With the eigenpairs, the quantities of interest in the PF are the real power injection at Node j p j , the reactive power injection at Node j q j , the real power flow over a line j-l at Node j � f j jÀ l , the reactive power flow over a line j-l at Node j � f j jÀ l , and the squared voltage's magnitude at Node j E j . There are quadratic relationships among them: Eq (1) expresses Kirchhoff's laws and the voltage magnitudes. The nodal variables α j , β j , γ j-l, j , δ j-l,j , and ω j are defined as follows: the α j are the voltages projected onto the eigenspaces spanned by the real power injection, and the β j cover the reactive power injection, the γ j-l,j cover the real power flow over j-l at Node j, the δ j-l,j cover the reactive power flow over j-l at Node j, and the ω j cover the voltage magnitudes at the j th node in the same manner. Although there is one each of α j , β j , and ω j at each node, the nl j of γ j-l,j and δ j-l,j are defined such that nl j represents the number of branches that are directly connected to Node j (i.e., l = 1, 2, . . ., nl j ). It is worth noting that the nodal variables are linearly associated with the voltages (the definitions of the nodal variables in Eq (1)), which indicates that PL = 1. Therefore, the nodal variables satisfy all the desired properties of the new network model. Next, we let nl j be the number of lines connected to the j th node. Then, the dimensions of α, β, γ, δ, and ω are 4, 4, 4nl j , 4nl j , and 2, respectively. Note that their dimensions do not depend on the system size. Let μ j be a nodal variable that integrates local generation variables as follows: where μ end = 1. The generalized PF formulation is listed in previous work [22]. Note that the cardinality of the variable is fixed-independent of the system size. Therefore, the new network model yields: 1) a fixed number of variables, and 2) a PL of 1.

Proposed network model: A star grid with two channels
In the definition of the local variable x j , it is noticeable that there are three variable types: 1) power injection and flow variables α, β, γ, and δ, 2) the voltage magnitude variable ω, and 3) quadratic variables, Q x j comprising power flows � f j jÀ k ; � f j jÀ k � � and power generations (g j ). The first two types of variables are linear, and the third type is quadratic, in terms of the voltages. Once x j (α, β, γ, and δ) is determined, it is straightforward to find the power generations at Node j with Eq (1) (i.e., x j defines the feasible region of μ j ). The voltage v p can be reconstructed from α, β, γ, and δ. The communication path used to collect their values is termed the power channel. Here, ω is directly associated with the voltage magnitudes, and the values are reported through another communication path termed the voltage channel to update voltages v m . Although the voltages v p and v m should be identical, conditions are relaxed so that they can be different. Fig 2 illustrates an example of the traditional network model (left) and the proposed network model (right) for a modified 4-bus system that has branches connecting 1-2, 1-3, 1-4, 2-4, and 3-4. In the proposed model, all the nodes are connected to the center node (black dot at the center) through the power channel (red lines) and the voltage channel (green lines), and all the nodes are of distance 1 from the center (PL = 1). The network is a complete bipartite graph of order 2 with a maximum-diameter star network. All communications between the variables take place in local nodes (α, β, γ, δ, and ω), and the center node is linear in terms of voltages.

Nodal OPF with nodal variables and its convex relaxation
Even though we describe the OPF problem in this paper, the same framework is applied to any other problem, including PF, OPF, state estimation, and probabilistic PF problems, because they have unified formulas [22]. The central OPF formulation as a nonconvex, quadratically constrained quadratic problem is given in Eq (3), its distributed optimization problem (nodal OPF) at each node I is given in Eq (3A), and the nodal OPF is formulated in terms of an ADMM algorithm in Eq (3B).  s:t: real and reactive power balance definitions for the real and the reactive power flows as shown in 2 ð2Þ i � 0 real and reactive power generation limits limits for the voltage magnitudes where A is a block matrix to collect an element and X j is the column space of Ф j . at Node j defined by the constraints listed in 1-5 in (3A).
Note that all A's are full column rank matrices that collect relevant parts from μ i , and that Ф i includes the linear relationships between the local variable x i and the central variable v and y. Using the definition of the function f j , the following observations are made: 1) f j is a nonconvex, but smooth, function that is C 1 on an open set containing X j (defined by the column space of Ф j ), and rf j is Lipschitz continuous on X j ; 2) X j is nonempty, closed, and convex; and 3) w j is coercive.
Even though Eqs (3A) and (3B) have low cardinalities in the decision variables due to the nonconvex nature of the problem, the uniqueness and existence of the solution are not guaranteed. To address the complexity issue, a surrogate function h j approximating Eq (3B) is defined as follows: Because u j is the SDP relaxation of f j , it is convex, Lipchitz continuous, and continuously differentiable on X j . Note that u j relaxes the rank constraint only; therefore, the first derivatives of u j and of f j regarding x j are identical (i.e., r x j f j x j Because the SDP is convex, the solution at the k th iteration is uniquely determined byx kþ1 j ≜ arg min given y k and z j k . Note that the relaxed nodal OPF is convex and has a small number of variables; therefore, its computational complexity remains manageable.

Properties of
The nodal problem f j is C 1 smooth, but nonconvex, and prox-regular at the Descent Lemma [23] with the properties of g j and of h j described below. Its SDP relaxation u j is strongly convex, Lipschitz continuous, and prox-regular, and its first derivative Whereas the variables considered in f j and u j are all local primal variables only, ADMM constraint g j depends on the central variable y and the multipliers z. With the properties of u j , f j , and g j , w j is C 1 smooth, but nonconvex, and prox-regular, whereas h j is strongly convex and Lipschitz continuous. Because all X j are bounded sets, w j and h j are coercive and lower-bounded over X j . Becausẽ x kþ1 j is the optimizer of convex h j (i.e.,x kþ1 j ≜ arg min , the following holds:

Proposed algorithm with provable convergence
Overview of proposed algorithm. Even though the number of variables in the nodal OPF is much smaller than that in the original problem, the complexity of the problem makes it difficult to solve the nodal problem due to its nonconvex nature. To manage the nonconvex nature of the problem, the nonconvex components are relaxed. Once identified, the relaxed solution is mapped onto the feasible region, to the closest point in the feasible region from the identified solution. There are three cases: 1) the convexified problem is not feasible, 2) the convexified problem is feasible, but the solution to the problem is too far from the feasible region, and 3) the solution to the convexified problem is sufficiently close. Only if the mapped point is close enough, the nodal variables are updated and fed to the center node to update voltages.
Proposed algorithm. We propose the following algorithm to solve w j in conjunction with h j in Eq 5: Distributed, Regulated, Optimally-Homogeneous, and Scalable (DROHS) Algorithm 1. Set k = 0, initialize all parameters, such as Δ k 2(0,1), y k , ρ j , 2. If x k j satisfies the termination criteria, terminate the algorithm. 3. Solve the local problems,x kþ1 j ≜ arg min • set the error bound asymptotically to vanish as iteration proceeds (i.e., lim  After the x-optimization is performed, y-optimization follows. The y-optimization at given x kþ1 j and z k is a least square problem-convex problem. It should be noted that the solution to the convex problem is always feasible. Rule 6 is a striking feature of the proposed algorithm that is critically different from the ADMM approaches. Instead of the locally independent update of the primal variables in the ADMM approaches, the nodal variables are updated 1) linearly for the computational efficiency, and 2) with respect to the central variable so that the update is reasonably agreeable among nodal variables. The proof of the convergence of the proposed algorithm is presented in this paper's S2 Appendix.

Implementation of Rule 4
The local SDP relaxation yields a matrix (known asZ To construct an equivalent criterion to kz kþ1 to determine whether or not the feasible projection is acceptable. The Let l m Z k j be the m th eigenvalue ofZ k j . Then, the criterion becomes l 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

Choice of parameters
In contrast to the ADMM approaches, only a few parameters are assigned in the proposed algorithm. Because z k j n o 2 Null F j � � , and the specific choice of z j is irrelevant, all z j 's are random vectors in the null space of known F j . It is important to hold X 1 k¼1 kε k k 2 < 1 to guarantee the convergence of the proposed algorithm [24]. Due to the fact that the maximum τ over all nodes at the k th iteration is set by t k max ¼ t 0 max =k: The parameters used were as follows: ρ j = 20 for v L and 200 for v M ; Δ 0 = 0.3; and the number of maximum iteration = 100; t 0 max ¼ 10 À 3 . The global variable y is randomly assigned (cold start) to test the algorithm's robustness. For comparison, a flat start and a warm start are also attempted. Here, a flat start refers a point at which real voltages are unity, and imaginary voltages are within [-0.1, 0.1]; and a warm start [20] is a solution to the PF. It is recognized that a faster convergence is obtained when the initial point is close to the feasible region, which is commonly observed in numerical iterative methods.

Termination criterion
The algorithm is terminated when no further progress is made. The progress is measured in terms of the global variable y, and the local variables x and z. After the solution is identified, the objective function W � is determined according to (3B) where W x; y; z ð Þ ≜ The interim value of W at the k th iteration is , and the progress measure η k at the k th iteration is defined: The criterion used in Rule 2 terminates the algorithm if the progress measure is less than the tolerance (10 −7 ) (i.e., η k � 10 −7 ).

Simulation environments
The model systems used in the simulations are available from MATPOWER [25] for 3-, 4-, 9-, 14-, 24-, 30-, 39-, 57-, 85-, 118-, 300-, and 2,000-bus systems. All simulations were performed using a Mac pro with 2 × 2.93 GHz 6-core Intel Xeon processors and 6 GB 1333 MHz DDR3 memory. The local SDP problems were solved using the CVX solver [26]. To compare the results from various systems, the lines in the figures are normalized so that all start at zero. The initial points for all the cases are cold starting points, and all the real and imaginary voltages are set to random numbers (not even a flat start). The qualities of the solutions are all numerically identical to those identified using MATPOWER for the model systems tested (all the solutions v � are less than 10 −4 from the MATPOWER solutions (v MATPOWER ), kv � − v MATPOWER k/kv � k � 10 −4 ). For comparison, we also attempted a flat start and observed that N iter decreases at least 30%, but the flat start also finds the same solutions.

Subsystems of the test cases
In the ADMM approaches referred to in the literature [14][15][16][17][18][19][20], the grid partitioning is performed, but the details of the partitioning are not well listed. In considering the computational coupling, the partitioning is based on spectral clustering [20,27] while each subsystem contains at least one generator [14,20]. The inclusion of a generator in a subsystem seems a natural choice to fulfill the power balance equality constraints (1 in (3)) within each substation. However, the inclusion of a generator does not serve the purpose well, because the generator may not be dispatched if its generation costs are excessively high. The spectral clustering is performed using two different algorithms: 1) unnormalized [27] and 2) normalized [28]. The algorithms yield different grid clustering. For example, the unnormalized algorithm results in two clusters for the IEEE 3-bus systems, while the normalized algorithm finds a single cluster for the same system [28].  [14,20]; path length affects the communication costs and N iter because the decision at each subsystem should be delivered to the rest of the system; the number of subsystems affects the number of computation cores and the communication costs unless PL equals one; and the maximum number of nodes affects the computation costs at each subproblem. Therefore, the decrease in the computation costs at each subsystem is achieved in exchange for the increased costs of communication among subsystems. Because the path length and number of subsystems increase rapidly with system size, the improvement in computational efficiency of the ADMM approach is questionable.

Maximum cardinality of nodal variables in the proposed algorithm
The nodal variables in the x-optimization μ j are Because μ end is unity, the cardinality of μ j is 10nl + 2ng + 10 where nl is the number of lines connected to the j th node, and ng is the number of generators located at the node. In determining the voltages through the power channel, the node with a high value for nl plays a key role.

PLOS ONE
relationship is not necessarily positive. The number of variables in the central OPF is typically 2Nb + 2Ng where Nb and Ng are the numbers of buses and generators in the system, respectively. For small systems, such as 3-, 4-, 9-, 14-, and 24-bus systems, the maximum cardinalities of the nodal variables are higher than the number of variables in the central OPF; therefore, it would be challenging to keep the computational costs of the distributed OPF for the small systems lower than those of the central OPF. It is noteworthy that the communication costs remain manageable, because each node directly communicates with the central node due to the fact that PL equals unity.

Distributed OPF for 3-, 4-, 24-, 39-, and 85-bus systems
To examine whether the choice of parameters affects the convergences, optimizations were executed with many randomly assigned initial points, and uniform convergences were consistently observed (Fig 6). In general, N iter depends on the initial guess of the voltages. As the distance of the initial point from the solution becomes closer, N iter decreases, which is commonly observed in numerical iterative methods. All the solutions identified are numerically the same as the solutions using MATPOWER.  ., IEEE 3-, 4-, 9-, 14-, 24-, 39-, and 85-bus systems).
A worst-case convergence was observed for the IEEE 39-bus system. The nodal OPF for the system involves large negative eigenvalues (i.e., the nodal OPFs are highly nonconvex). According to Rule 4, a relatively large number of solutions are rejected because they are not close to the feasible regions and, therefore, N iter becomes large (approximately 40 iterations). However, the solution identified is a local minimizer, and N iter is still much smaller than the number of iterations for the ADMM approaches. For the visual presentations of the results obtained for various systems, the progress is normalized to the 1 st iteration (i.e., the curves begin at 0 for the 1 st iteration).

Comparison of convergence to ADMM approaches
The convergence measures were compared to those from the ADMM approaches reported in multiple studies [14,[17][18][19]. It is worth noting that Erseghe [14] and Zhang et al. [17] do not take the flow limits into consideration; the approach of Engelmann et al. [18] involves a high communication cost; and the approach of Madani, Kalbat, and Lavaei [19] may yield a physically infeasible solution. We attempted the ADMM approaches with the flow constraints and feasibility, but our implemented solver failed to converge with any starting points. Instead, we compared the convergence behaviors of the proposed algorithm to those reported in the previous [14,[17][18][19]. For the visual presentation, the convergence behaviors are normalized so that all the curves begin at the same point (see Fig 7 for the comparison of the convergence).
In the previous research [14,[17][18][19], the size of the subsystems increased with that of the system. Therefore, the tradeoff between communication costs amongst the subsystems and the computation costs for each subsystem make it difficult to develop a scalable algorithm. As the system size increases, N iter increases significantly for the results in two of the previous studies [14,19]. The performance measure in another study [17] fluctuates consistently with various systems, indicating that the convergence may not be guaranteed. Because the final study [18] requires the information exchange regarding sensitivities as well as the primal variables, the computation and communication costs per iteration increase much more rapidly than do the computation costs of a central OPF solver.
Different from the clustering containing multiple nodes used in the ADMM approach, the proposed algorithm for each subsystem contains only 1 node regardless of the system size; each subsystem directly communicates with the rest of the entire system; and the communication involves only the primal variables. With this difference in mind, the fast convergences observed in the proposed algorithm are similar to those in two of the previous studies [17,18], which are consistently much faster than those in the other two previous studies [14,19]. Whereas N iter are like those reported by Engelmann et al. [18], the convergences of the proposed algorithm occur more uniformly and consistently.
Many approaches using ADMM do not consider the flow limits [14,17], and/or feasibility [19] and, therefore, it is not appropriate to compare the quality of the solutions. Instead, we compared our solutions with those using a central OPF solver, MATPOWER. For all the cases described above, they yield numerically identical solutions. It is not clear why the proposed algorithm finds the same solution as the central OPF solver. We tested with the central SDPrelaxation for a small-scale system to estimate the "global" solution. Due to the memory issue, our SDP tests are limited up to 118-bus systems. When the relaxation returns the rank-1 solution, the solution is global and physically feasible. For several cases, the test cases are of the global solutions that are also identified by MATPOWER and by the proposed algorithm. However, there are cases where the SDP relaxation returns a physically infeasible solution. For these cases, both MATPOWER and the proposed algorithm identify numerically the same local minimizers that may not be global solutions. MATPOWER finds a point that meets the  ., IEEE 9-, 14-, 30-, 57-, 118-, and 300-bus systems). The convergences of the ADMM approaches [14,[17][18][19], are compared. https://doi.org/10.1371/journal.pone.0251948.g007

PLOS ONE
Distributed OPF first-order necessary conditions for optimality [25]. Although it finds a minimizer in most cases, there is no guarantee that the identified solution is a minimizer-it can be either a maximizer or a saddle point. On the other hand, the proposed algorithm identifies a solution that meets the second-order optimality conditions, which guarantees that the solution is a minimizer [29]. From a practical point of view, there are two advantages of the proposed algorithms over MATPOWER: 1) they are numerically stable and, therefore, robust because all the subproblems are convex-no issues regarding the rank-deficient Hessian matrix, and 2) they use distributed computation and, therefore, manageable computation in each subproblem. However, there are disadvantages of the proposed algorithm: 1) the requirement of multiple cores to perform the distributed optimization, and 2) the number of nodal variables can be larger than that of the central OPF; for example, 3-, 4-, 9-, 14-, and 24-bus (See Fig 5).

Large-scale OPF: 2,000-bus system
We tested the algorithm for a large-scale system, the 2,000-bus system that is a synthetic grid on a footprint of Texas. Fig 8 presents a convergence pattern. The solution identified is numerically identical to the one found using MATPOWER. N iter remains small, as is the case for the small systems (See Figs 6 and 7). Note that the nodal OPF includes only 1 bus, and the communication costs remain small because PL equals unity. If a sufficient number of computation cores are provided, the proposed algorithm is scalable if the computation cost for solving a subproblem does not increase significantly as the system size increases.
Different from other algorithms reported in the literature (i.e., [14][15][16][17][18]), the proposed algorithm converges regardless of the initial points, but the number of iterations decreases by at least a factor of 2 if it uses a flat start. Another shortcoming of the algorithms reported in previous research [14][15][16][17][18], is the quality of their solutions. Whereas the existing OPFs find low-voltage [17] or suboptimal [16] solutions, the proposed algorithm finds the same solution as a central OPF solver. One previous study [16] claimed that there might be an optimal number of partitioned subsystems due to the tradeoff between communications and computational costs used to obtain a local solution. The proposed method keeps the PL = 1, and the central computation is a simple addition of x j through both the power and voltage channels. From the simulations with various parameter values such as ρ j , Δ k , t 1 j , and a, we obtained the same solutions, indicating that the proposed algorithm is robust as well as efficient. In addition, the proposed model does not require any partitioning.

Comparison to a heuristic central OPF solver, MATPOWER
In comparing the convergence between the proposed algorithm and the central OPF, the number of variables and N iter were examined. The maximum cardinality n max var of μ j increases with the system size in n max var / Nb 1=4 (See Fig 4). Fig 9 illustrates the number of variables in the central OPF as well as the maximum cardinality of the nodal variables. The blue line is the best-fit line for the central OPF (n OPF var ¼ 2Nb þ 2Ng / Nb 0:96 ), whereas the red line is the best-fit line for the distributed algorithm. It is clear that n OPF var increases with Nb in a much faster way than n max var does. The black line is the boundary at which n OPF var equals n max var . If enough computation cores are available, the proposed algorithm involves reduced computation costs per iteration for systems larger than or equal to IEEE 24-bus. Whereas n OPF var almost linearly increases with Nb because Ng � Nb in most systems, n max var does not necessarily increase theoretically. A key observation is that the proposed approach can be a scalable algorithm if N iter does not rapidly increase with the system size. From the comparisons (the number of variables and N iter ), we conclude that the computation cost of the distributed algorithm increases with Nb at a much slower rate than that of the central OPF if sufficient computational resources are available.

Scalability of the algorithm
The computation efficiency of a distributed computation depends on the number of iterations and the computational cost per iteration. The cost is determined by the local computation where CT represents the maximum nodal computation time. The total computation cost is bounded by the product between N iter and CT. The total computation cost is bounded by CT / ϑ(dNb/N core eNb 0.66 ) where N core is the number of cores. If the computational resource is sufficient, CT / ϑ(Nb 0.66 ). From this observation, the computation cost increases sub-linearly for a large-scale network. If a single core is available for computation, the computation cost is in W Nb 1þp SDP =4 ð Þ where p SDP is the computational complexity for a nodal SDP. For the state-of-art central heuristic OPF solvers, the computation cost is in ϑ(Nb 1.5 ). Therefore, when the computation resources are highly limited, the proposed algorithm would still be efficient with a convex problem solver that yields p SDP � 2. A potential improvement in p SDP of the SDP solver is to explore the sparse structure of the matrices [31] or to utilize a commercial solver such as MOSEK. problem without ignoring any constraints, 4) guaranteed convergence to a local minimum, rather than a maximum or saddle point, that meets the first-order necessary conditions for optimality, and 5) a completely distributed algorithm (i.e., a scalable algorithm, which has never been achieved before in the literature). The challenges that the proposed algorithm faces are 1) an increased number of nodal variables that may be higher than that of the variables in the central OPF for a small system, 2) an increased number of iterations when highly connected nodes involve solutions far from the feasible regions, and 3) a prolonged wait time for nodes with low cardinalities. Therefore, the proposed algorithm is an efficient alternative to the central OPF for a large-scale network. Future research directions include the development of 1) a way to accommodate the impact of the rejected solutions in updating the x-variables if the corresponding nodes are highly connected, 2) an efficient computation to solve the nodal SDP, particularly a way to explore the sparse structure of the nodal OPF, and 3) asynchronous distributed optimization for improving the computational efficiency where the scheduling of the distributed computation is identified in terms of a knapsack problem. We also present the proof showing that the surrogate function improves at every iteration and that the iteration converges to a fixed point of the nonconvex OPF. The numerical results exhibit rapid convergence, and the convergence behavior is discussed.