Distributed gradient descent for solving optimal power flow in radial networks

Abstract

Node controllers and power distribution networks in accordance with embodiments of the invention enable distributed power control. One embodiment includes a node controller comprising a memory containing: a plurality of node operating parameters; and a plurality of node operating parameters describing operating parameters for a set of at least one node selected from the group consisting of at least one downstream node and at least one upstream node; wherein the processor is configured by the node controller application to: receive and store in memory a plurality of coordinator parameters describing operating parameters of a node coordinator by the network interface; and calculate updated node operating parameters using an iterative gradient projection process to determine updated node parameters using node operating parameters that describe operating parameters of node and operating parameters of the set of at least one node, where each iteration is determined by the coordinator parameters.

Description

FIELD OF THE INVENTION

The present invention generally relates to optimal power flow and more specifically relates to distributed gradient projection processes for solving for optimal power flow.

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 users' 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. One embodiment includes a node controller comprising a network interface; a processor; a memory containing: a node controller application; a plurality of node operating parameters describing the operating parameters of a node; and a plurality of node operating parameters describing operating parameters for a set of at least one node selected from the group consisting of at least one downstream node and at least one upstream node; wherein the processor is configured by the node controller application to: receive and store in memory a plurality of coordinator parameters describing the operating parameters of a node coordinator by the network interface; and calculate a plurality of updated node operating parameters using an iterative gradient projection process to determine the updated node 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 is determined by the coordinator parameters.

In a further embodiment, the iterative gradient projection process is a distributed process.

In another embodiment, the iterative gradient projection process further comprises a backwards forwards sweep process.

In a still further embodiment, the iterative gradient projection process further comprises calculating gradient parameters.

In still another embodiment, the iterative gradient projection process further comprises calculating a gradient step size.

In a yet further embodiment, the backward forward sweep process further comprises a backward sweep process and a forward sweep process.

In yet another embodiment, the backward sweep process further comprises: receiving operating parameters from the one or more downstream nodes; calculating a plurality of updated node operating parameters using the operating parameters from the one or more downstream nodes; and sending the plurality of updated node operating parameters to the one or more upstream nodes.

In a further embodiment again, the forward sweep process further comprises: receiving operating parameters from the one or more upstream nodes; calculating a plurality of updated node operating parameters using the operating parameters from the one or more upstream nodes; and sending the plurality of updated node operating, parameters to the one or more downstream nodes.

In another embodiment again, calculating, gradient parameters further comprises exact calculation of the gradient parameters.

In a further additional embodiment, calculating gradient parameters further comprises an approximation of the gradient parameters.

In another additional embodiment, calculating gradient parameters is a distributed process.

In a still yet further embodiment, calculating a gradient step size is evaluated using a line search.

In still yet another embodiment, the node is part of a radial network topology.

In a still further embodiment again, node operating parameters include power injection, current, and impedance.

In still another embodiment again, the node is configured to control operating parameters as components in a single phase power distribution network.

Another further embodiment of the method of the invention includes: the node is configured to control operating parameters as components in a multiphase balanced power distribution network.

Still another further embodiment of the method of the invention includes: the node is configured to control operating parameters as components in a multiphase unbalanced power distribution network.

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 illustrating a coordinator controller in accordance with an embodiment of the invention.

FIG. 4 is a block diagram illustrating a node controller in accordance with an embodiment of the invention.

FIG. 5 is a diagram illustrating a radial network in accordance with an embodiment of the invention.

FIG. 6 is a diagram illustrating the relationship between nodes and operating parameters in accordance with an embodiment of the invention.

FIG. 7 is a flow chart illustrating a process for solving for optimal power flow in a radial network in accordance with an embodiment of the invention.

FIG. 8 is a flow chart illustrating a process to solve for optimal power flow utilizing a gradient projection solution in accordance with an embodiment of the invention.

FIG. 9 is pseudocode illustrating a gradient calculation process utilized in a gradient projection process for solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 10 is pseudocode illustrating a line search process utilized in a gradient projection process solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 11 is pseudocode illustrating a process to solve for optimal power flow a gradient descent process solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 12 is a diagram illustrating backwards sweep process and forwards sweep process utilized in a gradient decent process for solving for optimal power flow.

FIG. 13 is a flow chart illustrating a distributed gradient process utilized in solving for optimal power flow in accordance with an embodiment of the invention.

FIGS. 14A-14B are flow charts illustrating a distributed backwards sweep process and a distributed forwards sweep respectively utilized in a gradient projection process to solve for optimal power flow, in accordance with an embodiment of the invention.

FIG. 15 is pseudocode illustrating a distributed gradient process utilized in solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 16 is a flow chart illustrating a node controller implementing a distributed gradient descent process utilized in solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 17 is a flow chart illustrating a coordinator controller implementing a distributed gradient process utilized in solving for optimal power flow in accordance with an embodiment of the invention.

FIG. 18 is pseudocode illustrating a distributed gradient decent process utilized in solving for optimal power flow in accordance with an embodiment of the invention.

FIGS. 19A-19B are diagrams conceptually illustrating a 47-node (47-bus) network and a 56-node (56-bus) network respectively.

FIG. 20 is a table illustrating parameters of a 47-node (47-bus) network.

FIG. 21 is a table illustrating parameters of a 56-node (56-bus) network.

FIG. 22 is a table illustrating simulated results of a distributed gradient projection process.

DETAILED DESCRIPTION

Turning now to the drawings, systems and methods for distributed control of power distribution systems configured as 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 children. Ancestor nodes and any ancestors of ancestor nodes are typically upstream nodes. Similarly, children nodes and any children of children nodes are typically downstream nodes. Radial networks are commonly utilized in modeling the distribution side of the power grid. In many embodiments, processing nodes are distributed throughout the pourer 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 upstream and/or downstream nodes within the radial network. Hybrid processes use a combination of centralized and distributed processing 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 gradient calculation. In many embodiments, this gradient calculation can be part of a gradient projection process (sometimes referred to as a gradient descent process).

Gradient projection processes typically involve a gradient calculation and choosing a step size to update optimization parameters. Gradient calculations can involve exact and/or approximate calculations of gradients. In various embodiments, gradient calculations can be distributed processes. In many embodiments, a line search can be utilized to update step size. In many embodiments a line search may also be a distributed calculation.

Systems and methods for performing distributed control of radial power distribution networks to achieve OPF and solutions to the distributed 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 equally distributed in equally spaced regular intervals (typically 120 degrees for a commonly used three phase network), the power distribution network is described as a multiphase balanced network. Single phase problems and multiphase balanced network problems can often be solved with similar analysis. In the discussion to follow, networks that distribute power in a single phase or a multiphase balanced manner will both be referred to as single phase networks. An alternative method of distributing power is multiphase unbalanced. In this manner of power distribution, several voltage sources are unevenly spaced. Multiphase unbalanced networks often require different analysis than single phase networks and will be referred to as multiphase networks in the discussion to follow.

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, 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. Connected nodes can be nodes within the power distribution network that are connected by distribution and/or transmission lines and can be controlled by a node controller. 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 n 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 computer systems 208 to monitor calculations made by or to send instructions to multiple node controllers 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 various system configurations are described above with reference to FIG. 2, any number of systems can be utilized to achieve centralized control of nodes within a power distribution network as appropriate to the requirements of specific applications in accordance with embodiments of the invention. Coordinator controllers in accordance with various embodiments of the invention are discussed further below.

Coordinator Controllers

A coordinator controller in accordance with an embodiment of the invention is shown in FIG. 3. In many embodiments, coordinator controller 300 can perform calculations using one or more computer systems to determine coordinator parameters for nodes on a radial network. In the illustrated embodiment, the coordinator controller includes at least one processor 302, an I/O interface 304, and a 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 a coordinator application 308, as well as coordinator parameters 301, node parameters 312, and updated coordinator parameters 314. The coordinator controller can calculate coordinator parameters for the power distribution network using node parameters it receives through the I/O interface from nodes in the network and/or from node controllers similar to node controllers described below with respect to FIG. 4. The coordinator application 308 will be discussed in greater detail below and can enable the coordinator controller to broadcast coordinator parameters to nodes to enable the power distribution network to solve for optimal power flow using a distributed gradient projection process. In several embodiments, coordinator parameters can enable nodes to act in a concurrent manner while node controllers solve for optimal power flow. It should be readily apparent that the coordinator controller can be at a node (for example at the root node) and/or at a centralized computer system in the power distribution network. Although a number of different coordinator controllers are described above with reference to FIG. 3, any of a variety of computing systems can be utilized to coordinate control of nodes in a power distribution system as appropriate to the requirements of specific applications in accordance with many embodiments of the invention. Node controllers in accordance with several embodiments of the invention are discussed further below.

Node Controllers

