Patent Publication Number: US-10317970-B2

Title: Distributed optimal power flow processes for unbalanced radial distribution networks

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     The present invention claims priority to U.S. Provisional Patent Application Ser. No. 62/150,411 entitled “Distributed Algorithm for Optimal Power Flow Problem on a Radial Network” to Peng et al., filed Apr. 21, 2015. The disclosure of U.S. Provisional Patent Application Ser. No. 62/150,411 is herein incorporated by reference in its entirety. 
    
    
     STATEMENT OF FEDERALLY SPONSORED RESEARCH 
     This invention was made with government support under Grant No. DE-AR0000226 awarded by the Department of Energy. The government has certain rights in the invention. 
    
    
     FIELD OF THE INVENTION 
     The present invention generally relates to optimal power flow and more specifically to a closed form iterative processes for solving for optimal power flow in an unbalanced network. 
     BACKGROUND 
     An incredible amount of infrastructure is relied upon to transport electricity from power stations, where the majority of electricity is currently generated, to individual homes. Power stations can generate electricity in a number of ways including using fossil fuels or using renewable sources of energy such as solar, wind, and hydroelectric sources. Once electricity is generated it travels along transmission lines to substations. Substations typically do not generate electricity, but can change the voltage level of the electricity as well as provide protection to other grid infrastructure during faults and outages. From here, the electricity travels over distribution lines to bring electricity to individual homes. The infrastructure used to transport electricity through the power grid can be viewed as a graph comprised of nodes and lines. The power stations, substations, and any end user can be considered nodes within the graph. Transmission and distribution lines connecting these nodes can be represented by lines. 
     Distributed power generation, electricity generation at the point where it is consumed, is on the rise with the increased use of residential solar panels and is fundamentally changing the path electricity takes to many user&#39;s homes. The term “smart grid” describes a new approach to power distribution which leverages advanced technology to track and manage the distribution of electricity. A smart grid applies upgrades to existing power grid infrastructure including the addition of more renewable energy sources, advanced smart meters that digitally record power usage in real time, and bidirectional energy flow that enables the generation and storage of energy in additional locations along the electrical grid. 
     SUMMARY OF THE INVENTION 
     Node controllers and power distribution networks in accordance with embodiments of the invention enable distributed power control on an unbalanced network. One embodiment includes anode controller, comprising: a network interface; a processor; a memory containing: a distributed power control application; a plurality of node operating parameters describing the operating parameter of a node in an unbalanced network; and a plurality of node operating parameters describing operating parameters for a set of at least one node selected from the group consisting of an ancestor node and at least one child node; wherein the processor is configured by the distributed power control application to: send node operating parameters to nodes in the set of at least one node; receive operating parameters from the nodes in the set of at least one node; calculate a plurality of updated node operating parameters using an iterative process to determine the updated node operating parameters using the node operating parameters that describe the operating parameters of the node, and the operating parameters of the set of at least one node, where each iteration in the iterative process involves evaluation of a subproblem; and adjust the node operating parameters. 
     In a further embodiment, the subproblem further comprises a closed form solution. 
     In another embodiment, the subproblem further comprises an eigen-decomposition. 
     In a still further embodiment, the iterative process further comprises an alternating direction method of multipliers (ADMM) process. 
     In still another embodiment, the ADMM process further comprises an x-update process, wherein the x-update process comprises minimizing an augmented Lagrangian for an augmented Relaxed Optimal Power Flow (ROPF) expression. 
     In a yet further embodiment, the x-update process is subject to the following constraints: 
                 min           H     i   ⁢           ⁢   0       ⁡     (     x     i   ⁢           ⁢   0       )               over           x     i   ⁢           ⁢   0       =     {       v     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     S     i   ⁢           ⁢   0       (   x   )       ,     s     i   ⁢           ⁢   0       (   x   )         }                   s   .   t   .     ⁢                     (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     ∈     𝕊   +                                   (     s     i   ⁢           ⁢   0     ϕ     )       (   x   )       ∈       𝒥   i   ϕ     ⁢           ⁢   ϕ     ∈     Φ   i       ,             min           H     i   ⁢           ⁢   1       ⁡     (     x     i   ⁢           ⁢   1       )               over           x     i   ⁢           ⁢   1       =     {     v     i   ⁢           ⁢   1       (   x   )       }                 s   .   t   .                   v   _     i   ϕ     ≤       (     v     i   ⁢           ⁢   1     ϕϕ     )       (   x   )       ≤         v   _     i   ϕ     ⁢           ⁢   ϕ       ∈     Φ   i       ,               
where H is a hermitian matrix, i is the node, v, l, S, and s are constraints,   is a set of hermitian matrices, and ϕ is a phase of the unbalanced network.
 
     In yet another embodiment, the ADMM process further comprises a y-update process, where the y-update process comprises minimizing an augmented Lagrangian for an augmented Relaxed Optimal Power Flow (ROPF) expression. 
     In a further embodiment again, the y-update process is subject to the following constraints: 
                     ⁢         min           G   i     ⁡     (     {       y   ji     ❘     j   ∈     N   i         }     )               over         {       y   ji     ❘     j   ∈     N   i         }               s   .   t   .               𝒫   i     ⁡     (     v       A   i     ⁢   i       (   y   )       )       =       v   ii     (   y   )       -         z   i     ⁡     (     S   ii     (   y   )       )       H     -       S   ii     (   y   )       ⁢     z   i   H       +       z   i     ⁢     ℓ   ii     (   y   )       ⁢     z   i   H                                     s   ii     (   y   )       =       -   diag     ⁢     (         ∑     i   ∈     C   i                 ⁢           ⁢       𝒫   i     ⁡     (       S   ji     (   y   )       -       z   j     ⁢     ℓ   ji     (   y   )           )         -     S   ii     (   y   )         )         ,                 
where G and N are matrices, i is the node, j is the child node, A i  is the ancestor node, C i  is the set of child nodes, and v, l, S, and s are constraints.
 
     In another embodiment again, the ADMM process further comprises a Lagrange multiplier update process, where the Lagrange multiplier expression comprises a set of Lagrange multipliers. 
     In a further additional embodiment, each Lagrange multiplier in the set of Lagrange multipliers is evaluated by the processor using the following expression:
 
λ k+1 =λ k +ρ( x   k+1   −y   k+1 )
 
where λ is a Lagrange multiplier in the set of Lagrange multipliers, k is a current iteration, k+1 is a next iteration, ρ is a constant, and x−y is a constraint.
 
     In still yet another embodiment, the updated node operating parameters are further calculated using the node operating parameters that describe a set of operating parameters of at least one node selected from the group consisting of an ancestor node of the ancestor node and at least one child node of the at least one child node. 
     In another additional embodiment, the node operating parameters include power injection, voltage, branch current to the ancestor node, and branch power flow to the ancestor node. 
     In a still yet further embodiment, a power distribution network, comprising: one or more centralized computing systems; a communications network; a plurality of node controllers, wherein each node controller in the plurality of node controllers contains: a network interface; a node processor; and a memory containing: a distributed power control application; a plurality of node operating parameters describing the operating parameters of a node in an unbalanced network; and a plurality of node operating parameters describing operating parameters for a set of at least one node selected from the group consisting of an ancestor node and at least one child node; wherein the node processor is configured by the distributed power control application to: send node operating parameters to nodes in the set of at least one node; receive operating parameters from the nodes in the set of at least one node; calculate a plurality of updated node operating parameters using an iterative process to determine the updated node operating parameters using the node operating parameters that describe the operating parameters of the node, and the operating parameters of the nodes in the set of at least one node, where each iteration in the iterative process involves evaluation of a subproblem; and adjust the node operating parameters. 
     In still yet another embodiment, the subproblem further comprises a closed form solution. 
     In a still further embodiment again, the subproblem further comprises an eigen-decomposition. 
     In still another embodiment again, the iterative process is part of a distributed process for achieving Optimal Power Flow (OPF) that is simplified using a convex relaxation. 
     Another further embodiment of the method of the invention includes: the convex relaxation is a semidefinite program (SDP). 
     Still another further embodiment of the method of the invention includes: the node controllers are modeled in the centralized computing system as an unbalanced radial network. 
     In a further embodiment again, the iterative process further comprises an alternating direction method of multipliers (ADMM) process. 
     In a further embodiment again, the ADMM process further comprises an x-update process, wherein the x-update process comprises minimizing an augmented Lagrangian for an augmented Relaxed Optimal Power Flow (ROPF) expression. 
     In another embodiment again, the x-update process is subject to the following constraints: 
                 min         H     i   ⁢           ⁢   0   ⁢     (     x     i   ⁢           ⁢   0       )                 over           x     i   ⁢           ⁢   0       =     {       v     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     S     i   ⁢           ⁢   0       (   x   )       ,     s     i   ⁢           ⁢   0       (   x   )         }                 s   .   t   .                 (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )           ∈     𝕊   +                     (     s     i   ⁢           ⁢   0     ϕ     )       (   x   )       ∈     𝒥   i   ϕ               ϕ   ∈     Φ   i       ,                               min           H     i   ⁢           ⁢   1       ⁡     (     x     i   ⁢           ⁢   1       )               over           x     i   ⁢           ⁢   1       =     {     v     i   ⁢           ⁢   1       (   x   )       }                 s   .   t   .                   v   _     i   ϕ     ≤       (     v     i   ⁢           ⁢   1     ϕϕ     )       (   x   )       ≤         v   _     i   ϕ     ⁢           ⁢   ϕ       ∈     Φ   i       ,               
where H is a hermitian matrix, i is the node, v, l, S, and s are constraints,   is a set of hermitian matrices, and ϕ is a phase of the unbalanced network.
 
     In a further additional embodiment, the ADMM process further comprises a y-update process, where the y-update process comprises minimizing an augmented Lagrangian for an augmented Relaxed Optimal Power Flow (ROPF) expression. 
     In another additional embodiment, the y-update process is subject to the following constraints: 
             min   ⁢           ⁢       G   i     ⁡     (     {       y   ji     ❘     j   ∈     N   i         }     )                   over   ⁢           ⁢     {       y   ji     ❘     j   ∈     N   i         }                     s   .   t   .           ⁢     i       ⁢     (     v       A   i     ⁢   i       (   y   )       )       =       v   ii     (   y   )       -         z   i     ⁡     (     S   ii     (   y   )       )       H     -       S   ii     (   y   )       ⁢     z   i   H       +       z   i     ⁢     ℓ   ii     (   y   )       ⁢     z   i   H                         s   ii     (   y   )       =     -     diag   (         ∑     i   ∈     C   i         ⁢       i     ⁢     (       S   ji     (   v   )       -       z   j     ⁢     ℓ   ji     (   y   )           )         -     S   ii     (   y   )         )         ,         
where G and N are matrices, i is the node, j is the child node, A i  is the ancestor node, C i  is the set of child nodes, and v, l, S, and s are constraints.
 
     In a still yet further embodiment, the ADMM process further comprises a Lagrange multiplier update process, where the Lagrange multiplier expression comprises a set of Lagrange multipliers. 
     In still yet another embodiment, each Lagrange multiplier in the set of Lagrange multipliers is evaluated by the processor using the following expression:
 
λ k+1 =λ k +ρ( x   k+1   −y   k+1 )
 
where λ is a Lagrange multiplier in the set of Lagrange multipliers, k is a current iteration, k+1 is a next iteration, ρ is a constant, and x-y is a constraint.
 
     In a still further embodiment again, the updated node operating parameters are further calculated using the node operating parameters that describe a set of operating parameters of at least one node selected from the group consisting of an ancestor node of the ancestor node and at least one child node of the at least one child node. 
     In still another embodiment again, the node operating parameters include power injection, voltage, branch current to the ancestor node, and branch power flow to the ancestor node. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a diagram conceptually illustrating a power distribution network in accordance with an embodiment of the invention. 
         FIG. 2  is a diagram conceptually illustrating nodes utilizing node controllers connected to a communication network in accordance with an embodiment of the invention. 
         FIG. 3  is a block diagram of a node controller in accordance with an embodiment of the invention. 
         FIG. 4  is diagram illustrating a radial network in accordance with an embodiment of the invention. 
         FIG. 5  is a diagram illustrating a relationship between nodes and operating parameters in a branch flow model in accordance with an embodiment of the invention. 
         FIG. 6  is is a diagram illustrating a relationship between a two node network and operating parameters in accordance with an embodiment of the invention. 
         FIG. 7  is a flow chart illustrating a process to solve for optimal power flow on an unbalanced power distribution network in accordance with an embodiment of the invention. 
         FIG. 8  is a flow chart illustrating a process to solve for optimal power flow using alternating direction method of multipliers in an iterative solution in accordance with an embodiment of the invention. 
         FIG. 9  is a flow chart illustrating a process for a node to perform a distributed power control application in accordance with an embodiment of the invention. 
         FIGS. 10A-10B  are diagrams illustrating a set of nodes before and after an x-update step, respectively, in accordance with an embodiment of the invention. 
         FIG. 11  is a table illustrating multipliers associated with parameters in accordance with an embodiment of the invention. 
         FIG. 12  is pseudocode illustrating a process for initializing distributed power control for an unbalanced network in accordance with an embodiment of the invention. 
         FIG. 13  is pseudocode illustrating a process for distributed power control on an unbalanced network in accordance with an embodiment of the invention. 
         FIG. 14  is a table illustrating statistics of different simulated networks in accordance with an embodiment of the invention. 
         FIGS. 15A-15B  are diagrams illustrating the topologies for a line and a fat tree network respectively, in accordance with an embodiment of the invention. 
         FIG. 16  is a table illustrating statistics of different simulated line and fat tree networks in accordance with an embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION 
     Turning now to the drawings, systems and methods for distributed control of power distribution systems configured as an unbalanced radial networks in accordance with embodiments of the invention are illustrated. Radial networks have a tree topology where each node is connected to a single unique ancestor and a set of unique children and radial networks are commonly utilized in modeling the distribution side of the power grid. In many embodiments, processing nodes are distributed throughout the power distribution network that control power load, distributed power generation, and remote battery storage. In several embodiments, the processing nodes control the operational parameters of aspects of the power distribution network in an effort to achieve what is often referred to as Optimal Power Flow (OPF). Achieving OPF involves optimizing the operation of a power system with respect to one or more objectives. These objectives can include (but are not limited to) minimizing the amount of power lost during the transmission of power to a user, minimizing the cost of generating the power needed for the system, and/or seeking to optimize other general operational constraints. 
     In a number of embodiments, the processing nodes within the power distribution network perform centralized, distributed, or hybrid processes that coordinate the control of the power distribution network. Centralized processes can use a centralized processing unit to calculate optimal power flow of all nodes within the network. Distributed processes can be based upon messages passed between the processing node and its ancestor and/or child nodes within the radial network. Hybrid processes use a combination of centralized and distributed process steps. In several embodiments, individual processing nodes determine the voltage, power injection, current, and/or impedance of a given power load, distributed power generation, or remote battery storage within the power distribution network by performing a closed form calculation using information received by power and/or child nodes. In many embodiments, the specific calculation utilized by the processing nodes is selected based upon a distributed solution for optimal power flow of the power distribution network obtained using alternating direction method of multiplies (ADMM). 
     In various embodiments, an ADMM process can reduce an optimal power flow problem to either a closed form solution or an eigen-decomposition of a hermitian matrix of a small size. Use of an expression that can either be closed form or an eigen-decomposition of a hermitian matrix obtained in the manner described below can be particularly advantageous relative to techniques for performing distributed control of a power distribution network to achieve OPF requiring the use of an iterative optimization process at each processing node. Specifically, performing a single set of calculations to obtain the control parameters as opposed to repeatedly iterating a set of calculations to obtain the control parameters can be significantly more computationally efficient and can enable the power distribution network to adapt to changes with much lower latency. Additional computational efficiencies may be obtained by a conic relaxation of the distribution network, for example by using a second order cone program (SOCP) relaxation. 
     Systems and methods for performing distributed control of unbalanced radial power distribution networks to achieve OPF and solutions to the distributed unbalanced OPF problem that can be utilized in the implementation of such systems and methods in accordance with embodiments of the invention are discussed further below. 
     Radial Power Distribution Networks 
     A power distribution network in accordance with an embodiment of the invention is shown in  FIG. 1 . Electricity is generated in power generator  102 . Power transmission lines  104  can transmit electricity between the power generator and power substation  106 . Power substation  106  additionally can connect to large storage battery  108  which temporarily stores electricity, as well as power distribution lines  110 . The power distribution lines  110  can transmit electricity from the power substation to homes  112 . The homes can include solar panels  114 , a house battery  116 , and/or an electric car  118 . Power distribution networks can transmit electricity in many ways. When alternating current is used, voltage reverses direction at regular intervals. When only one voltage source is used, the power distribution network is described as single phase. When several sources are used, and the sources are distributed in equally spaced regular intervals (typically 120 degrees for a commonly used three phase network), the power distribution network is described as multiphase balanced network. Single phase and multiphase balanced network problems can often be solved with similar analysis. 
     In the discussions to follow, networks that distribute power in a multiphase unbalanced manner will be referred to as multiphase networks. Single phase and multiphase balanced networks will be referred to as single phase networks. 
     The power generator  102  can represent a power source including those using fossil fuels, nuclear, solar, wind, or hydroelectric power. Substation  106  changes the voltage of the electricity for more efficient power distribution. Solar panels  114  are distributed power generation sources, and can generate power to supply the home as well as generate additional power for the power grid. House battery  116  can store excess electricity from the solar panels to power the home when solar energy is unavailable, or store electricity from the power grid to use at a later time. Substations  106 , large storage batteries  108 , homes  112 , solar panels  114 , house batteries  116 , and electric cars  118  can all be considered to be nodes within the power distribution network and the distribution lines  110  can be considered to be lines within the power distribution network. In combination, the nodes and lines form a radial network. In many embodiments, node controllers are located at nodes throughout the network to control the operating parameters of different nodes to achieve OPF. The type of control utilized can depend on the specifics of the network and may include distributed, centralized, and/or hybrid power control. Although many different systems are described above with reference to  FIG. 1 , any of a variety of power distribution network including node controllers may be utilized to perform power distribution as appropriate to the requirements of specific applications in accordance with embodiments of the invention. Nodes utilizing node controllers connected to a communication network in accordance with various embodiments of the invention are discussed further below. 
     Node Controller Architectures 
     Nodes utilizing node controllers connected to a communication network in accordance with an embodiment of the invention are shown in  FIG. 2 . Nodes  202  can connect to communication network  204  using a wired and/or wireless connection  206 . In some embodiments the power distribution network can act in place of the communication network. The communication network may also be connected to one or more centralized computers  208  to monitor calculations made by or to send instructions to multiple nodes to, for example, control power distribution in the network at a global level. Additionally, in many embodiments a database management system  210  can be connected to the network to track node data which, for example, may be used to historically track power usage at various locations over time. Although many systems are described above with reference to  FIG. 2 , any of a variety of systems can be utilized to perform connecting nodes to a network as appropriate to the requirements of specific applications in accordance with embodiments of the invention. Node controllers in accordance with various embodiments of the invention are discussed further below. 
     A node controller in accordance with an embodiment of the invention is shown in  FIG. 3 . In various embodiments, node controller  300  can perform calculations in a distributed manner at a node in the radial network. The node controller includes at least one processor  302 , an I/O interface  304 , and memory  306 . The at least one processor  302 , when configured by software stored in memory, can perform calculations on and make changes to data passing through the I/O interface as well as data stored in memory. In many embodiments, the memory  306  includes software including the distributed power control application  308  as well as operating parameters  310 , operating parameters of ancestor and children nodes  312 , and updated operating parameters  314 . A node can calculate updated operating parameters by using a combination of its own node operating parameters, and/or operating parameters received through the I/O interface from its ancestor and children nodes, and/or operating parameters received from a centralized computer in the case of a centralized or hybrid approach. The distributed power control application  308  will be discussed in greater detail below and can enable the node to perform calculations to solve for optimal power flow in a distributed manner. In many embodiments, the distributed power control application can be applied to unbalanced radial networks. These distributed calculations preformed on the current operating parameters can specifically involve evaluating an OPF problem to solve for either a closed form expression or an eigen-decomposition of a hermitian matrix expression when solving for optimal power flow. Various operating parameters of a node that can be controlled by a node controller are also discussed further below, and may include (but are not limited to) node voltage, current, impedance, and power injection. Although a number of different node controller implementations are described above with reference to  FIG. 3 , any of a variety of computing systems can be utilized to control a node within a power distribution system as appropriate to the requirements of specific applications in accordance with various embodiments of the invention. As noted above, node controllers in accordance with many embodiments of the invention can control the operation of nodes within a radial power distribution network in such a way as to approach OPF. Use of node controllers to implement OPF in a distributed manner within a radial network in accordance with various embodiments of the invention are discussed further below. 
     Use of Node Controllers to Achieve Optimal Power Flow 
     Node controllers in accordance with many embodiments of the invention utilize processes that control nodes in a manner that attempts to achieve OPF in a computationally efficient manner. In order to do this, a closed form expression has been developed enabling calculations to preform at each node concurrently. The term ‘closed form expression’ refers to a calculation that can be performed in a finite number of operations. Overall the closed form solution for OPF can be more computationally efficient and predictable to compute than the use of an iterative process by a node controller to determine operating parameters. Various models can be used to develop a closed form solution that can be utilized to achieve OPF in a distributed manner. 
     The branch flow model (BFM) and the bus injection model (BIM) can be used for solving the OPF problem. The BFM focuses on the current and power in the branches of the model, and the BIM focuses of current, voltage, and power injection at the nodes of the model. Although the BFM and the BIM are generated with different sets of equations and variables, they produce related solutions since they are both modeled based on Kirchhoffs laws. The process utilized by the processing nodes in accordance with various embodiments of the invention utilizes calculations determined by the BFM. Many network shapes can be used to construct the BFM, such as a radial network. In certain cases the structure of a radial network can simplify the computations of the power equations in the OPF problem. Additionally, a convex relaxation of the model can further simplify the calculations. An approach to solve OPF of an unbalanced network in a distributed manner using a semidefinite program (SDP) relaxation is described in detail below. As can readily be appreciated, any of a variety of techniques that can be utilized to solve the OPF in a distributed closed form manner can be utilized as the basis for configuring node controllers as appropriate to the requirements of specific applications in accordance with various embodiments of the invention. Therefore, the inventions described herein should not be considered to be limited to the specific closed form expressions discussed below. 
     Use of Node Controllers to Achieve Optimal Power Flow 
     Node controllers in accordance with many embodiments of the invention utilize processes that control nodes in a manner that attempts to achieve OPF in a computationally efficient manner. In order to do this, a closed form expression has been developed enabling calculations to performed at each node concurrently. The term ‘closed form expression’ refers to a calculation that can be performed in a finite number of operations. Overall the closed form solution for OPF can be more computationally efficient and predictable to compute than the use of an iterative process by a node controller to determine operating parameters. Various models can be used to develop a closed form solution that can be utilized to achieve OPF in a distributed manner. 
     The branch flow model (BFM) and the bus injection model (BIM) can be used for solving the OPF problem. The BFM focuses on the current and power in the branches of the model, and the BIM focuses of current, voltage, and power injection at the nodes of the model. Although the BFM and the BIM are generated with different sets of equations and variables, they produce related solutions since they are both modeled based on Kirchhoffs laws. The process utilized by the processing nodes in accordance with various embodiments of the invention utilizes calculations determined by the BFM. Many network shapes can be used to construct the BFM, such as a radial network. In certain cases the structure of a radial network can simplify the computations of the power equations in the OPF problem. Additionally, a convex relaxation of the model can further simplify the calculations. An approach to solve OPF in a distributed closed form manner using a second order cone program (SOCP) conic relaxation is described in detail below. As can readily be appreciated, any of a variety of techniques that can be utilized to solve the OPF in a distributed closed form manner can be utilized as the basis for configuring node controllers as appropriate to the requirements of specific applications in accordance with various embodiments of the invention. Therefore, the inventions described herein should not be considered to be limited to the specific closed form expressions discussed below. 
     Branch Flow Model 
     A radial network in accordance with an embodiment of the invention is shown in  FIG. 4 . In various embodiments, radial network  400  includes a node  402 . Node  402  has an ancestor node  404  and one or more children nodes  406 . A radial network also has a unique root node  408 . A detailed discussion of these nodes is discussed further below. Although many radial networks are described above with reference to  FIG. 3 , any of a variety of network configurations can be utilized as the network shape as appropriate to the requirements of specific applications in accordance with embodiments of the invention. The relationship between nodes and operation parameters in a BFM is discussed further below. 
     The relationship between nodes and operating parameters in a BFM is discussed in accordance with an embodiment of the invention is shown in  FIG. 5 . A node  502  has a unique ancestor node  504 . Node  502  and unique ancestor node  504  are connected by line  514 . Both node  502  and unique ancestor node  504  have a series of operating parameters. In many embodiments of the invention, example operating parameters for node  502  include voltage  506  and power injection  510 . Unique ancestor node  504  has corresponding voltage  508  and power injection  512 . The line  514  also has operating parameters which for example include an impedance value as well as a current and/or power injection  516 . In various embodiments, node  502  can be connected to one or more child nodes  518  by line  520 . The relationship between node  502  and child node  518  is similar to the relationship between node  502  and ancestor node  504 , where node  502  acts like the ancestor node. The relationship between nodes and operating values will be discussed in greater detail below. Although many operating parameters between two nodes in a BFM are described above with respect to  FIG. 5 , any of a variety of operating parameters can be utilized as operating parameters as appropriate to the requirements of specific applications in accordance with embodiments of the invention. 
     The OPF problem on unbalanced radial distribution networks will be described in greater detail below. In many embodiments OPF can be solved though SDP relaxation. In many embodiments, the set of complex numbers can be denoted with  , the set of set of n-dimensional complex numbers can be denoted with    n , and the set of m×n complex matrix can be denoted with    m×n . The set of hermitian (positive semidefinite) matrix can be denoted by   ( +). The hermitian transpose of a vector (matrix) x can be denoted by x H . The trace of a square matrix x∈   n×n  can be denoted by tr(x):=Σ i=1   n x ii . The inner product of two matrices (vectors) x,y∈   m×n  can be denoted by  x,y :=Re(tr(x H y)). In several embodiments, the Frobenius (Euclidean) norm of a matrix (vector) x∈   m×n  can be defined as ∥x∥ 2 :=√{square root over (x,x)}. Diag(x)∈   n×1  can denote the vector composed of x&#39;s diagonal elements when given x∈   n×n    
     In several embodiments, a distribution network can be modeled by a directed tree graph  :=( ,E) where  :={0, . . . , n} represents the set of buses (nodes) and E represents the set of distribution lines connecting the buses (nodes) in  . The root of the tree can be indexed by 0 and in several embodiments,    + := \{0} can denote the other buses (nodes). For each bus (node) i, there can be a unique ancestor A i , and a set of children buses (nodes) which can be denoted by C i . In various embodiments the graph orientation where every line points towards the root can be adopted. Each directed line connects a bus (node) i and its unique ancestor A i . Hence, the lines can be labeled by ε:={1, . . . , n} where each i∈ε denotes a line from i to A i . Note that ε=   +  and    +  can be used to represent the lines set for convenience. 
     In many embodiments, a, b, c can denote the three phases of the network. For each bus (node) i∈ , Φ i   ⊂  {a,b,c} can denote the set of phases. In a typical distribution network, the set of phases for bus (node) i is a subset of the phases of its parent and superset of the phases of its children. i.e. Φ i    ⊂ Φ A     i    and Φ j    ⊂ Φ i  for j∈C i . On each phase ϕ∈Φ i , V i   ϕ ∈  can denote the complex voltage and s i   ϕ :=p i   ϕ +jq i   ϕ  can denote the complex power injection. In various embodiments, let V i :=(V i   ϕ ,ϕ∈Φ i )∈   |Φ     i     | , s i :=(s i   ϕ ,ϕ∈Φ i )∈   |Φ     i     |  and v i :=V i   H V i ∈   |Φ     i     |×|Φ     i     | . For each line i∈   +  connecting bus (node) i and its ancestor A i , the set of phases can be Φ i  ∩Φ A     i   =Φ i  since Φ i   ⊂ Φ A     i   . On each phase ϕ∈Φ i , let I i   ϕ ∈  denote the complex branch current. In several embodiments, let I i :=(I i   ϕ ,ϕ∈Φ i )∈   |Φ     i     | , l i :=I i I i   H ∈   |Φ     i     |×|Φ     i     |  and S i :=V i I i   H ∈   |Φ     i     |×|Φ     i     | . Two bus (node) network  600  in accordance with an embodiment of the invention is illustrated in  FIG. 6 . In many embodiments, notations can include a variable without a subscript which can denote the set of variables with appropriate components, as summarized below. 
     
       
         
           
               
               
               
             
               
                   
                   
               
             
            
               
                   
                 v := (v i , i ∈    ) 
                 s := (s i , i ∈    ) 
               
               
                   
                 l := (l i , i ∈    + ) 
                 S := (S i , i ∈    + ) 
               
               
                   
                   
               
            
           
         
       
     
     In many embodiments, given a radial network  , the branch flow model for unbalanced network can be defined by: 
                             i     ⁢     (     v     A   i       )       =       v   i     -       z   i     ⁢     S   i   H       -       S   i     ⁢     z   i   H       +       z   i     ⁢     ℓ   i     ⁢     z   i   H                 i   ∈     +                   (     3   ⁢   a     )                       s   i     =     -     diag   (         ∑     i   ∈     C   i         ⁢       i     ⁢     (       S   j     -       z   j     ⁢     ℓ   j         )         -     S   i       )               i   ∈                 (     3   ⁢   b     )                       (           v   i           S   i               S   i   H           ℓ   i           )     ∈     +             i   ∈     +                   (     3   ⁢   c     )                       rank   ⁡     (           v   i           S   i               S   i   H           ℓ   i           )       =   1           i   ∈     +                   (     3   ⁢   d     )               
where    i  (v A     i   ) can denote projecting va, to the set of phases on bus (node) i and    i  (S j −z j l j ) can denote lifting the result of S j −z j l j  to the set of phases Φ i  and filling the missing phase with 0, e.g. if Φ A     i   ={a,b,c}, Φ i ={a,b} and Φ j ={a}, then
 
     
       
         
           
             
               
                 i 
               
               ⁢ 
               
                 ( 
                 
                   v 
                   
                     A 
                     i 
                   
                 
                 ) 
               
             
             := 
             
               ( 
               
                 
                   
                     
                       v 
                       
                         A 
                         i 
                       
                       aa 
                     
                   
                   
                     
                       v 
                       
                         A 
                         i 
                       
                       ab 
                     
                   
                 
                 
                   
                     
                       v 
                       
                         A 
                         i 
                       
                       ba 
                     
                   
                   
                     
                       v 
                       
                         A 
                         i 
                       
                       bb 
                     
                   
                 
               
               ) 
             
           
         
       
       
         
           
             
               
                 i 
               
               ⁢ 
               
                 ( 
                 
                   
                     S 
                     j 
                   
                   - 
                   
                     
                       z 
                       j 
                     
                     ⁢ 
                     
                       ℓ 
                       j 
                     
                   
                 
                 ) 
               
             
             := 
             
               ( 
               
                 
                   
                     
                       
                         S 
                         j 
                         aa 
                       
                       - 
                       
                         
                           z 
                           j 
                           aa 
                         
                         ⁢ 
                         
                           ℓ 
                           j 
                           aa 
                         
                       
                     
                   
                   
                     0 
                   
                 
                 
                   
                     0 
                   
                   
                     0 
                   
                 
               
               ) 
             
           
         
       
     
     In several embodiments, given a vector (v,s,l,S) that satisfies (3), it can be proved that the bus (node) voltages V i  and branch currents I i  can be uniquely determined if the network is a tree. Hence this model (3) can be equivalent to a full unbalanced AC power flow model. 
     OPF and SDP Relaxation 
     In many embodiments, the OPF problem seeks to optimize certain objectives. e.g. total power loss or generation cost, subject to unbalanced power flow equations (3) and various operational constraints. An objective function of the following form can be considered: 
     
       
         
           
             
               
                 
                   
                     F 
                     ⁡ 
                     
                       ( 
                       s 
                       ) 
                     
                   
                   := 
                   
                     
                       ⁢ 
                       
                         
                           f 
                           i 
                         
                         ⁡ 
                         
                           ( 
                           
                             s 
                             i 
                           
                           ) 
                         
                       
                     
                     := 
                     
                       ⁢ 
                       
                         
                           ∑ 
                           
                             ϕ 
                             ∈ 
                             
                               Φ 
                               i 
                             
                           
                         
                         ⁢ 
                         
                           
                             
                               f 
                               i 
                               ϕ 
                             
                             ⁡ 
                             
                               ( 
                               
                                 s 
                                 i 
                                 ϕ 
                               
                               ) 
                             
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     For instance, to minimize total line loss, set ϕ∈Φ i , i∈ ,
 
ƒ i   ϕ ( s   i   ϕ )= p   i   ϕ .  (5)
 
or in other various embodiments, to minimize generation cost, set i∈ ,
 
                         f   i   ϕ     ⁡     (     s   i   ϕ     )       =     (           α   i   ϕ     2     ⁢       (     p   i   ϕ     )     2       +       β   i   ϕ     ⁢     p   i   ϕ         )       ,           (   6   )               
where α i   ϕ ,β i   ϕ &gt;0 depend on the load type on bus (node) i, e.g. α i   ϕ =0 and β i   ϕ =0 for bus (node) i where there is no generator and for generator bus (node) i, the corresponding α i   ϕ , β i   ϕ  depends on the characteristic of the generator.
 
     In many embodiments, for each bus (node) i∈ , there can be two operational constraints on each phase ϕ∈Φ i . First, the power injection s i   ϕ  is constrained to be in an injection region    i   ϕ , i.e.
 
 s   i   ϕ ∈   i   ϕ  for ϕ∈Φ i  and  i∈         (7)
 
The feasible power injection region    i   ϕ  can be determined by the controllable loads attached to phase ϕ on bus (node) i. Second, the voltage magnitude can be maintained within a prescribed region. Note that the diagonal element of v i  describes the voltage magnitude square on each phase ϕ∈Φ i . Thus the constraints can be written as
 
   v     i   ϕ   ≤v   i   ϕϕ   ≤ v     i   ϕ   i∈         (8)
 
where v i   ϕϕ  denotes the ϕ th  diagonal element of v i . Typically the voltage magnitude at substation buses (nodes) is assumed to be fixed at a prescribed value. i.e.  v   0   ϕ = v   0   ϕ  for ϕ∈Φ 0 . At other load buses (nodes) i∈   + , the voltage magnitude is typically allowed to deviate by 5% from its nominal value, i.e.,  v   i   ϕ =0.95 2  and  v   i   ϕ =1.05 2  for ϕ∈Φ i .
 
     To summarize, in many embodiments the OPF problem for unbalanced radial distribution networks can be: 
     
       
         
           
             
               
                 
                   OPF 
                   ⁢ 
                   
                     : 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         min 
                       
                       
                         
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 ϕ 
                                 ∈ 
                                 
                                   Φ 
                                   i 
                                 
                               
                             
                             ⁢ 
                             
                               
                                 f 
                                 i 
                                 ϕ 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   s 
                                   i 
                                   ϕ 
                                 
                                 ) 
                               
                             
                           
                         
                       
                     
                     
                       
                         over 
                       
                       
                         
                           v 
                           , 
                           s 
                           , 
                           S 
                           , 
                           ℓ 
                         
                       
                     
                     
                       
                         
                           s 
                           . 
                           t 
                           . 
                         
                       
                       
                         
                           
                             
                               ( 
                               3 
                               ) 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             and 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ( 
                               7 
                               ) 
                             
                           
                           - 
                           
                             ( 
                             8 
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     The OPF problem (9) can be nonconvex due to the rank constraint (3d). In several embodiments, an SDP relaxation for (9) can be obtained by removing the rank constraint (3d), resulting in a semidefinite program (SDP): 
     
       
         
           
             
               
                 
                   ROPF 
                   ⁢ 
                   
                     : 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         min 
                       
                       
                         
                           ⁢ 
                           
                             
                               ∑ 
                               
                                 ϕ 
                                 ∈ 
                                 
                                   Φ 
                                   i 
                                 
                               
                             
                             ⁢ 
                             
                               
                                 f 
                                 i 
                                 ϕ 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   s 
                                   i 
                                   ϕ 
                                 
                                 ) 
                               
                             
                           
                         
                       
                     
                     
                       
                         over 
                       
                       
                         
                           v 
                           , 
                           s 
                           , 
                           S 
                           , 
                           ℓ 
                         
                       
                     
                     
                       
                         
                           s 
                           . 
                           t 
                           . 
                         
                       
                       
                         
                           
                             ( 
                             
                               3 
                               ⁢ 
                               a 
                             
                             ) 
                           
                           ⁢ 
                           
                               
                           
                           - 
                           
                             
                               ( 
                               
                                 3 
                                 ⁢ 
                                 c 
                               
                               ) 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             and 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ( 
                               7 
                               ) 
                             
                           
                           - 
                           
                             ( 
                             8 
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     In many embodiments, the relaxation ROPF (10) provides a lower bound for the original OPF problem (9) since the original feasible set is enlarged. The relaxation can be called exact if every optimal solution of ROPF satisfies the rank constraint (3d) and hence can also be optimal for the original OPF problem. Distributed processes to solve optimal power flow in an unbalanced network are described below. 
     Solving for OPF 
     An overview of a process for solving for optimal power flow in an unbalanced radial network utilizing a distributed solution is illustrated in  FIG. 7 . The OPF model is initialized  702 . The nonconvex model is relaxed  704  using a semidefinite program relaxation. A distributed solution of a relaxed OPF is calculated  706 . In many embodiments, this distributed solution can be either a closed form expression or an eigen-distribution of a hermitian matrix. Finally, results are generated  708 . A detailed discussion of the process follows. 
     An overview of the ADMM process to solve for OPF, where each distributed calculation involves evaluation of either closed form solutions or an eigen-distribution of a hermitian matrix is illustrated in  FIG. 8 . Nodes are initialized  802  prior to performing  804  the x-update calculation. Pseudocode for a process to initialize nodes similar to node initialization  804  is described below with respect to  FIG. 11 . To perform the x-update calculation, a node will request required operating parameters from its ancestor and children nodes to enable performing this calculation. The remaining calculations do not require communication with neighboring nodes. The y-update calculation is performed  806 , and then the Lagrange multiplier update is calculated  808 . The end condition is checked  810 . If the end condition has not been satisfied, the ADMM process begins performing the x-update calculation again. Otherwise, the process is complete. Although the process of ADMM is described with reference to  FIG. 8 , any of a variety of ADMM variants can be utilized as appropriate to the requirements of specific applications in accordance with embodiments of the invention. A detailed discussion of each of the update steps utilized in ADMM process similar to those described above with reference to  FIG. 8  follows. Obtaining closed form expressions or an eigen-distribution of a hermitian matrix for each of the update steps utilized in the distributed processes involves rewriting the relaxed OPF problem in a manner that lends itself to solving as a series of subproblems. 
     In many embodiments, it can be assumed that the SPD relaxation is exact. A distributed process that solves the ROPF problem follows. A distributed process for a broad class of optimization problems can be developed through alternating direction method of multipliers (ADMM). The process can be applied to the ROPF problem, and it can be shown that the optimization subproblems can be solved efficiently either through closed form solutions or eigen-decomposition of a 6×6 matrix. 
     ADMM 
     ADMM blends the decomposability of dual decomposition with the superior convergence properties of the method of multipliers. It can solve optimization problems of the form: 
                         min     x   ,   y               f   ⁡     (   x   )       +     g   ⁡     (   y   )                   s   .   t   .             x   ∈     x       ,     y   ∈     y                               x   =   y                 (   11   )               
where ƒ(x),g(y) are convex functions and    x ,   y  are convex sets. Let λ denote the Lagrange multiplier for the constraint x=y. Then the augmented Lagrangian is defined as
 
                         L   ρ     ⁡     (     x   ,   y   ,   λ     )       :=       f   ⁡     (   x   )       +     g   ⁡     (   y   )       +     〈     λ   ,     x   -   y       〉     +       ρ   2     ⁢            x   -   y          2   2           ,           (   12   )               
where ρ≥0 is a constant. When ρ=0, the augmented Lagrangian degenerates to the standard Lagrangian. At each iteration k, ADMM consists of the iterations:
 
                     x     k   +   1       ∈     arg   ⁢   ⁢       L   ρ     ⁡     (     x   ,     y   k     ,     λ   k       )                 (     13   ⁢   a     )                 y     k   +   1       ∈     arg   ⁢   ⁢       L   ρ     ⁡     (       x     k   +   1       ,   y   ,     λ   k       )                 (     13   ⁢   b     )                 λ     k   +   1       =       λ   k     +       ρ   ⁡     (       x     k   +   1       -     y     k   +   1         )       .               (     13   ⁢   c     )               
Specifically, at each iteration, ADMM first updates x based on (13a), then updates y based on (13b), and after that updates the multiplier λ based on (13c). Compared to dual decomposition. ADMM is guaranteed to converge to an optimal solution under less restrictive conditions. Let
 
 r   k   :=∥x   k   −y   k ∥ 2   (14a)
 
 s   k   :=ρ∥y   k   −y   k−1 ∥ 2 ,  (14b)
 
which can be viewed as the residuals for primal and dual feasibility, respectively. They converge to 0 at optimality and are usually used as metrics of convergence in the experiment.
 
     Referring back to  FIG. 8 , the process performed by node controllers in accordance with many embodiments of the invention involves performing the x-update calculation  804 , the z-update calculation  806 , and the Lagrange multiplier update  808  until an end condition is reached  810 . As noted above, the calculations performed in each of the x-update, y-update, and Lagrange multiplier update calculations involve the evaluation of closed form expressions or eigen-decomposition of a 6×6 matrix. The information exchange that occurs during the update calculations is illustrated in  FIGS. 9, 10A and 10B . A node requests and makes local copies  902  of current operating parameters from its ancestor and children nodes. Updated operating parameters are calculated  904  based upon current operating parameters of ancestor and children nodes. After this calculation is performed, the updated operating parameters are sent  906  to ancestor and children nodes. The end condition is checked  908  and if it has not been met the process begins again by requesting and making local  902  copies of current operating parameters from ancestor and children nodes. Otherwise the process is complete and it ends. The specific operating parameters passed between nodes and their ancestors and/or children nodes before and after the x-update process is illustrated in  FIGS. 10A and 10B . As can readily be appreciated, additional communication and/or information can occur between nodes during the execution of the distributed processes by the node controllers. 
     The below discussion provides a detailed explanation of the manner in which closed form solutions or eigen-decomposition of a 6×6 matrix solution can be developed for subproblems that can be solved by individual node controllers to achieve a distributed solution to achieve OPF in a radial power distribution network. The manner in which the solutions outlined in the above discussion can be utilized to implement processes that are executed in node controllers in accordance with various embodiments of the invention is discussed further below. 
     In many embodiments, the above ADMM process can be generalized such that the optimization subproblems can be solved efficiently for the ROPF problem. Instead of using the quadratic penalty term ρ/2∥x−y∥ 2   2  in (12), a more general quadratic penalty term can be used: ∥x−y∥ Λ   2 , where ∥x−y∥ Λ   2 :=(x−y) H Λ(x−y) and Λ is a positive diagonal matrix. In many embodiments the augmented Lagrangian becomes 
                       L   ρ     ⁡     (     x   ,   y   ,   λ     )       :=       f   ⁡     (   x   )       +     g   ⁡     (   y   )       +     〈     λ   ,     x   -   y       〉     +       ρ   2     ⁢              x   -   y          Λ   2     .                 (   15   )               
The convergence result can carry over directly to this general case.
 
Distributed ADMM Processes
 
     A distributed ADMM process for a broad class of optimization problems can be developed, of which the ROPF problem can be a special case. Consider the following optimization problem: 
                       min         ⁢       f   i     ⁡     (     x   i     )                     (     16   ⁢   a     )                   over         {       x   i     ❘     i   ∈       }                 (     16   ⁢   b     )                     s   .   t   .               ∑     j   ∈     N   i         ⁢       A   ij     ⁢     x   i         =   0             for   ⁢           ⁢   i     ∈                 (     16   ⁢   c     )                                         ⁢       x   i     ∈       ⋂     r   =   0       R   i       ⁢     ir                     for   ⁢           ⁢   i     ∈     ,                 (     16   ⁢   d     )               
where for each i∈ , x i  can be a complex vector, ƒ i (x i ) can be a convex function,    ir  can be a convex set, and A ij  (j∈N i , i∈ ) can be matrices with appropriate dimensions. In many embodiments, a broad class of graphical optimization problems (including ROPF) can be formulated as (16). Each node i∈  can be associated with some local variables stacked as x i , which belongs to an intersection of R i +1 local feasible sets    ir  and has a cost objective function ƒ i (x i ). Variables in node i can be coupled with variables from their neighbor nodes in N i  through linear constraints (16c). The objective can then be to solve a minimal total cost across all the nodes.
 
     In various embodiments, a distributed process that solves (16) can be developed such that each node i solves its own subproblem and only exchanges information with its neighbor nodes N i . In order to transform (16) into the form of ADMM (11), two sets of variables x and y can be used. In many embodiments, two sets of slack variables can be utilized:
         1. x ir  can represent a copy of the original variable x i  for 1≤r≤R i . For convenience, denote the original x i  by x i0 .   2. y ij  can represent the variables in node i observed at node j, for j∈N i .
 
In various embodiments, (16) can be reformulated as
       

                       min         ⁢       f   i     ⁡     (     x     i   ⁢           ⁢   0       )                 over         x   =     {         x   ir     ❘     0   ≤   r   ≤     R   i         ,     i   ∈       }                             y   =     {         y   ij     ❘     j   ∈     N   i         ,     i   ∈       }                   (     17   ⁢   a     )                     s   .   t   .               ∑     j   ∈     N   i         ⁢       A   ij     ⁢     y   ji         =       0   ⁢           ⁢   for   ⁢           ⁢   i     ∈                   (     17   ⁢   b     )                                   x   ir     ∈         ir     ⁢           ⁢   for   ⁢           ⁢   0     ≤   r   ≤       R   i     ⁢           ⁢   i       ∈                 (     17   ⁢   c     )                                   x   ir     ∈         y   ii     ⁢           ⁢   for   ⁢           ⁢   1     ≤   r   ≤       R   i     ⁢           ⁢   i       ∈                 (     17   ⁢   d     )                                     x     i   ⁢           ⁢   0       =         y   ij     ⁢           ⁢   for   ⁢           ⁢   j     ∈       N   i     ⁢           ⁢   i     ∈       ,                 (     17   ⁢   e     )               
where x and y represent the two groups of variables in ADMM. Note that the consensus constraints (17d) and (17e) can force all the duplicates x ir  and y ij  to be the same. Thus its solution x i0  can also be optimal to the original problem (16). (17) can fall into the general ADMM form (11), where (17b) corresponds to    y , (17c) corresponds to    x , and (17d) and (17e) are the consensus constraints that relates x and y.
 
     Following the ADMM procedure, the consensus constraints can be relaxed (17d) and (17e), whose Lagrangian multipliers are denoted by λ ir  and μ ij , respectively. The generalized augmented Lagrangian then can be written as 
                       L   ρ     ⁡     (     x   ,   y   ,   λ   ,   μ     )       =     ⁢       (         ∑     r   =   1       R   i       ⁢     (       〈       λ   ir     ,       x   ir     -     y   ii         〉     +       ρ   2     ⁢              x   ir     -     y   ii              Λ   ir     2         )       +       f   i     ⁡     (     x     i   ⁢           ⁢   0       )       +       ∑     j   ∈     N   i         ⁢     (       〈       μ   ij     ,       x     i   ⁢           ⁢   0       -     y   ij         )     +       ρ   2     ⁢              x     i   ⁢           ⁢   0       -     y   ij              M   ij     2         )         )     .               (   18   )               