A node controller in accordance with an embodiment of the invention is shown in FIG. 4. In various embodiments, node controller 400 can control the operating parameters of one or more nodes in a (radial) power distribution network. In the illustrated embodiment, the node controller includes at least one processor 402, an I/O interface 404, and memory 406. The at least one processor 402, 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 406 includes software including a node controller application 408 as well as node parameters 410, neighboring node parameters 412, coordinator parameters 414, and updated parameters 416. A node can calculate updated parameters by using a combination of its own node parameters, and/or neighboring node parameters received through the I/O interface from upstream and/or downstream nodes, and/or coordinator parameters received from a coordinator controller similar to a coordinator controller described above with reference to FIG. 2. The distributed optimal power flow application 408 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. These distributed calculations performed on the node parameters can specifically involve a gradient projection process to solve for optimal power flow in a distributed manner. 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) power injection, current, resistance, impedance, and/or voltage. Although a number of different node controller implementations are described above with respect to FIG. 4, 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. Radial networks in accordance with embodiments of the invention are discussed further below.

Radial Network Models

A radial network in accordance with an embodiment of the invention is shown in FIG. 5. In various embodiments, radial network 500 includes a node 502. Node 502 has an ancestor node 504 and one or more children nodes 506. Ancestor nodes are upstream from node 502. Children nodes are downstream from node 502. A radial network also has a unique root node 508. A detailed discussion of these nodes is provided further below. Although many radial networks are described above with reference to FIGS. 3 and 4, 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 relationships between nodes and operation parameters in an OPF model are discussed further below.

The relationship between nodes and operating parameters in an OPF model in accordance with an embodiment of the invention is shown in FIG. 6. A node 602 has a unique ancestor node 604. Node 602 and unique ancestor node 604 are connected by line 614 and have a unique root node 606. Both node 602 and unique ancestor node 604 have a series of operating parameters. In many embodiments of the invention, operating parameters for node 602 include voltage 608 and power injection 616. Unique ancestor node 608 has corresponding voltage 610 and unique root node has a corresponding voltage 612. The line 614 also has operating parameters which for example include an impedance value as well as a current and/or power injection 618. The relationship between nodes and operating values will be discussed in greater detail below. Although various node operating parameters are described above with respect to FIG. 6, any of a variety of operating parameters can be controlled to achieve OPF as appropriate to the requirements of specific applications in accordance with embodiments of the invention. The OPF problem is discussed further below.

OPF Problems

To highlight the main components of the distributed process without worrying about the complication of multiple phases (which will be discussed later below), the network will be assumed to be single-phase and radial. In particular, the node (bus) voltages V_i, node (bus) power injections s_i, and branch power flows S_ij are scalars.

A distribution network can be modeled as the following: For each node (bus) i∈N, let ν_i=|V_i|² denote the square of its complex voltage magnitude, e.g., if the voltage is V_i=1.05∠120° per unit, then ν_i=1.05². Note that ν₀ is fixed and given, p_i=Re(s_i) and q_i=Im(s_i) can denote the real and reactive power injections respectively. _i can denote the unique path from node (bus) 0 to node (bus) i. Since the network is radial, the path _i is well-defined. For each line (i,j)∈ε, let l_ij=|I_ij|² denote the square of the complex current magnitude, e.g., if the current is I_ij=0.5∠10°, then l_ij=0.5². Let P_ij=Re(S_ij) and Q_ij=Im(S_ij) denote the real and reactive power flows respectively. Let r_ij=Re(z_ij) and x_ij=Im(z_ij) denote the resistance and impedance respectively.

In many embodiments, the following simplified statement of the OFF problem can be utilized:

$\begin{array}{cc}\min \sum _{i=0}^n{a}_i{p}_i^2+b_i{p}_i\mathrm{over}{p}_i\mathrm{and}{q}_i\mathrm{for}i\in p_0,q_0,v_i\mathrm{for}i\in {+}^{},P_{\mathrm{ij}},Q_{\mathrm{ij}}\mathrm{and}{\ell }_{\mathrm{ij}}\mathrm{for}i\to j;& \left(1a\right)\\ \sum _{i:i\to j}\left(P_{\mathrm{ij}}-r_{\mathrm{ij}}{\ell }_{\mathrm{ij}}\right)+p_j=\sum _{k:j\to k}{P}_{\mathrm{jk}},j\in ;& \left(1b\right)\\ \sum _{i:i\to j}\left(Q_{\mathrm{ij}}-x_{\mathrm{ij}}{\ell }_{\mathrm{ij}}\right)+q_j=\sum _{k:j\to k}=Q_{\mathrm{jk}},j\in {+}^{};& \left(1c\right)\\ v_i-v_j=2\left(r_{\mathrm{ij}}{P}_{\mathrm{ij}}+x_{\mathrm{ij}}{Q}_{\mathrm{ij}}\right)-{z_{\mathrm{ij}}}^2{\ell }_{\mathrm{ij}},i\to j;& \left(1d\right)\\ {\ell }_{\mathrm{ij}}=\frac{P_{\mathrm{ij}}^2+Q_{\mathrm{ij}}^2}{v_i},i\to j;& \left(1e\right)\\ {\underset{\_}{v}}_i\le v_i\le {\stackrel{\_}{v}}_i,i\in {+}^{};& \left(1f\right)\\ p_i\in \left[{\underset{\_}{p}}_i,{\stackrel{\_}{p}}_i\right],q_i\in \left[{\underset{\_}{q}}_i,{\stackrel{\_}{q}}_i\right],i\in {+}^{}.& \left(1g\right)\end{array}$

The objective function (1a) is assumed to be separable, quadratic, and purely a function of p. Equations (1b)(1e) describe the physical laws that govern the power flow. Equation (1f) is the voltage regulation constraints, and equation (1g) is the power injection constraints.

Basic and Derived Variables

Given the substation voltage ν₀ and the branch node (bus) power injections s_i for i∈⁺, there exists a unique practical power flow solution (, [l_ij]_i→j, [S_ij]_i→j, s₀). Hence, the optimization variables can be classified into two categories: (1) Basic variables that are controllable. These variables include the branch node (bus) real and reactive power injections p_i and q_i for i∈³⁰. (2) Derived variables that are uniquely determined by power flow equations (1b)-(1e) after specifying the basic variables. These variables include the substation power injection p₀, q₀, the branch node (bus) voltages ν_i for i∈⁺, and line flows P_ij, Q_ij, l_ij for i→j.

In some embodiments, x=(p_l, . . . , p_n, q_l, . . . q_n) can denote all basic variables. Then derived variables can be functions of x, i.e.,p₀=p₀(x), q₀=q₀(x);ν_i=ν_i(x), i∈⁺;P_ij=P_ij(x), Q_ij(x), l_ij=l_ij(x), i→j. Removing the derived variables can transform the statement of the OFF problem (1) into the following form:

$\begin{array}{cc}\min a_0{p}_0^2\left(x\right)+b_0{p}_0\left(x\right)+\sum _{i=1}^n\left(a_i{p}_i^2+b_i{p}_i\right)\mathrm{over}xs.t.{\underset{\_}{v}}_i\le v_i\left(x\right)\le {\stackrel{\_}{v}}_i,i\in {𝒩}^{+};& \left(2a\right)\\ {\underset{\_}{p}}_i\le p_i\le {\stackrel{\_}{p}}_i,{\underset{\_}{q}}_i\le q_i\le {\stackrel{\_}{q}}_i,i\in {𝒩}^{+}.& \left(2b\right)\end{array}$

While (2) is equivalent to (1) for radial networks, (2) has much fewer optimization variables than (1) and in various embodiments can be more efficient to compute. Furthermore, the derived variables (p₀, +) will be automatically computed by power flow physics once the basic variable x is computed. This motivates an iterative procedure for solving OPF: in each iteration, first update the basic variable x and then let the derived variables be automatically computed by computations based upon power flow physics.

Gradient Projection Processes

A process for achieving optimal power flow is illustrated in FIG. 7. The process 700 includes receiving 702 optimal power flow parameters at one or more nodes and/or one or more coordinators. Optimal power how can then be calculated 704 using a gradient projection process. The calculations for gradient projection will be discussed in detail further below. In various embodiments, gradient projection includes calculating derivatives and/or choosing a step size for updating optimization parameters. Although a number of processes for achieving optimal power flow are described with reference to FIG. 7 any of a variety of gradient projection processes can be utilized to solve for optimal power flow as appropriate to the requirements of specific applications in accordance with embodiments of the invention. Gradient decent processes in accordance with embodiments of the invention are described further below.

A process to solve for optimal power flow utilizing a gradient projection solution is illustrated in FIG. 8. The process 800 includes receiving 802 optimal power flow parameters. Parameters including, an iterative parameter are used to solve 804 for optimal power flow. A backwards forwards sweep can be calculated 806 to obtain power flow parameters. A backwards forwards sweep will be discussed in detail further below. Gradients are calculated 808 using, power flow parameters. Gradient calculations can be exact and/or approximations, and various ways to calculate gradients including in distributed processes will be described further below. Update parameters can be calculated 810 using a line search 810. Parameters including the iterative parameter are updated 812 using gradients and update parameters. In various embodiments, some or all of these updated parameter values can come from a coordinator. End conditions are checked 814. If an end condition 814 has not been reached, the parameters including the iterative parameter are used to once again solve 804 for optimal power flow. If an end condition 814 has been reached, optimal power flow is returned 816 and the process completes. Although a variety of processes for gradient descent are illustrated with respect to FIG. 8, any of a variety of processes to iteratively calculate optimal power flow in a power distribution network can be utilized as appropriate to the requirements of specific applications in accordance with embodiments of the invention.

The process proposed to solve the OPF problem (2) is a gradient projection process. In each iteration of the process, derivatives of the modified objective function with respect to the basic variables are estimated, and then used as the negative direction of updating the basic variables. Line search is implemented to determine the step sizes of basic variable updates so as to ensure the convergence of the process.

Exact Derivatives

One of the reasons that the gradient projection process proposed can be adept is that the derivatives can be estimated efficiently. First a process for exactly computing derivatives will be illustrated. A process for estimating derivatives will be described further below.

As noted above, gradients computed during OPF calculations can be computed exactly. Since the network is radial, (1b)-(1e) imply

${\partial }_x{P}_{\mathrm{ij}}=r_{\mathrm{ij}}{\partial }_x{\ell }_{\mathrm{ij}}-{\partial }_x{p}_j+\sum _{k:j\to k}{\partial }_x{P}_{\mathrm{jk}},{\partial }_x{Q}_{\mathrm{ij}}=x_{\mathrm{ij}}{\partial }_x{\ell }_{\mathrm{ij}}-{\partial }_x{q}_j+\sum _{k:j\to k}{\partial }_x{Q}_{\mathrm{jk}},{\partial }_x{v}_j={\partial }_x{v}_i-2\left(r_{\mathrm{ij}}{\partial }_x{P}_{\mathrm{ij}}+x_{\mathrm{ij}}{\partial }_x{Q}_{\mathrm{ij}}\right)+{z_{\mathrm{ij}}}^2{\partial }_x{\ell }_{\mathrm{ij}},{\partial }_x{\ell }_{\mathrm{ij}}=\frac{2P_{\mathrm{ij}}}{v_i}{\partial }_x{P}_{\mathrm{ij}}+\frac{2Q_{\mathrm{ij}}}{v_i}{\partial }_x{Q}_{\mathrm{ij}}-\frac{{\ell }_{\mathrm{ij}}}{v_i}{\partial }_x{v}_i$ for i→j. I can denote the 2×2 identity matrix and ∂_xl_ij can be removed to obtain

$\left[I-\frac{2}{v_i}\left(\begin{array}{c}{r}_{\mathrm{ij}}\\ x_{\mathrm{ij}}\end{array}\right)\left(\begin{array}{cc}{P}_{\mathrm{ij}}& Q_{\mathrm{ij}}\end{array}\right)\right]\left(\begin{array}{c}{\partial }_x{P}_{\mathrm{ij}}\\ {\partial }_x{Q}_{\mathrm{ij}}\end{array}\right)=\sum _{k:j\to k}\left(\begin{array}{c}{\partial }_x{P}_{\mathrm{jk}}\\ {\partial }_x{Q}_{\mathrm{jk}}\end{array}\right)-\left(\begin{array}{c}{\partial }_x{p}_j\\ {\partial }_x{q}_j\end{array}\right)-\left(\begin{array}{c}{r}_{\mathrm{ij}}\\ x_{\mathrm{ij}}\end{array}\right)\frac{{\ell }_{\mathrm{ij}}}{v_i}{\partial }_x{v}_i,{\partial }_x{v}_j=\left(1-{z_{\mathrm{ij}}}^2\frac{{\ell }_{\mathrm{ij}}}{v_i}\right){\partial }_x{v}_i-2\left(r_{\mathrm{ij}}-{z_{\mathrm{ij}}}^2\frac{P_{\mathrm{ij}}}{v_i}\right){\partial }_x{P}_{\mathrm{ij}}-2\left(x_{\mathrm{ij}}-{z_{\mathrm{ij}}}^2\frac{Q_{\mathrm{ij}}}{v_i}\right){\partial }_x{Q}_{\mathrm{ij}}$ for i→j. Pseudocode that can be utilized to solve for optimal power flow utilizing a gradient calculation process to perform an exact derivative calculation in accordance with many embodiments of the invention is illustrated in FIG. 9. Hence, the gradients ∂_x(P, Q, ν, p₀, q₀) can be computed by a process similar to process 1 as illustrated in FIG. 9.Approximated Derivatives

To avoid the iterative procedure used in process 1 described above and improve computational efficiency, the gradients can be estimated as follows. Note that since the current terms in (1b)-(1d) are much smaller than the other terms in practice, (P, Q, ν) can be estimated by ({circumflex over (P)}, {circumflex over (Q)}, {circumflex over (ν)}) defined as

$\sum _{i:i\to j}{\hat{P}}_{\mathrm{ij}}+p_j=\sum _{k:j\to k}{\hat{P}}_{\mathrm{jk}},j\in 𝒩;$ $\sum _{i:i\to j}{\hat{Q}}_{\mathrm{ij}}+q_j=\sum _{k:j\to k}{\hat{Q}}_{\mathrm{jk}},j\in 𝒩;$ ${\hat{v}}_i-{\hat{v}}_j=2\left(r_{\mathrm{ij}}{\hat{P}}_{\mathrm{ij}}+x_{\mathrm{ij}}{\hat{Q}}_{\mathrm{ij}}\right),i\to j.$

Note that the equations are linear and in many embodiments straightforward to obtain

${\partial }_x{\hat{P}}_{\mathrm{ij}}=\sum _{k:j\to k}{\partial }_x{\hat{P}}_{\mathrm{jk}}-{\partial }_x{p}_j,i\to j;$ ${\partial }_x{\hat{Q}}_{\mathrm{ij}}=\sum _{k:j\to k}{\partial }_x{\hat{Q}}_{\mathrm{jk}}-{\partial }_x{q}_j,i\to j;$ ${\partial }_x{\hat{v}}_j={\partial }_x{\hat{v}}_i-2r_{\mathrm{ij}}{\partial }_x{\hat{P}}_{\mathrm{ij}}-2x_{\mathrm{ij}}{\partial }_x{\hat{Q}}_{\mathrm{ij}},i\to j.$

In several embodiments ij can denote the joint of node (bus) i and j for i, j∈, and R_i can denote the total resistance from node (bus) 0 to node (bus) i for i∈. Then ∂_x({circumflex over (P)}, {circumflex over (Q)}, {circumflex over (ν)}) has the following closed-form expression:∂_pk{circumflex over (P)}_ij=−, ∂_qk{circumflex over (P)}_ij=0, k=1,2, . . . n, i→j;  (3a)∂_pk{circumflex over (Q)}_ij=0, ∂_qk{circumflex over (Q)}_ij=−, k=1,2, . . . n,i→j;  (3b)∂_pkν_i=2, ∂_qkν_i=2, k=1,2, . . . n, i→j;  (3c)∂_x(P, Q, ν) can be approximated by ∂_x({circumflex over (P)}, {circumflex over (Q)}, {circumflex over (ν)}), i.e.,∂_x(P,Q,ν)≈∂_x({circumflex over (P)},{circumflex over (Q)},{circumflex over (ν)}).Note that ∂_x({circumflex over (P)}, {circumflex over (Q)}, {circumflex over (ν)}) is a constant that does not depend on (P, Q, ν), and therefore can be computed once ahead of time.

Finally, ∂_x(p₀, q₀) can be approximated as follows, (1b) can be summed up for j∈ to obtain

$\sum _{i=0}^n{p}_i=\sum _{i\to j}{r}_{\mathrm{ij}}{\ell }_{\mathrm{ij}}=\sum _{i\to j}{r}_{\mathrm{ij}}\frac{P_{\mathrm{ij}}^2+Q_{\mathrm{ij}}^2}{v_i}.$ Hence,

$\begin{array}{c}{\partial }_{p_i}{p}_0=-1+\sum _{k\to l}{r}_{\mathrm{kl}}{\partial }_{p_i}\left(\frac{P_{\mathrm{kl}}^2+Q_{\mathrm{kl}}^2}{v_k}\right)\\ =-1+\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{P}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{Q}_{\mathrm{kl}}-\frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{v}_k\right)\\ \approx -1+\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{P}}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{Q}}_{\mathrm{kl}}-\frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{v}}_k\right)\end{array}$ $\begin{array}{c}{\partial }_{q_i}{p}_0=\sum _{k\to l}{r}_{\mathrm{kl}}{\partial }_{q_i}\left(\frac{P_{\mathrm{kl}}^2+Q_{\mathrm{kl}}^2}{v_k}\right)\\ =\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{P}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{Q}_{\mathrm{kl}}-\frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{v}_k\right)\\ \approx \sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{\hat{P}}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{\hat{Q}}_{\mathrm{kl}}-\frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{q_i}{\hat{v}}_k\right)\end{array}$ $\mathrm{for}i=1,2,\dots ,n.$ Modified Objective Functions