where the parameter Λ ir  and M ij  depend on the process describe further below.
 
     In several embodiments, it can be shown that both the x-update (13a) and y-update (13b) can be solved in a distributed manner, i.e. both of them can be decomposed into local subproblems that can be solved in parallel by each node i with only neighborhood communications. 
     A set of local variables can be defined for each node i, denoted by    i , which includes its own duplicates x ir  and the associated multiplier λ ir  for 0≤r≤R i , and the “observations” y ji  of variables from its neighbor N i  and the associated multiplier μ ji , i.e.
 
   i   :={x   ir ,λ ir |0≤ r≤R   i   }∪{y   ji ,μ ji   |j∈N   i }.  (19)
 
It can be shown how each node i can update {x ir |0≤r≤R i } in the x-update and {y ji |j∈N i } in the y-update.
 
     In the x-update at each iteration k, the optimization subproblem that updates x k+1  can be 
                       min     x   ∈     𝒦   x         ⁢       L   ρ     ⁡     (     x   ,     y   k     ,     λ   k     ,     μ   k       )         ,           (   20   )               
where the constraint    x  is the Cartesian product of    ir , i.e.
 
   x :=   i∈       r=0   R     i       ir .
 
The objective can be written as a sum of local objectives as shown below
 
                   L   ρ     ⁡     (     x   ,     y   k     ,     λ   k     ,     μ   k       )       =         ∑     i   ∈   𝒩       ⁢     (         ∑     r   =   1       R   i       ⁢     (       〈       λ   ir   k     ,       x   ir     -     y   ii   k         〉     +       ρ   2     ⁢              x   ir     -     y   ii   k              Λ     ir   ⁢               2         )       +       f   i     ⁡     (     x     i   ⁢           ⁢   0       )       +       ∑     j   ∈     N   i         ⁢     (       〈       μ   ij   k     ,       x     i   ⁢           ⁢   0       -     y   ij   k         〉     +       ρ   2     ⁢              x     i   ⁢           ⁢   0       -     y   ij   k              M   ij     2         )         )       =         ∑     i   ∈   𝒩       ⁢       ∑     r   =   0       R   i       ⁢       H   ir     ⁡     (     x   ir     )           -       ∑     i   ∈   𝒩       ⁢     (         ∑     r   =   0       R   i       ⁢     〈       λ   ir   k     ,     y   ii   k       〉       +       ∑     j   ∈     𝒩   i         ⁢     〈       μ   ij   k     ,     y   ij   k       〉         )             ,         