To enable a distributed process, i.e., each node (bus) i updates its own (p_i, q_i) locally, the constraints of OFF can be decoupled, i.e., constraints of the form (2b) are easy to handle while constraints of the form (2a) should be avoided.

To avoid coupled constraints (2a), a log-barrier function can be added to the objective as

$L\left(x;\mu \right)=a_0{p}_0^2\left(x\right)+b_0{p}_0\left(x\right)+\sum _{i=1}^n\left(a_i{p}_i^2+b_i{p}_i\right)-\underset{\_}{\mu }\sum _{i=1}^n\ln \left({\stackrel{\_}{v}}_i-{\underset{\_}{v}}_i\right)-\stackrel{\_}{\mu }\sum _{i=1}^n\ln \left({\stackrel{\_}{v}}_i-v_i\left(x\right)\right)$

where μ=(μ, μ)>0 component wise. Note that ν_i≤ν_i(x)≤ν_i for i∈⁺ is enforced since

$\underset{t\downarrow {\underset{\_}{v}}_i}{\lim }-\underset{\_}{\mu }\ln \left(t-{\underset{\_}{v}}_i\right)=\infty ,\underset{t\uparrow {\stackrel{\_}{v}}_i}{\lim }-\stackrel{\_}{\mu }\ln \left({\stackrel{\_}{v}}_i-t\right)=\infty ,i\in {𝒩}^{+}$ and OPF seeks to minimize L. Besides,

$\underset{\mu \downarrow 0}{\lim }L\left(x;\mu \right)=a_0{p}_0^2\left(x\right)+b_0{p}_0\left(x\right)+\sum _{i=1}^n\left(a_i{p}_i^2+b_i{p}_i\right)$ and therefore solving OPF is similar to minimize L(x;μ) with small enough μ.

To summarize, in various embodiments the distributed process can solveOPF(μ):min L(x;μ)over p₁, . . . , p_n,q₁, . . . , q_n;s.t. p_i≤p_i≤p_i,q_i≤q_i≤q_i, i=1, 2, . . . ,n for a decreasing sequence of μ. OPF(μ) will be solved for a specific μ below.Gradient Projection

In many embodiments, there are two key steps in a gradient projection process: 1) compute or approximate) the gradients; and 2) choose a step size to update the optimization variables.

The gradients can be approximated using (3). In particular,

$\begin{array}{cc}\begin{array}{c}{\partial }_{p_i}L=\left(2a_0{p}_0+b_0\right){\partial }_{p_i}{p}_0+\left(2a_i{p}_i+b_i\right)-\sum _{k=1}^n\left[\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right]\\ {\partial }_{p_i}{v}_k\\ =\left(2a_0{p}_0+b_0\right)[-1+\sum _{k\to l}{r}_{\mathrm{kl}}(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{P}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{Q}_{\mathrm{kl}}-\\ \frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{v}_k)]+\\ \left(2a_i{p}_i+b_i\right)-\sum _{k=1}^n\left[\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right]{\partial }_{p_i}{v}_k\\ \approx \left(2a_0{p}_0+b_0\right)[-1+\sum _{k\to l}{r}_{\mathrm{kl}}(\frac{2P_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{P}}_{\mathrm{kl}}+\frac{2Q_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{Q}}_{\mathrm{kl}}-\\ \frac{{\ell }_{\mathrm{kl}}}{v_k}{\partial }_{p_i}{\hat{v}}_k)]+\left(2a_i{p}_i+b_i\right)-\sum _{k=1}^n\left[\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right]{\partial }_{p_i}{\hat{v}}_k\\ =-\left(2a_0{p}_0+b_0\right)\left[1+\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{𝕝}_{l\in {𝒫}_i}+\frac{2{\ell }_{\mathrm{kl}}}{v_k}{R}_{i\bigwedge k}\right)\right]+\\ \left(2a_i{p}_i+b_i\right)-\sum _{k=1}^n\left[\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right]2R_{i\bigwedge k}\end{array}\mathrm{for}i=1,2,\dots ,n.\mathrm{Similarly},& \left(4\right)\\ {\partial }_{q_i}L\approx -\left(2a_0{p}_0+b_0\right)\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2Q_{\mathrm{kl}}}{v_k}{𝕝}_{l\in {𝒫}_i}+\frac{2{\ell }_{\mathrm{kl}}}{v_k}{X}_{i\bigwedge k}\right)-\sum _{k=1}^n\left[\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right]2X_{i\bigwedge k}\mathrm{for}i=1,2,\dots ,n.& \left(5\right)\end{array}$

The step size can be determined by doing a line search along the direction of −∂_(p,q)L, i.e., back off the step size until the modified Objective function can be well-approximated by its linearization around the current solution point. Three parameters α (determine the back off speed, set to 0.5 in the current implementation), β (criteria for the linearization of the objective to be accurate enough, set to 0.5 in the current implementation), and ε (criteria for the progress to be too slow, set to 1e-4 in the current implementation) are needed in the line search step. Pseudocode that can be utilized to implement line search processes utilized for solving for optimal power flow by node controllers to achieve gradient descent in accordance with many embodiments of the invention is illustrated in process 2 in FIG. 10.

The introduction of ε in the “if” branch in Step 7) of process 2 illustrated in FIG. 10 is to stop the iterations when progress gets too slow, i.e., ∥Δp∥+∥Δq∥≤ε. When this happens, stopFlag is set to 1 and iterations are stopped. Otherwise, a large number of iterations will run without updating (p,q) significantly.

With the “if” branch in Step 7) of process 2 illustrated in FIG. 10, it is possible that L(p′,q′)>L(p,q). In this case, (p′,q′) is set to (p,q) to ensure that new point (p′,q′) does not have a larger objective value than (p,q).

In many embodiments, a point (p,q) is a local optimum for minimizing L if∂_piL,{tilde over (p)}_i−p_i≥0, ∀{tilde over (p)}_i∈(p_i,p_i), ∀i∈⁺;∂_qiL,{tilde over (q)}_i−q_i≥0, ∀{tilde over (q)}_i∈(q_i, q_i), ∀i∈⁺.

It can be proven that if (p,q)=(p′,q′), then (p,q) is a local optimum. Additionally, it can be proven that if ε=0, then the input (p,q) and output (p′,q′) of process 2 satisfyL(p′,q′)≤L(p,q).if and only if (p′,q′)=(p,q). These imply that L(p′,q′)<L(p,q) unless (p,q) is a local optimum, in which case (p′,q′)=(p,q).Distributed Gradient Projection Processes

Pseudocode that can be utilized to implement distributed gradient projection processes utilized to solve for optimal power flow in accordance with various embodiments of the invention is illustrated in FIG. 11. The distributed gradient projection process is summarized in process 3 illustrated in FIG. 11, which solves OPF(μ) with different values of μ. In particular, μ₁, μ₂, . . . , μ_K can denote a sequence of μ that approaches 0 and x=(p,q). Given a feasible initial point x⁽⁰⁾, process 3 solves OPF(μ₁) with initial point x⁽⁰⁾ to obtain x⁽¹⁾, then solves OPF(μ₂) with initial point x⁽¹⁾ to obtain x⁽²⁾, . . . , and finally solves OPF(μ_K) with initial point x^(K-1) to obtain the final output solution x^(K).

To solve each OPF(μ), process 3 repeatedly calculates the gradient ∂_xL(x;μ) and does a line search along the direction of −∂_xL to update x, until stopFlag=1, which indicates numerically no further improvements can be made. It can be proven that a local minimum of L(x;μ) is obtained in each inner loop of process 3.

Distributed Implementation

An important advantage of process 3 is that it can be implemented in a distributed way. The infrastructure required to implement process 3 in a distributed way is described as follows. (1) There is a node controller (agent) at each branch node (bus) i∈⁺, and that there is a centralized coordinator at the substation node (bus) 0. (2) Each node controller (agent) i has access to its local power injection (p_i,q_i), local power flow (P_ij,Q_ij), and local voltage ν_i. It can directly communicate with its neighbors, i.e., node controller (agent) i can communicate with node controller (agent) j if and only if i˜j. (3) The centralized coordinator has access to the substation power injection (p₀,q₀) and can communicate with all node controllers (agents) in the network.

In various embodiments, to implement process 3 in a distributed way, only nodes (buses) i whose power injections (p_i, q_i) can be controlled need to have node controllers (agents). For example, if p_i=p_i and q_i=q_i, then node (bus) i does not need to have a node controller (agent).

There are two key components in the distributed implementation of process 3: 1) compute gradient ∂_(p,q)L(p*,q*;μ) in a distributed way; and 2) run line search (process 2 in a distributed way.

Distributed Gradient Computation

The approximated gradients (see equations (4)-(5)) can be computed in a backward forward sweep fashion as described below.

In several embodiments. Down(i):={j∈|i∈_i} can denote the downstream nodes (buses) of node (bus) i∈. Define

$c_i=\sum _{k=1}^{n}2R_{i\bigwedge k}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right),i\in 𝒩;$ $d_i=\sum _{k=1}^{n}2X_{i\bigwedge k}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right),i\in 𝒩;$ $e_i=\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2P_{\mathrm{kl}}}{v_k}{𝕝}_{l\in {𝒫}_i}+\frac{{\ell }_{\mathrm{kl}}}{v_k}2R_{i\bigwedge k}\right),i\in 𝒩;$ $f_i=\sum _{k\to l}{r}_{\mathrm{kl}}\left(\frac{2Q_{\mathrm{kl}}}{v_k}{𝕝}_{l\in {𝒫}_i}+\frac{{\ell }_{\mathrm{kl}}}{v_k}2X_{i\bigwedge k}\right),i\in 𝒩.$ Then the approximated gradients (4)-(5) can be simplified as∂_piL=−(2a₀p₀+b₀)(1+e_i)+(2a_ip_i+b_i)−c_i, i∈⁺;  (6a)∂_qiL=−(2a₀p₀+b₀)f_i−d_i, i∈⁺.  (6b)The centralized coordinator has access to p₀ can broadcast 2α₀p₀+b₀ to all branch nodes (buses). Each node controller (agent) i knows p_i and can easily compute 2α_ip_i+b_i. Hence, the main challenge is to compute c_i, d_i, e_i, f_i in a distributed manner.

The quantities c_i, d_i, e_i, f_i can be computed recursively. To derive the recursive equations, note that for each i→j, one hasR_i∧k−R_j∧k=−r_{ijk∈Down(j)} and therefore

$\begin{array}{c}{c}_i-c_j=\sum _{k=1}^{n}2\left(R_{i\bigwedge k}-R_{j\bigwedge k}\right)\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right)\\ =-\sum _{k=1}^{n}2r_{\mathrm{ij}}{𝕝}_{k\in \mathrm{Down}\left(j\right)}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right)\\ =-2r_{\mathrm{ij}}\sum _{k\in \mathrm{Down}\left(j\right)}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right).\end{array}$

Similarly,

$d_i-d_j=-2x_{\mathrm{ij}}\sum _{k\in \mathrm{Down}\left(j\right)}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right).$

Besides,−=and therefore

$\begin{array}{c}{e}_i-e_j=\sum _{k\to l}{r}_{\mathrm{kl}}\left[\frac{2P_{\mathrm{kl}}}{v_k}\left({𝕝}_{l\in {𝒫}_i}-{𝕝}_{l\in {𝒫}_j}\right)+\frac{{\ell }_{\mathrm{kl}}}{v_k}2\left(R_{i\bigwedge k}-R_{j\bigwedge k}\right)\right]\\ =-\sum _{k\to l}{r}_{\mathrm{kl}}\left[\frac{2P_{\mathrm{kl}}}{v_k}{𝕝}_{l=j}+\frac{{\ell }_{\mathrm{kl}}}{v_k}2r_{\mathrm{ij}}{𝕝}_{k\in \mathrm{Down}\left(j\right)}\right]\\ =-\frac{2r_{\mathrm{ij}}{P}_{\mathrm{ij}}}{v_i}-r_{\mathrm{ij}}\sum _{k\in \mathrm{Down}\left(j\right)}\frac{2r_{\mathrm{kl}}{\ell }_{\mathrm{kl}}}{v_k}.\end{array}$

Similarly,

$f_i-f_j=-\frac{2r_{\mathrm{ij}}{Q}_{\mathrm{ij}}}{v_i}-x_{\mathrm{ij}}\sum _{k\in \mathrm{Down}\left(j\right)}\frac{2r_{\mathrm{kl}}{\ell }_{\mathrm{kl}}}{v_k}.$ Hence, g_i and h_i can be defined as:

$g_i=\sum _{k\in \mathrm{Down}\left(i\right)}\left(\frac{\underset{\_}{\mu }}{v_k-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{v_k-{\stackrel{\_}{v}}_k}\right),i\in {𝒩}^{+};$ $h_i=\sum _{k\in \mathrm{Down}\left(i\right)}\frac{2r_{\mathrm{kl}}{\ell }_{\mathrm{kl}}}{v_k},i\in {𝒩}^{+},$ then c_i, d_i, e_i, f_i can be computed recursively as

$\begin{array}{cc}{c}_j=c_i+2r_{\mathrm{ij}}{g}_j,i\to j;& \left(7a\right)\\ d_j=d_i+2x_{\mathrm{ij}}{g}_j,i\to j;& \left(7b\right)\\ e_j=e_i+\frac{2r_{\mathrm{ij}}{P}_{\mathrm{ij}}}{v_i}+r_{\mathrm{ij}}{h}_j,i\to j;& \left(7c\right)\\ f_j=f_i+\frac{2r_{\mathrm{ij}}{Q}_{\mathrm{ij}}}{v_i}+x_{\mathrm{ij}}{h}_j,i\to j.& \left(7d\right)\end{array}$ Distributed Gradient Projection Processes

A distributed gradient process utilized in solving for optimal power flow to implement distributed gradient projection is illustrated in FIG. 13. The process 1300 includes performing 1302 a backwards sweep to calculate parameters utilized in OPF calculations. Additionally, a forwards sweep is performed to calculate parameters. Various distributed backwards and forwards sweep processes are described in detail below with reference to FIGS. 14A and 14B. Gradients are calculated 1306. Gradients can be calculated with exact derivative calculations and/or with approximations as described above. Pseudocode to implement a distributed gradient projection process similar to process 1300 is described further below in FIG. 15. Although a variety of gradient projection processes are described above with respect to FIG. 13, any of a variety of distributed gradient projection processes to solve for optimal power flow in a power distribution network can be utilized as appropriate to the requirements of specific applications in accordance with embodiments of the invention.

A process to implement a distributed backwards sweep utilized in a gradient projection process to solve for optimal power flow is illustrated in FIG. 14A. The process 1400 includes a node receiving 1402 parameters from downstream neighbor nodes. The node calculates 1404 updated parameters with parameters received from downstream neighbor nodes. Updated parameters are sent 1406 to upstream neighbor nodes, in some embodiments, parameters can include (but are not limited to) voltage, power injection, current, resistance, and/or impedance. In many embodiments, backwards sweep begins at the leaf nodes and completes at the root node, but it should be readily appreciated that it may be performed over a subsection of nodes in a power distribution network as specific applications require.

A process to implement a distributed forwards sweep utilized in a gradient projection process to solve for optimal power flow is illustrated in FIG. 14B. The process 1450 includes a node receiving 1452 parameters from upstream neighbor nodes. The node calculates 1454 updated parameters using parameters from upstream neighbor nodes. Updated parameters are sent 1456 to downstream neighbor nodes, in some embodiments, parameters can include (but are not limited to) voltage, power injection, current, resistance, and/or impedance. In various embodiments, forwards sweep generally starts after a backwards sweep process is finished. In many embodiments, the forwards sweep begins at the root node and completes at the leaf nodes. However, it should be readily appreciated that it may be performed over a subsection of nodes in a power distribution network as specific applications require.

Although a variety of processes for performing backwards and forwards sweeps have been described above with respect to FIGS. 14A and 14B, any of a variety of distributed processes to propagate information throughout nodes in a power distribution network may be utilized as appropriate to the requirements of specific applications in accordance with embodiments of the invention.

Pseudocode that can be utilized to implement distributed gradient projection processes utilized in solving for optimal power flow executed by node controllers is illustrated in FIG. 15. The illustrated process 4 is similar to those outlined above with respect to FIG. 13. To summarize, the distributed gradient projection process is given in process 4, where:={i∈⁺|∃j such that i→j}denotes the set of leaf nodes (buses). It includes three main steps:

• • S1 Backward sweep to compute g, h: for each node controller (agent) i∈⁺, when all its downstream neighbors j have computed g_j, h_j, computes g_i, h_i, according to

$\begin{array}{cc}{e}_i=\frac{\underset{\_}{\mu }}{v_i-{\underset{\_}{v}}_i}+\frac{\stackrel{\_}{\mu }}{v_i-{\stackrel{\_}{v}}_i}+\sum _{j:i\to j}{e}_j;f_i=\sum _{j:i\to j}\left(\frac{2r_{\mathrm{ij}}{\ell }_{\mathrm{ij}}}{v_i}+f_j\right),& \left(8\right)\end{array}$

• • S2 Forward sweep to compute c, d, e, f: for each node controller (agent) h∈⁺, when all its upstream neighbors i have computed c_i, d_i, e_i+2r_ijP_ij/ν_i, f_i+2r_ijQ_ij/ν_i, computes c_j, d_j, e_j, f_j as in (7).
  • S3 Compute ∂_(p,q)L: each node controller (agent) i∈⁺ computes ∂_piL and ∂_qiL according to (6).Distributed Gradient Descent Processes

A distributed gradient decent process similar to process 3 in FIG. 11 can be implemented in a distributed manner by calling a distributed implementation of the inner loop for a sequence μ₁, μ₂, . . . , μ_K) of decreasing μ. In various embodiments, a node controller and a coordinator controller can act in concert to achieve optimal power flow in a distributed manner utilizing a gradient decent process, which will be described further below with respect to the processes conceptually illustrated in FIG. 16 and FIG. 17. A distributed gradient decent process can also include a distributed line search process similar to line search process 2 illustrated in FIG. 10.