where the last term is independent of x and
 
     
       
         
           
             
               
                 
                   
                     
                       H 
                       ir 
                     
                     ⁡ 
                     
                       ( 
                       
                         x 
                         ir 
                       
                       ) 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               
                                 
                                   f 
                                   i 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     x 
                                     
                                       i 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       0 
                                     
                                   
                                   ) 
                                 
                               
                               + 
                               
                                 
                                   ∑ 
                                   
                                     j 
                                     ∈ 
                                     
                                       N 
                                       i 
                                     
                                   
                                 
                                 ⁢ 
                                 
                                   ( 
                                   
                                     
                                       〈 
                                       
                                         
                                           μ 
                                           ij 
                                           k 
                                         
                                         , 
                                         
                                           x 
                                           
                                             i 
                                             ⁢ 
                                             
                                                 
                                             
                                             ⁢ 
                                             0 
                                           
                                         
                                       
                                       〉 
                                     
                                     + 
                                     
                                       
                                         ρ 
                                         2 
                                       
                                       ⁢ 
                                       
                                         
                                            
                                           
                                             
                                               x 
                                               
                                                 i 
                                                 ⁢ 
                                                 
                                                     
                                                 
                                                 ⁢ 
                                                 0 
                                               
                                             
                                             - 
                                             
                                               y 
                                               ij 
                                               k 
                                             
                                           
                                            
                                         
                                         
                                           M 
                                           ij 
                                         
                                         2 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                           
                             
                               r 
                               = 
                               0 
                             
                           
                         
                         
                           
                             
                               
                                 〈 
                                 
                                   
                                     λ 
                                     ir 
                                     k 
                                   
                                   , 
                                   
                                     x 
                                     ir 
                                   
                                 
                                 〉 
                               
                               + 
                               
                                 
                                   ρ 
                                   2 
                                 
                                 ⁢ 
                                 
                                   
                                      
                                     
                                       
                                         x 
                                         ir 
                                       
                                       - 
                                       
                                         y 
                                         ii 
                                         k 
                                       
                                     
                                      
                                   
                                   
                                     Λ 
                                     ir 
                                   
                                   2 
                                 
                               
                             
                           
                           
                             
                               r 
                               &gt; 
                               0 
                             
                           
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     In many embodiments, the problem (20) in the x-update can be written explicitly as 
                     min   ⁢       ∑     i   ∈   𝒩       ⁢       ∑     r   =   0       R   i       ⁢       H   ir     ⁡     (     x   ir     )             ⁢     
     ⁢       over   ⁢           ⁢   x     =     {         x   ir     |     0   ≤   r   ≤     R   i         ,     i   ∈   𝒩       }       ⁢     
     ⁢         s   .   t   .           ⁢     x   ir       ∈         𝒦   ir     ⁢           ⁢   for   ⁢           ⁢   0     ≤   r   ≤     R   i         ,     i   ∈   𝒩     ,             (   22   )               
where both the objective and constraint are separable for 0≤r≤R i  and i∈ . Thus it can be decomposed into Σ i∈   (R i +1) independent problems that can be solved in parallel. There are R i +1 problems associated with each node i and the r th  (0≤r≤R i ) one can be simply written as
 
                     min       x   ir     ∈     𝒦   ir         ⁢       H   ir     ⁡     (     x   ir     )               (   23   )               
whose solution is the new update of variables x ir  for node i. In the above problem, the constants y ij   k ,μ ij   k ∈   j  are not local to i and stored in i&#39;s neighbors j∈N i . Therefore, each node i needs to collect (y ij ,μ ij ) from all of its neighbors prior to solving (23). The message exchanges are illustrated in  FIG. 10A .
 
     In the y-update, the optimization problem that updates y k+1  can be 
                     min     y   ∈     𝒦   y         ⁢       L   ρ     ⁡     (       x     k   +   1       ,   y   ,     λ   k     ,     μ   k       )               (   24   )               
where the constraint set    y  can be represented as a Cartesian product of | | disjoint sets, i.e.
 
     
       
         
           
             
               𝒦 
               y 
             
             := 
             
               
                 ⊗ 
                 
                   i 
                   ∈ 
                   𝒩 
                 
               
               ⁢ 
               
                 
                   { 
                   
                     
                       y 
                       ji 
                     
                     , 
                     
                       
                         
                           j 
                           ∈ 
                           
                             N 
                             i 
                           
                         
                         | 
                         
                           
                             ∑ 
                             
                               j 
                               ∈ 
                               
                                 N 
                                 i 
                               
                             
                           
                           ⁢ 
                           
                             
                               A 
                               ij 
                             
                             ⁢ 
                             
                               y 
                               ji 
                             
                           
                         
                       
                       = 
                       0 
                     
                   
                   } 
                 
                 . 
               
             
           
         
       
     
     The objective can be written as a sum of local objectives as below. 
                   L   ρ     ⁡     (     x   ,     y   k     ,     λ   k     ,     μ   k       )       =         ∑     i   ∈   𝒩       ⁢     (         ∑     r   =   1       R   i       ⁢     (       〈       λ   ir   k     ,       x   ir     k   +   1       -     y   ii         〉     +       ρ   2     ⁢              x   ir     k   +   1       -     y   ii              Λ     ir   ⁢               2         )       +       f   i     ⁡     (     x     i   ⁢           ⁢   0       k   +   1       )       +       ∑     j   ∈     N   i         ⁢     (       〈       μ   ji   k     ,       x     j   ⁢           ⁢   0       k   +   1       -     y   ji         〉     +       ρ   2     ⁢              x     j   ⁢           ⁢   0       k   +   1       -     y   ji              M   ji     2         )         )       =         ∑     i   ∈   𝒩       ⁢       G   i     ⁡     (     {       y   ji     |     j   ∈     N   i         }     )         +       ∑     i   ∈   𝒩       ⁢     (         f   i     ⁡     (     x     i   ⁢           ⁢   0       k   +   1       )       +       ∑     r   =   0       R   i       ⁢     〈       λ   ir   k     ,     y   ir     k   +   1         〉       +       ∑     j   ∈     N   i         ⁢     〈       μ   ji   k     ,     y     j   ⁢           ⁢   0     k       〉         )             ,         
where the last term is independent of y and
 
                 G   i     ⁡     (     {       y   ji     |     j   ∈     N   i         }     )       =         ∑     r   =   1       R   i       ⁢     (       -     〈       λ   ir   k     ,     y   ii       〉       +       ρ   2     ⁢              x   ir     k   +   1       -     y   ii              Λ   ir     2         )       +       ∑     j   ∈     N   i         ⁢       (       -     〈       μ   ji   k     ,     y   ji       〉       +       ρ   2     ⁢              x     j   ⁢           ⁢   0       k   +   1       -     y   ji              M   ji     2         )     .               
In various embodiments, the problem (24) in the y-update can be written explicitly as
 
                     min   ⁢       ∑     i   ∈   𝒩       ⁢       G   i     ⁡     (     {       y   ji     |     j   ∈     N   i         }     )           ⁢     
     ⁢       over   ⁢           ⁢   y     =     {       {       y   ji     |     j   ∈     N   i         }     |     i   ∈   𝒩       }       ⁢     
     ⁢         s   .   t   .           ⁢       ∑     j   ∈     N   i         ⁢       A   ij     ⁢     y   ji           =   0     ,     i   ∈   𝒩                               
which can be decomposed into | | subproblems and the subproblem for node i is
 
                     min   ⁢           ⁢       G   i     ⁡     (     {       y   ji     |     j   ∈     N   i         }     )         ⁢     
     ⁢       over   ⁢           ⁢   y     =     {       y   ji     |     j   ∈     N   i         }       ⁢     
     ⁢         s   .   t   .           ⁢       ∑     j   ∈     N   i         ⁢       A   ij     ⁢     y   ji           =   0     ,             (   25   )               
whose solution is the new update of {y ji |j∈N i }∈   i . In (25), the constants x j0 ∈   j  are stored in i&#39;s neighbor j∈N i . Hence, each node i needs to collect x j0  from all of its neighbor prior to solving (25). The message exchanges in the y-update is illustrated in  FIG. 10B .
 
     The problem (25) can be solved with closed form solution. The real and imaginary part of the variables {y ji |j∈N i } can be stacked in a vector with appropriate dimensions and can denoted as {tilde over (y)}. Then (25) can take the following form: 
                       min   ⁢     1   2     ⁢       y   ~     T     ⁢   M   ⁢     y   ~       +       c   T     ⁢     y   ~         ⁢     
     ⁢     over   ⁢           ⁢     y   ~       ⁢     
     ⁢           s   .   t   .           ⁢     A   ~       ⁢           ⁢     y   ~       =   0     ,             (   26   )               
where M is a positive diagonal matrix, Ã is a full row rank real matrix, and c is a real vector. M,c,A are derived from (25). There exists a closed form expression for (26) given by
 
{tilde over ( y )}=( M   −1   Ã   T ( ÃM   −1   Ã   T ) −1   ÃM   −1   −M   −1 ) c.   (27)
 
     In summary, the original problem (16) is decomposed into local subproblems that can be solved in a distributed manner using ADMM. At each iteration, each node i solves (23) in the x-update and (26) in the y-update. There exists a closed form solution to the subproblem (26) in the y-update as shown in (27), and hence whether the original problem (16) can be solved efficiently in a distributed manner depends on the existence of efficient solutions to the subproblems (23) in the x-update, which depends on the realization of both the objectives ƒ i (x i ) and the constraint sets    ir . 
     Next, it can be shown the ROPF problem (10) is a special case of (16), and hence can be solved in a distributed manner using the above method. In particular, it can be shown that the corresponding subproblems in the x-update can be solved efficiently. 
     ADMM Applied to the OPF Problem 
     In various embodiments, the SDP relaxation can be assumed to be exact. A distributed algorithm for solving ROPF (10) can be derived. Using the ADMM based process developed above, the global ROPF problem can be decomposed into local subproblems that can be solved in a distributed manner with only neighborhood communication. Note that the subproblems in the y-update for each node i can always be solved with closed form solution, and therefore an efficient solution for the for the subproblems (23) in the x-update for the ROPF problem is needed. In particular, a sufficient condition can be provided for the existence of efficient solutions to all the optimization subproblems. Compared with existing methods that use generic iterative optimization solver to solve each subproblem, the computation time is improved by more than 100 times. 
     The ROPF problem defined in (10) can be written explicitly as 
                     min   ⁢       ∑     i   ∈   𝒩       ⁢       ∑     ϕ   ∈     Φ   i         ⁢       f   i   ϕ     ⁡     (     s   i   ϕ     )             ⁢     
     ⁢       over   ⁢           ⁢   v     ,   s   ,   S   ,   ℓ             (     28   ⁢   a     )                 s   .   t   .           ⁢       𝒫   i     ⁡     (     v     A   i       )         =         v   i     -       z   i     ⁢     S   i   H       -       S   i     ⁢     z   i   H       +       z   i     ⁢     ℓ   i     ⁢     z   i   H     ⁢   i       ∈   𝒩             (     28   ⁢   b     )                 s   i     =         -     diag   (         ∑     i   ∈     C   i         ⁢       𝒫   i     ⁡     (       S   j     -       z   j     ⁢     ℓ   j         )         -     S   i       )       ⁢   i     ∈   𝒩             (     28   ⁢   c     )                 (           v   i           S   i               S   i   H           ℓ   i           )     ∈       𝕊   +     ⁢   i     ∈   𝒩           (     28   ⁢   d     )                   s   i   ϕ     ∈       𝒥   i   ϕ     ⁢   ϕ     ∈     Φ   i       ,     i   ∈   N             (     28   ⁢   e     )                         v   _     i   ϕ     ≤     v   i   ϕϕ     ≤         v   _     i   ϕ     ⁢   ϕ       ∈     Φ   i       ,     i   ∈   𝒩       ⁢     
     ⁢   Denote           (     28   ⁢   f     )                 x   i     :=     {       v   i     ,     s   i     ,     S   i     ,     ℓ   i       }             (   29   )                 𝒦     i   ⁢           ⁢   0       :=     {         x   i     |       (           v   i           S   i               S   i   H           ℓ   i           )     ∈     𝕊   +         ,     {         s   i   ϕ     ∈     𝒥   i   ϕ       |     ϕ   ∈     Φ   i         }       }             (   30   )                 𝒦     i   ⁢           ⁢   1       :=     {         x   i     |         v   _     i   ϕ     ≤     v   i   ϕϕ     ≤       v   _     i   ϕ         ,     ϕ   ∈     Φ   i         }             (   31   )               
Then (28) can take the form of (16) with R i =1 for all i∈ , where (28b)-(28c) correspond to (16c) and (28d)-(28f) correspond to (16d). A sufficient condition for the existence of an efficient solution to (23) can be proven.
 
     Recall that there is always a closed form solution to the optimization subproblem (25) in the y-update, if the objective function ƒ i   ϕ  (s ϕ ) and injection region    i   ϕ  satisfy the sufficient conditions, and all the subproblems can be solved efficiently. 
     Following the procedure described above, two sets of slack variables can be introduced: x ir  and y ij . Then the counterpart of (17) can be 
                     min   ⁢       ∑     i   ∈   𝒩       ⁢       ∑     ϕ   ∈     Φ   i         ⁢       f   i   ϕ     ⁡     (       (     s     i   ⁢           ⁢   0     ϕ     )       (   x   )       )             ⁢     
     ⁢       over   ⁢           ⁢   x     :=     {         x   ir     |     0   ≤   r   ≤   1       ,     i   ∈   𝒩       }       ⁢     
     ⁢     y   :=     {         y   ji     |     j   ∈     N   i         ,     i   ∈   𝒩       }               (     32   ⁢   a     )                 s   .   t   .           ⁢       𝒫   i     ⁡     (     v       A   i     ⁢   i       (   y   )       )         =         v   ii     (   y   )       -         z   i     ⁡     (     S   ii     (   y   )       )       H     -       S   ii     (   y   )       ⁢     z   i   H       +       z   i     ⁢     ℓ   ii     (   y   )       ⁢     z   i   H     ⁢   i       ∈   𝒩             (     32   ⁢   b     )                 s   ii     (   y   )       =         -     diag   (         ∑     i   ∈     𝒞   i         ⁢       𝒫   i     ⁡     (       S   ji     (   y   )       -       z   j     ⁢     ℓ   ji     (   y   )           )         -     S   ii     (   y   )         )       ⁢   i     ∈   𝒩             (     32   ⁢   c     )                 (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     ∈       𝕊   +     ⁢   i     ∈   𝒩           (     32   ⁢   d     )                   (     s     i   ⁢           ⁢   0     ϕ     )       (   x   )       ∈       𝒥   i   ϕ     ⁢   ϕ     ∈       Φ   i     ⁢           ⁢   and   ⁢           ⁢   i     ∈   𝒩           (     32   ⁢   e     )                     v   _     i   ϕ     ≤       (     v     i   ⁢           ⁢   1     ϕϕ     )       (   x   )       ≤       v   _     i   ϕ       ⁢     
     ⁢     ϕ   ∈       Φ   i     ⁢           ⁢   and   ⁢           ⁢   i     ∈   𝒩             (     32   ⁢   f     )                     x   ir     -     y   ii       =   0     ⁢           ⁢     
     ⁢     r   =       1   ⁢           ⁢   and   ⁢           ⁢   i     ∈   𝒩               (     32   ⁢   g     )                     x     i   ⁢           ⁢   0       -     y   ij       =   0     ⁢     
     ⁢       j   ∈       N   i     ⁢           ⁢   and   ⁢           ⁢   i     ∈   𝒩     ,             (     32   ⁢   h     )               
where superscript (⋅) (x)  and (⋅) (y)  can be put on each variable to denote whether the variable is updated in the x-update or y-update step. Note that each node i does not need full information of its neighbor. In many embodiments each node i, only voltage information v A     i     i   (y)  is needed from its parent A i  and branch power S ji   (y)  and current l ji   (y)  information from its children j∈C i  based on (32). Thus, y ij  can contains only partial information about x i0 , i.e.
 
               y   ij     :=     {             (       S   ii     (   y   )       ,     ℓ   ii     (   y   )       ,     v   ii     (   y   )       ,     s   ii     (   y   )         )           j   =   i               (       S     iA   i       (   y   )       ,     ℓ     iA     i   ⁢                 (   y   )         )           j   =     A   i                 (     v   ij     (   y   )       )           j   ∈     C     i   ⁢                       .             
On the other hand, only x i0  needs to hold all the variables and it suffices for x i1  to only have a duplicate of v i , i.e.
 
               x   ir     :=     {             (       S     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     v     i   ⁢           ⁢   0       (   x   )       ,     s     i   ⁢           ⁢   0       (   x   )         )           r   =   0               (     v     i   ⁢           ⁢   1       (   x   )       )           r   =   i           .             
As a result, x ir , y ii  in (32g) and x i0 , y ij  in (32h) may not consist of the same components. Notations can be abused in both (32g) and (32h), which are composed of components that appear in both items, i.e.
 
     
       
         
           
             
               
                 x 
                 
                   i 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   0 
                 
               
               - 
               
                 y 
                 ij 
               
             
             := 
             
               { 
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 
                                   S 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   S 
                                   ii 
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                               , 
                               
                                 
                                   ℓ 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   ℓ 
                                   ii 
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                               , 
                               
                                 
                                   v 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   v 
                                   ii 
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                               , 
                               
                                 
                                   s 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   s 
                                   ii 
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                             
                           
                         
                         
                           
                             j 
                             = 
                             i 
                           
                         
                       
                       
                         
                           
                             ( 
                             
                               
                                 
                                   S 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   S 
                                   iA 
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                               , 
                               
                                 
                                   ℓ 
                                   
                                     i 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     0 
                                   
                                   
                                     ( 
                                     x 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   ℓ 
                                   
                                     iA 
                                     i 
                                   
                                   
                                     ( 
                                     y 
                                     ) 
                                   
                                 
                               
                             
                             ) 
                           
                         
                         
                           
                             j 
                             = 
                             
                               A 
                               i 
                             
                           
                         
                       
                       
                         
                           
                             ( 
                             
                               
                                 v 
                                 
                                   i 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   0 
                                 
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                               - 
                               
                                 v 
                                 ij 
                                 
                                   ( 
                                   y 
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         
                           
                             j 
                             ∈ 
                             
                               C 
                               i 
                             
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                       x 
                       ir 
                     
                   
                   - 
                   
                     y 
                     ii 
                   
                 
                 := 
                 
                   { 
                   
                     
                       
                         
                           
                             ( 
                             
                               
                                 v 
                                 
                                   i 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   1 
                                 
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                               - 
                               
                                 v 
                                 ii 
                                 
                                   ( 
                                   y 
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         
                           
                             r 
                             = 
                             1 
                           
                         
                       
                     
                     . 
                   
                 
               
             
           
         
       
     
     Let λ denote the Lagrangian multiplier for (32g) and μ the Lagrangian multiplier for (32g). The detailed mapping between constraints and those multipliers are illustrated in  FIG. 11 . Next, the efficient solution for the subproblems in the x-update can be derived. For notational convenience, the iteration number k on the variables will be skipped. In the x-update, there are 2 subproblems (23) associated with each bus (node) i. The first problem, which updates x i0 , can be written explicitly as: 
                       min           H     i   ⁢           ⁢   0       ⁡     (     x     i   ⁢           ⁢   0       )                   (     33   ⁢   a     )                   over           x     i   ⁢           ⁢   0       =     {       v     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     S     i   ⁢           ⁢   0       (   x   )       ,     s     i   ⁢           ⁢   0       (   x   )         }                   (     33   ⁢   b     )                     s   .   t   .             (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     ∈     +                   (     33   ⁢   c     )                         (     s     i   ⁢           ⁢   0     ϕ     )       (   x   )       ∈     𝒥   i   ϕ               ϕ   ∈     Φ   i       ,                 (     33   ⁢   d     )               
where H i0 (x i0 ) is defined in (21) and for our application. ∥x i0 −y ij ∥ M     ij     2  is chosen to be
 
                              x     i   ⁢           ⁢   0       -     y   ij              M   ij     2     =     {                     (       2   ⁢          C   i            +   3     )     ⁢              S     i   ⁢           ⁢   0       (   x   )       -     S   ii     (   y   )              2   2       +              s     i   ⁢           ⁢   0       (   x   )       -     s   ii     (   y   )              2   2     +                 2   ⁢              v     i   ⁢           ⁢   0       (   x   )       -     v   ii     (   y   )              2   2       +       (            C   i          +   1     )     ⁢              ℓ     i   ⁢           ⁢   0       (   x   )       -     ℓ   ii     (   y   )              2   2                     j   =   i                          S     i   ⁢           ⁢   0       (   x   )       -     S     iA   j       (   y   )              2   2     +              ℓ     i   ,     A   i         (   x   )       -     ℓ   i     (   y   )              2   2             j   =   Ai                        x     i   ⁢           ⁢   0       -     y   ij            2   2           j   ∈     C   i                       (   34   )               
By using (34), H i   (1) (S i0   (x) , l i0   (x) , v i0   (x) ), which is defined below, can be written as the Euclidean distance of two Hermitian matrix. In many embodiments this can be a key reason that can lead to an efficient process.
 
     Therefore, H i0 (x i0 ) can be further decomposed as 
                               H     i   ⁢           ⁢   0       ⁡     (     x     i   ⁢           ⁢   0       )       =       ⁢         f   i     ⁡     (     x     i   ⁢           ⁢   0       )       +       ∑     j   ∈     N   i         ⁢     (       〈       μ   ij     ,     x     i   ⁢           ⁢   0         〉     +       ρ   2     ⁢              x     i   ⁢           ⁢   0       -     y   ij              M   ij     2         )                       =       ⁢           ρ   ⁡     (            C   i          +   2     )       2     ⁢       H   i     (   1   )       ⁡     (       S     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     v     i   ⁢           ⁢   0       (   x   )         )         +       H   i     (   2   )       ⁡     (     s     i   ⁢           ⁢   0       (   x   )       )       +   constant       ,           ⁢           
     ⁢           ⁢       where   ⁢     
     ⁢           ⁢       H   i     (   1   )       ⁡     (       S     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     v     i   ⁢           ⁢   0       (   x   )         )         =                  (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     -     (             v   ^     i             S   ^     i                 S   ^     i   H             ℓ   ^     i           )            2   2     ⁢     
     ⁢           ⁢       H   i     (   2   )       ⁡     (     s     i   ⁢           ⁢   0       (   x   )       )         =         f   i     ⁡     (     s     i   ⁢           ⁢   0       (   x   )       )       +       ρ   2     ⁢                s     i   ⁢           ⁢   0       (   x   )       -       s   ^     i            2   2     .                       (   35   )               
In various embodiments, the last step in (35) can be obtained using square completion and the variables labeled with hat are some constants.
 
     Hence, the objective (33a) in (42) can be decomposed into two parts, where the first part H (1) (S i0   (x) , l i0   (x) , v i0   (x) ) involves variables (S i0   (x) , l i0   (x) , v i0   (x) ) and the second part H (2) (s i0   (x) ) involves s i0   (x) . Note that the constraint (42a)-(42a) can also be separated into two independent constraints. Variables (S i0   (x) , l i0   (x) , v i0   (x) ) only depend on (42a) and s i0   (x)  depends on (42a). Then (42) can be decomposed into two subproblems, where the first one (36) solves the optimal (S i0   (x) , l i0   (x) , v i0   (x) ) and the second one (37) solves the optimal s i0   (x) . The first subproblem can be written explicitly as 
                       min           H   i     (   1   )       ⁡     (       S     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     v     i   ⁢           ⁢   0       (   x   )         )               over           S     i   ⁢           ⁢   0       (   x   )       ,     ℓ     i   ⁢           ⁢   0       (   x   )       ,     v     i   ⁢           ⁢   0       (   x   )                   s   .   t   .               (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     ∈     +       ,                 (   36   )               
which can be solved using eigen-decomposition of a 6×6 matrix. In various embodiments:
 
             W   :=     (             v   ^     i             S   ^     i                 S   ^     i   H             ℓ   ^     i           )               and             X   :=       (           v     i   ⁢           ⁢   0       (   x   )             S     i   ⁢           ⁢   0       (   x   )                   (     S     i   ⁢           ⁢   0       (   x   )       )     H           ℓ     i   ⁢           ⁢   0       (   x   )             )     .           
Then (36) can be written abbreviately as
 
               min   X     ⁢            X   -   W          2   2                   s   .   t   .     
     ⁢   X     ∈     +           
which can be solved efficiently using eigen-decomposition. The second problem is
 
                       min             f   i     ⁡     (     s     i   ⁢           ⁢   0       (   x   )       )       +       ρ   2     ⁢              s     i   ⁢           ⁢   0       (   x   )       -       s   ^     i            2   2                             over           s     i   ⁢           ⁢   0       (   x   )       ∈     𝒥   i   ϕ             ϕ   ∈       Φ   i     .                   (   37   )               
Recall that if ƒ i (s i0   (x) )=Σ ϕ∈Φ     i   ƒ i   ϕ ((s i0   ϕ ) (x) ), then both the objective and constraint are separable for each phase ϕ∈Φ i . Therefore, (37) can be further decomposed into |Φ i | number of subproblems as below.
 
min ƒ i   ϕ (( s   i0   ϕ ) (x) )
 
over ( s   i0   ϕ ) (x) ∈   i   ϕ ,  (38)
 
which takes the same form as a closed form expression and thus can be solved with closed form solution based on the assumptions.
 
     For the problem (43) that updates x i1 , which consists of only one component v i1   (x) , it can be written explicitly as
 
min  H   i1 ( x   i1 )
 
over  x   i1   ={v   i1   (x) }
 
 s.t.  v     i   ϕ ≤( v   i1   ϕϕ ) (x)   ≤ v     i   ϕ  ϕ∈Φ i ,  (39)
 
where H i1 (x i1 ) is defined in (21) and for many embodiments of the present invention, ∥x ir −y ii ∥ Λ     ir     2  is chosen to be
 
∥ x   ir   −y   ii ∥ Λ     ir     2   =∥x   ir   −y   ii ∥ 2   2  
 
Then the closed form solution is given as:
 
                 (     v     i   ⁢           ⁢   1         ϕ   1     ⁢     ϕ   2         )       (   x   )       =     {               [         λ     i   ⁢           ⁢   1         ϕ   1     ⁢     ϕ   2         ρ     +       (     v   ii       ϕ   1     ⁢     ϕ   2         )       (   y   )         ]         v   _     i     ϕ   1           v   _     i     ϕ   1                 ϕ   1     =     ϕ   2                     λ     i   ⁢           ⁢   1         ϕ   1     ⁢     ϕ   2         ρ     +       (     v   ii       ϕ   1     ⁢     ϕ   2         )       (   y   )                 ϕ   1     ≠     ϕ   2             .             
To summarize, the subproblems in the x-update for each bus (node) i can be solved either through a closed form solution or a eigen-decomposition of a 6×6 matrix.
 
     In the y-update, the subproblem solved by each node i takes the form of (25) and can be written explicitly as 
                           min   ⁢           ⁢       G   i     ⁡     (     {       y   ij     ❘     j   ∈     N   i         }     )                   over   ⁢           ⁢     {       y   ij     ❘     j   ∈     N   i         }                     s   .   t   .           ⁢     i       ⁢     (     v       A   i     ⁢   i       (   y   )       )       =       v   ii     (   y   )       -         z   i     ⁡     (     S   ii     (   y   )       )       H     -       S   ii     (   y   )       ⁢     z   i   H       +       z   i     ⁢     ℓ   ii     (   y   )       ⁢     z   i   H                 ⁢     
     ⁢         s   ii     (   y   )       =     -     diag   (         ∑     i   ∈     C   i         ⁢       i     ⁢     (       S   ji     (   y   )       -       z   j     ⁢     ℓ   ii     (   y   )           )         -     S   ii     (   y   )         )         ,             (   40   )               
which has a closed form solution given in (27).
 
     Finally, the initialization and stopping criteria for the process can be specified. A good initialization usually reduces the number of iterations for convergence. In various embodiments, the following initialization can be used. The auxiliary variables {V i |i∈ } and {I i |i∈E} can first be initialized, which represent the complex nodal voltage and branch current, respectively. Then these auxiliary variables can be used to initialize the variables in (32). Intuitively, the above initialization procedure can be interpreted as finding a solution assuming zero impedance on all the lines. The procedure is formally stated in in  FIG. 12  which illustrates pseudocode for initialization of a process in accordance with an embodiment of the invention. 
     For the stopping criteria, there is no general rule for ADMM based processes and stopping criteria usually hinge on the particular problem. It can be argued that a reasonable stopping criteria is that both the primal residual r k  defined in (14a) and the dual residual s k  defined in (14b) are below 10 −4 √{square root over (||)}. Various embodiments of the present invention adopt criteria and the empirical results show that the solution is accurate enough. The pseudocode for a complete process in accordance with many embodiments of the invention are illustrated in  FIG. 13 . It should readily be apparent that many other stopping criteria for ADMM based processes can be used. 
     Simulated Results 
     The scalability of the distributed process can be demonstrated by testing it on the standard IEEE test feeders. To show the efficiency of an implementation of the process in accordance with an embodiment of the invention, the computation time of solving the subproblems in both the x and y-update with off-the-shelf solver (CVX) can also be compared. A simulated process in accordance with several embodiments of the invention can be run on networks of different topology to understand the factors that affect the convergence rate. The simulated process is implemented in Python and run on a Macbook pro 2014 with i5 dual core processor. 
     Simulations on IEEE Test Feeders 
     The IEEE 13, 34, 37, 123 bus (node) distribution systems can test a simulated process. In a simulation of the present invention, all the networks are unbalanced three phase networks. The substation is modeled as a fixed voltage bus (node) (1 p.u.) with infinite power injection capability. The other buses (nodes) are modeled as load buses (nodes) whose voltage magnitude at each phase can vary within [0.95, 1.05] p.u. and power injections are specified in the test feeder. There is no controllable device in the original IEEE test feeders, and hence the OPF problem degenerates to a power flow problem. To demonstrate the effectiveness of the process, all the capacitors can be replaced with inverters, whose reactive power injection ranges from 0 to the maximum ratings specified by the original capacitors. The objective is to minimize power loss across the network, namely ƒ i   ϕ (s i   ϕ )=p i   ϕ  for ϕ∈Φ i  and i∈ . 
     In various simulations, the focus can be the time of convergence (ToC) for the distributed process. Furthermore, the process is simulated on a single machine. To roughly estimate the ToC (excluding communication overhead) if the process is run on multiple machines in a distributed manner, the total time by can be divided by the number of buses (nodes). 
       FIG. 14  illustrates the simulated IEEE 13, 34, 37, 123 bus (node) distribution systems number of iterations to converge, total computation time required to run on a single machine and the average time required for each node if the process is run on multiple machines excluding communication overhead. From the simulation results, a process in accordance with many embodiments of the invention converges within 2.5 second for all the standard IEEE test networks if the process is run in a distributed manner. 
     Moreover, the advantage of using a process in accordance with an embodiment of the invention can be shown by comparing the computation time of solving the subproblems between off-the-shelf solver (CVX) and a simulated process in accordance with various embodiments of the invention. In particular, the average computation time of solving the subproblem in both the x and y update can be compared. In the x-update, the average time required to solve the subproblem (40) is 9.8×10 −5  s for a simulated process but 0.13 s for CVX. In the y-update, the average time required to solve the subproblems (42)-(43) are 3.7×10 −3  s for our algorithm but 0.45 s for CVX. Thus, each ADMM iteration takes about 3.8×10 −3  s for a simulated process in accordance with many embodiments of the invention but 5.8×10 −1  s for using an iterative algorithm, a more than 100× speedup. 
     Simulated Network Topology 
     A distributed process in accordance with an embodiment of the invention can dramatically reduce the computation time within each iteration as described above. The time of convergence (ToC) can be determined by both the computation time required within each iteration and the number of iterations. The number of iterations, namely the rate of convergence can be simulated. 
     Rate of convergence can be determined by many different factors. In several embodiments of the present invention, the rate of convergence from two factors, network size N, and diameter D can be considered, i.e. given the termination criteria similar to pseudocode for a distributed OPF process as illustrated in  FIG. 13  above, the impact of network size and diameter on the number of iterations. It should readily be appreciated, the impact from other factors, e.g. form of objective function and constraints can additionally be simulated. 
     To illustrate the impact of network size N and diameter D on the rate of convergence, the process of two extreme cases can be simulated: 1) Line network as illustrated in  FIG. 15A , whose diameter is the largest given the network size, and 2) Fat tree network as illustrated in  FIG. 15B , whose diameter is the smallest given the network size.  FIG. 16  illustrates simulated results from processes using a line network and a fat tree network including the number of iterations for both line and fat tree network of different sizes. For the line network, the number of iterations increases notably as the size increases. For the fat tree network, the trend is less obvious compared to line network. In many embodiments of the invention, network diameter has a stronger impact than the network size on the rate of convergence. 
     Although the present invention has been described in certain specific aspects, many additional modifications and variations would be apparent to those skilled in the art. It is therefore to be understood that the present invention can be practiced otherwise than specifically described without departing from the scope and spirit of the present invention. Thus, embodiments of the present invention should be considered in all respects as illustrative and not restrictive. Accordingly, the scope of the invention should be determined not by the embodiments illustrated, but by the appended claims and their equivalents.