A node controller implementing a distributed gradient descent process utilized in solving for optimal power flow is illustrated in FIG. 16. The process 1600 includes a node receiving 1602 an iterative parameter from the coordinator. In many embodiments, this iterative parameter can be a value from a decreasing sequence that approaches 0. The node calculates 1604 the old power injection parameters and sends them to the coordinator. The gradients are calculated 1608. A step size is received 1610 from the coordinator. New power injection parameters are calculated 1612 using information including gradients and step size, and the new power injection parameters are sent to the coordinator. The node checks if the voltage parameters are within a specified range 1612. If the voltage parameters are not within a specified range 1612, a signal is sent 1616 to the coordinator for a new step size. A new step size is received 1610 and the process repeats. If the voltage parameter is within a specified range 1612, end conditions 1618 are checked. If end conditions 1618 have not been met, the values of old power injection parameters are updated 1620 such that they now equal the values of the new power injection parameters. These parameters are used to calculate 1604 remaining old power injection parameters to be sent to the coordinator and the process repeats. If the end conditions 1618 have been met old or new power injection parameters are returned 1622 based on the end conditions and the process completes. It should readily be apparent that once complete, process 1600 will repeat again for additional iterative parameters until the system has (but is not limited to) solved for a local optimum of the optimal power flow problem, in several embodiments, a node controller similar to node controller 400 illustrated in FIG. 4 can implement this process.

A a distributed gradient descent process utilized in solving for optimal power flow that can be implemented by a coordinator controller is illustrated in FIG. 17. The process 1700 includes broadcasting 1702 an iterative parameter to one or more nodes in the power distribution network. Old power injection parameters are received 1704 from one or more nodes. An objective value parameter is calculated 1706. A step size is initialized 1708 and broadcast 1710 to one or more nodes. New power injection parameters are received 1712 from one or more nodes. If a voltage parameter out of range signal 1714 is received from one or more nodes, the step size is modified 1716, the modified step size is broadcast 1710 and the process repeats. If a voltage parameter out of range signal 1714 is received from one or more nodes, termination parameters are calculated 1718. End conditions 1720 are checked. If end conditions 1720 are not met, a signal is broadcast 1722 to one or more nodes to update power injection parameters. If end conditions 1720 are met, end conditions are broadcast 1720 and the process completes. It should readily be apparent that once complete, process 1700 will repeat again for additional iterative parameters until the system has (but is not limited to) solved for a local optimum of the optimal power flow problem. In various embodiments, a coordinator controller similar to coordinator controller 300 illustrated in FIG. 3 can implement this process.

It should be appreciated that node controller process 1600 and coordinator process 1700 can run concurrently to solve for optimal power flow. Distributed gradient descent processes will be discussed further in detail below. Although a variety of distributed gradient descent processes have been described with respect to FIG. 16 and FIG. 17, a variety of processes can be utilized to implement distributed gradient descent on nodes and coordinator controllers in a power distribution network as appropriate to the requirements of specific applications in accordance with embodiments of the invention. Pseudocode that can be utilized to implement distributed gradient descent processes utilized in solving for optimal power flow is illustrated in FIG. 18. The illustrated process 5 is similar to those outlined above with respect to FIG. 16 and FIG. 17.

Distributed Line Searches

Line search processes similar to the line search process 2 illustrated in FIG. 10 can also be implemented in a distributed manner, which can give rise to a distributed implementation of the inner loop of process 3 illustrated in FIG. 11.

In process 5 illustrated in FIG. 18, each node controller (agent) i∈⁺ keeps track of its last approved power injection (p_i^old, q_i^old) and proposes tentative power injections (p_i^new, q_i^new) to the coordinator; the coordinator decides whether to approve the tentative power injections.

Tentative power injection (p^new, q^new) is computed by gradient projection

$p_i^{\mathrm{new}}\leftarrow \prod _{\left[{\underset{\_}{p}}_i,{\stackrel{\_}{p}}_i\right]}\left(p_i^{\mathrm{old}}-\eta {\partial }_{p_i}L\right);$ $q_i^{\mathrm{new}}\leftarrow \prod _{\left[{\underset{\_}{q}}_i,{\stackrel{\_}{q}}_i\right]}\left(q_i^{\mathrm{old}}-\eta {\partial }_{q_i}L\right);$ where the gradient is computed as in process 5 and the step size η is controlled by the coordinator. The coordinator initializes η=1 and reduces η by a fraction of 1−α until voltage constraints are satisfied, i.e., ν_i≤ν_i≤ν_i for i∈⁺, and the modified objective function is well-approximated by its linearization, i.e., L^new<L^thres.

The coordinator decides whether to approve the tentative power injection, i.e., set (p^old, q^old)←(p^new, q^new), and when to terminate an inner loop, i.e., set (p′, q′). In other cases, the coordinator reduces the current to η to αη ask for the submission of new tentative power injection (p^new, q^new).

The coordinator makes these decisions based on (Δs_i, ΔL_i) computed by the node controllers (agents) i∈⁺. The quantity Δs_i captures the update size of (p_i, q_i), and ΔL_i is the product of gradient ∂_{(pi, qi)}L and power injection update. In particular, if Δs_i≤ε, the coordinator decides to stop the inner loop; else if L^new≤L^old+βL_i, the coordinator decides to approve the tentative power injection, i.e., (p^old, q^old)←(p^new, q^new); else the coordinator reduces the step size η by a fraction of 1−α, i.e., η←αη.

Multiphase Networks

In various embodiments, gradient projection processes can also be applied to multiphase networks. Solving for OPF for multiphase networks using distributed gradient projection processes will be discussed further below.

In several embodiments, the objective function to minimize is

$L\left(x;\mu \right)=\sum _{\varphi \in {\Phi }_0}{a\left(p_0^{\varphi }\left(x\right)\right)}^2+{\mathrm{bp}}_0^{\varphi }\left(x\right)-\sum _{k=1}^n\sum _{ \varphi \in {\Phi }_k}\left(\underset{\_}{\mu }\ln \left(v_k^{\varphi }-{\underset{\_}{v}}_k\right)+\stackrel{\_}{\mu }\ln \left({\stackrel{\_}{v}}_k-v_k^{\varphi }\right)\right).$ Hence, to compute ∂_xL, it suffices to compute ∂_xp₀^ϕ and ∂_xν_k^ϕ for each k∈⁺ and each ϕ∈{a, b, c}.

∂_xp₀^ϕ and ∂_xν_k^ϕ can be estimated for each k and each k∈⁺ and each ϕ∈{a, b, c} using the linearized power flow equations

$\sum _{i:i\to j}{\Lambda }_{\mathrm{ij}}=s_j+\mathrm{diag}\left(\mathrm{Diag}\left(v_j\right){\gamma }^{{\Phi }_j}{y}_j^H\right)+\sum _{k:j\to k}{\Lambda }_{\mathrm{jk}}^{{\Phi }_j},j\in 𝒩;$ $S_{\mathrm{ij}}={\gamma }^{{\Phi }_{\mathrm{ij}}}\mathrm{Diag}\left({\Lambda }_{\mathrm{ij}}\right),i\to j;$ $v_j=v_i^{{\Phi }_{\mathrm{ij}}}-\mathrm{diag}\left(S_{\mathrm{ij}}{z}_{\mathrm{ij}}^H+z_{\mathrm{ij}}{S}_{\mathrm{ij}}^H\right),i\to j$ $\mathrm{where}$ $\alpha =e^{-2\pi /3},\beta =\left[\begin{array}{c}1\\ \alpha \\ {\alpha }^2\end{array}\right],\gamma ={\mathrm{\beta \beta }}^H.$

The above equations can be written in scalar format to obtain

$\begin{array}{cc}\sum _{i:i\to j}{\Lambda }_{\mathrm{ij}}^{\varphi }=s_j^{\varphi }+\sum _{k:j\to k}{\Lambda }_{\mathrm{jk}}^{\varphi }+v_j^{\varphi }\sum _{\phi \in {\Phi }_j}{\alpha }^{\varphi -\phi }{\stackrel{\_}{y}}_j^{\mathrm{\varphi \phi }},j\in 𝒩,\varphi \in {\Phi }_j;& \left(10a\right)\\ v_j^{\varphi }=v_i^{\varphi }-\sum _{\phi \in {\Phi }_j}\left({\Lambda }_{\mathrm{ij}}^{\phi }{\alpha }^{\varphi -\phi }{\stackrel{\_}{\stackrel{\_}{z}}}_{\mathrm{ij}}^{\mathrm{\varphi \phi }}+{\Lambda }_{\mathrm{ij}}^{-\phi }{\alpha }^{\phi -\varphi }{z}_{\mathrm{ij}}^{\mathrm{\varphi \phi }}\right),i\to j,\varphi \in {\Phi }_{\mathrm{ij}}.& \left(10b\right)\end{array}$

Additionally, ∂_x(Λ_ij^ϕ, ν_i^ϕ) can be estimated in two rounds as follows. In the first round, the V_j^ϕ . . . term can be ignored in (10a) to obtain∂_piϕΛ_kl^φ=_{ϕ=φi∈Down(l)}, ∂_qiφΛ_kl^ϕ=_{ϕ=φi∈Down(l)}.

This expression can be used to obtain the estimates of ∂_xν_k^ϕ as

$\begin{array}{c}{\partial }_x{v}_k^{\varphi }={\partial }_x{v}_j^{\varphi }-\sum _{\phi \in {\Phi }_k}\left({\partial }_x{\Lambda }_{\mathrm{jk}}^{\phi }{\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{jk}}^{\mathrm{\varphi \phi }}+{\partial }_x{\stackrel{\_}{\Lambda }}_{\mathrm{jk}}^{\phi }{\alpha }^{\phi -\varphi }{z}_{\mathrm{jk}}^{\mathrm{\varphi \phi }}\right)\\ =-\sum _{\left(i,j\right)\in {𝒫}_k}\sum _{\phi \in {\Phi }_j}\left({\partial }_x{\Lambda }_{\mathrm{ij}}^{\phi }{\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{ij}}^{\mathrm{\varphi \phi }}+{\partial }_x{\stackrel{\_}{\Lambda }}_{\mathrm{ij}}^{\phi }{\alpha }^{\phi -\varphi }{z}_{\mathrm{ij}}^{\mathrm{\varphi \phi }}\right)\end{array}$ which simplifies to

$\begin{array}{cc}{\partial }_{p_i^{\phi }}{v}_k^{\varphi }=-2\sum _{\left(s,t\right)\in {𝒫}_{k\bigwedge i}}\mathrm{Re}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right),& \left(11a\right)\\ {\partial }_{q_i^{\phi }}{v}_k^{\varphi }=2\sum _{\left(s,t\right)\in {𝒫}_{k\bigwedge i}}\mathrm{Im}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)& \left(11b\right)\end{array}$ for k∈⁺, ϕ∈Φ_k, and x=p_i^φ or q_i^φ for i∈⁺ and φ∈Φ_i. At last, we estimate ∂_xΛ_kl^ϕ as

$\begin{array}{cc}\begin{array}{c}{\partial }_x{\Lambda }_{\mathrm{kl}}^{\varphi }={\partial }_x{s}_l^{\varphi }+\sum _{m:l\to m}{\partial }_x{\Lambda }_{lm}^{\varphi }+{\partial }_x{v}_l^{\varphi }\sum _{\phi \in {\Phi }_l}{\alpha }^{\varphi -\phi }{\stackrel{\_}{y}}_l^{\mathrm{\varphi \phi }}\\ =\sum _{t\in \mathrm{Down}\left(l\right)}\left[{\partial }_x{s}_t^{\varphi }+{\partial }_x{v}_t^{\varphi }\sum _{\xi \in {\Phi }_t}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_t^{\mathrm{\varphi \xi }}\right],\end{array}k\to l,\varphi \in {\Phi }_l.& \left(12\right)\end{array}$

This can be used to obtain

${\partial }_{p_i^{\phi }}{p}_0^{\varphi }={𝕝}_{\varphi =\phi }+\sum _{t=1}^n{\partial }_{p_i^{\phi }}{v}_t^{\varphi }\mathrm{Re}\left(\sum _{\xi \in {\Phi }_t}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_t^{\mathrm{\varphi \xi }}\right),{\partial }_{q_i^{\phi }}{p}_0^{\varphi }=\sum _{t=1}^n{\partial }_{q_i^{\phi }}\mathrm{Re}\left(\sum _{\xi \in {\Phi }_t}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_t^{\mathrm{\varphi \xi }}\right)$ $\mathrm{for}\varphi \in {\Phi }_0,i\in {𝒩}^{+},\mathrm{and}\phi \in {\Phi }_i.$

Furthermore, the gradient ∂_xL can be estimated as follows:

$\begin{array}{c}{\partial }_{p_i^{\phi }}L=\sum _{\varphi \in {\Phi }_0}\left(2{\mathrm{ap}}_0^{\varphi }+b\right){\partial }_{p_i^{\phi }}{p}_0^{\varphi }-\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}-\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]{\partial }_{p_i^{\phi }}{v}_k^{\varphi }\\ =\sum _{\varphi \in {\Phi }_0}\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\left[{𝕝}_{\varphi =\phi }+\sum _{k=1}^n{\partial }_{p_i^{\phi }}{v}_k^{\varphi }\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)\right]-\\ \sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}-\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]{\partial }_{p_i^{\phi }}{v}_k^{\varphi }\\ =2{\mathrm{ap}}_0^{\phi }+b+\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\\ \frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}]{\partial }_{p_i^{\phi }}{v}_k^{\varphi }\\ =2{\mathrm{ap}}_0^{\phi }+b-2\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\\ \frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}]\sum _{\left(s,t\right)\in {𝒫}_{k\bigwedge i}}\mathrm{Re}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\\ =2{\mathrm{ap}}_0^{\phi }+b-2\sum _{\left(s,t\right)\in {𝒫}_i}\sum _{\varphi \in {\Phi }_t}\mathrm{Re}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\sum _{k\in \mathrm{Down}\left(t\right)}\\ \left[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]\end{array}$ $\mathrm{and}$ $\begin{array}{c}{\partial }_{p_i^{\phi }}L=\sum _{\varphi \in {\Phi }_0}\left(2{\mathrm{ap}}_0^{\varphi }+b\right){\partial }_{q_i^{\phi }}{p}_0^{\varphi }-\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}-\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]{\partial }_{q_i^{\phi }}{v}_k^{\varphi }\\ =\sum _{\varphi \in {\Phi }_0}\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\sum _{k=1}^n{\partial }_{q_i^{\phi }}{v}_k^{\varphi }\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\\ \sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}-\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]{\partial }_{q_i^{\phi }}{v}_k^{\varphi }\\ =\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]\\ {\partial }_{q_i^{\phi }}{v}_k^{\varphi }\\ =2\sum _{k=1}^n\sum _{\varphi \in {\Phi }_k}\left[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]\\ \sum _{\left(s,t\right)\in {𝒫}_{k\bigwedge i}}\mathrm{Im}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\\ =2\sum _{\left(s,t\right)\in {𝒫}_i}\sum _{\varphi \in {\Phi }_t}\mathrm{Im}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\sum _{k\in \mathrm{Down}\left(t\right)}\\ \left[\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\mathrm{\varphi \xi }}\right)-\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}\right]\end{array}$ $\mathrm{for}$ $i\in {𝒩}^{+}$ $\mathrm{and}$ $\phi \in {\Phi }_i.$

In many embodiments, g_k^ϕ can be defined as

$g_k^{\varphi }\text{:}=\left(2{\mathrm{ap}}_0^{\varphi }+b\right)\mathrm{Re}\left(\sum _{\xi \in {\Phi }_k}{\alpha }^{\varphi -\xi }{\stackrel{\_}{y}}_k^{\varphi \xi }\right)-\frac{\underset{\_}{\mu }}{v_k^{\varphi }-{\underset{\_}{v}}_k}+\frac{\stackrel{\_}{\mu }}{{\stackrel{\_}{v}}_k-v_k^{\varphi }}$ $\mathrm{for}k\in {𝒩}^{+}$ $\mathrm{and}$ $\varphi \in {\Phi }_k,\mathrm{then}$ ${\partial }_{p_i^{\phi }}L=2{\mathrm{ap}}_0^{\phi }+b-2\sum _{\left(s,t\right)\in {𝒫}_i}\sum _{\varphi \in {\Phi }_t}\mathrm{Re}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\sum _{k\in \mathrm{Down}\left(t\right)}{g}_k^{\varphi },{\partial }_{q_i^{\phi }}L=2\sum _{\left(s,t\right)\in {𝒫}_i}\sum _{\varphi \in {\Phi }_t}\mathrm{Im}\left({\alpha }^{\varphi -\phi }{\stackrel{\_}{z}}_{\mathrm{st}}^{\mathrm{\varphi \phi }}\right)\sum _{k\in \mathrm{Down}\left(t\right)}{g}_k^{\varphi }$ $\mathrm{for}$ $i\in {𝒩}^{+}$ $\mathrm{and}$ $\phi \in {\Phi }_i.$

g_k can be written in a more compact form:g_k=Re[Diag(β^Φk)·y_k·β^Φk][2ap₀+b]−μ(ν−ν)+μ(ν−ν)for k∈⁺, where denotes pair-wise multiplication, and denotes component-wise division, i.e., for a, b∈ⁿ and c∈,(ab)_k=a_k·b_k,(ca)_k=c/a_k, k=1,2, . . . ,n.

Then, to compute ∂_piL and ∂_qiL, first do a backward sweep to compute sumDown_k and then do a forward sweep to compute sumUp_i, more explicitly as follows:

${\mathrm{sumDown}}_t=\sum _{k\in \mathrm{Down}\left(t\right)}{g}_k^{{\Phi }_t},t\in {𝒩}^{+};$ $h_{\left(s,t\right)}=\mathrm{Diag}\left({\beta }^{-{\Phi }_{\mathrm{st}}}\right)·z_{\mathrm{st}}^H·\mathrm{Diag}\left({\beta }^{{\Phi }_{\mathrm{st}}}\right)·{\mathrm{sumDown}}_t,s\to t;$ ${\mathrm{sumUp}}_i=\sum _{\left(s,t\right)\in {\Phi }_i}{h}_{\left(s,t\right)}^{{\Phi }_i},i\in {𝒩}^{+}.$

The backward-forward procedure woks as

${\mathrm{sumDown}}_t=g_t+\sum _{u:t\to u}{\mathrm{sumDown}}_u^{{\Phi }_t},t\in {𝒩}^{+};$ ${\mathrm{sumUp}}_i={\mathrm{sumUp}}_h^{{\Phi }_i}+h_{\left(h,i\right)},h\to i.$

At last, the expressions for ∂_piL and ∂_qiL is as follows:∂_piL=2ap₀^Φi+b−2Re(sumUp_i),∂_qiL=2Im(sumUp_i)for i∈⁺.Simulated Numerical Results

The accuracy and efficiency of a process similar to process 3 illustrated in FIG. 11 can be evaluated for a number of simulated balanced networks. In particular, the convex relaxation approach can be used to obtain the global optimal value of (1). The amount the objective value obtained by process 3 deviates from the global optimal value can be checked. The convex relaxation is solved by CVX, a convex relaxation approach, and its execution time is used as a benchmark to investigate the simulated efficiency of process 3. All simulations use Matlab 7.9.0.529 (64-bit) with toolbox cvx 1.21 on Mac OS X 103.5 with 2.66 GHz Intel Core 2 Due CPU and 4 GB 1067 MHz DDR3 memory.

Simulated Test Networks

The simulated test networks include a 47-node (47-bus), a 56-node (56-bus) network, and subnetworks of a 2065-node (2065-bus) network. These networks are all in the service territory of Southern California Edison (SCE), a utility company in California, USA. Topologies of the 47-node (47-bus) network and the 56-node (56-bus) network are illustrated in FIGS. 19A and 19B respectively. Table I illustrated in FIG. 20 shows parameters of a simulated 47-node (47-bus) network including line impedances, peak spot load, and nameplate ratings of capacitors and PV generators. Table II illustrated in FIG. 21 shows parameters of a simulated 56-node (56-bus) network including Line impedances, peak spot load, and nameplate ratings of capacitors and PV generators.

Simulated OPF Setup

The following OPF setup is used throughout this section.

• 1. The objective is to minimize power loss in the network,
• 2. The power injection constraints are as follows. For each node (bus) i∈, there may be multiple devices including loads, capacitors, and PV panels. Assume that there is a total of D_i such devices and label them by 1, 2, . . . , D_i. Let s_i,d=p_i,d+iq_i,d denote the power injection of device d for d=1, 2, . . . , D_i. If device d is a load with given real and reactive power consumptions p and q, then we imposes_i,d=−p−iq.  (13)If device d is a load with given peak apparent power s_peak, then we imposes_i,d=−s_peakexp(jθ)  (14)where θ=cos⁻¹(0.9), i.e, power injection s_i,d is considered to be a constant, obtained by assuming a power factor of 0.9 at peak apparent power. If device d is a capacitor with nameplate q, then we imposeRe(s_i,d)=0 and 0≤Im(s_i,d)≤q.  (15)If device d is a PV panel with nameplate s and real power generation p_i, then we imposeRe(s_i,d)=p_i and |s_i,d|≤s.  (16)The power injection at node (bus) i is

$s_i=\sum _{d=1}^{D_i}{s}_{i,d}$

where s_i,d satisfies one of (13)-(16).

• 3. The voltage regulation constraint is considered to be 0.95²≤ν_i≤1.05² for i∈⁺.Simulated Results

Numerical results are summarized in Table III illustrated in FIG. 22. Table III shows objective values and CPU times of CVX and IPM. It can be seen that process 3 obtains 70× speed up over the CVX/sedumi approach, at the cost of around 2e-7 suboptimality ratio for large-scale networks. Note that process 3 is run in series rather than parallel due to the limitation of platform on which the simulation was performed. The speed up can be even more significant when parallel implementation is utilized.

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 may be practiced otherwise than specifically described. Thus, embodiments of the present invention should be considered in all respects as illustrative and not restrictive.

Claims

1. A node controller comprising: a network interface:a processor;a memory containing: a node controller application;a plurality of node operating parameters describing operating parameters of a node; anda plurality of node operating parameters describing operating parameters for a set of at least one connected node selected from a group consisting of at least one downstream node and at least one upstream node, where the node is connected to the set of at least one connected node in a power distribution network;wherein the processor is directed by the node controller application to: receive a plurality of coordinator parameters describing operating parameters of a node coordinator by the network interface;calculate a plurality of updated node operating parameters of the node using an iterative gradient projection process using the plurality of node operating parameters describing operating parameters of the node and the plurality of node operating parameters describing operating parameters for the set of at least one connected node, wherein: the iterative gradient projection process comprises a backwards sweep process comprising: receiving operating parameters from the one or more downstream nodes;calculating a plurality of updated node operating parameters using the operating parameters from the one or more downstream nodes; andsending the plurality of updated node operating parameters to the one or more upstream nodes;each iteration in the iterative process is determined by the coordinator parameters; andthe node operating parameters comprise power injection parameters calculated based on the operating parameters, the plurality of node operating parameters, and the coordinator parameters; andadjust a power flow in the power distribution network by controlling the plurality of node operating parameters describing operating parameters of the node based on the plurality of updated node operating parameters of the node when the power injection parameters are within a specified range.

2. The node controller of claim 1, wherein the iterative gradient projection process is a distributed process.

3. The node controller of claim 1, wherein the iterative gradient projection process further comprises calculating gradient parameters.

4. The node controller of claim 1, wherein the iterative gradient projection process further comprises calculating a gradient step size.

5. The node controller of claim 1, wherein the iterative gradient projection process further comprises a forward sweep process.

6. A node controller comprising: a network interface;a processor;a memory containing: a node controller application,a plurality of node operating parameters describing operating parameters of a node; anda plurality of node operating parameters describing operating parameters for a set of at least one connected node selected from a group consisting of at least one downstream node and at least one upstream node, where the node is connected to the set of at least one connected node in a power distribution network;wherein the processor is directed by the node controller application to: receive a plurality of coordinator parameters describing operating parameters of a node coordinator by the network interface;calculate a plurality of updated node operating parameters of the node using an iterative gradient projection process using the plurality of node operating parameters describing operating parameters of the node and the plurality of node operating parameters describing operating parameters for the set of at least one connected node, wherein: the iterative gradient projection process comprises a forward sweep process comprising: receiving operating parameters from the one or more upstream nodes;calculating a plurality of updated node operating parameters using the operating parameters from the one or more upstream nodes; andsending the plurality of updated node operating parameters to the one or more downstream nodes;each iteration in the iterative process is determined by the coordinator parameters; andthe node operating parameters comprise power injection parameters calculated based on the operating parameters, the plurality of node operating parameters and the coordinator parameters; andadjust a power flow in the power distribution network by controlling the plurality of node operating parameters describing operating parameters of the node based on the plurality of updated node operating parameters of the node when the power injection parameters are within a specified range.

7. The node controller of claim 3, wherein calculating gradient parameters further comprises exact calculation of the gradient parameters.

8. The node controller of claim 3, wherein calculating gradient parameters further comprises an approximation of the gradient parameters.

9. The node controller of claim 3, wherein calculating gradient parameters is a distributed process.

10. The node controller of claim 4, wherein calculating a gradient step size is evaluated using a line search.

11. The node controller of claim 1, wherein the node is part of a radial network topology.

12. The node controller of claim 1, wherein node operating parameters include power injection, current, and impedance.

13. The node controller of claim 1, wherein the node is configured to control operating parameters as components in a single phase power distribution network.

14. The node controller of claim 1, wherein the node is configured to control operating parameters as components in a multiphase balanced power distribution network.

15. The node controller of claim 1, wherein the node is configured to control operating parameters as components in a multiphase unbalanced power distribution network.