Method for optimizing power flows in electric power networks

Abstract

Power flow in an electric power network is optimized by first decomposing an optimization problem into a set of disjoint parameterized optimization problems. The disjoint optimization problems are independent of each other, and the decomposition is based on dualized coupled constraints having corresponding multipliers. Each optimization problem is solved independently to obtain a corresponding solution, and a sensitivity of each solution to changes in the parameter. The parameters are updated using the corresponding solutions and the sensitivities, and iterated until reaching a convergence.

Description

FIELD OF THE INVENTION

The present invention relates generally to electric power networks, and more particularly to optimizing the power flows in the networks.

BACKGROUND OF THE INVENTION

An electric power network includes buses connected to transmission lines. The buses are locally connected to generators and loads. Optimal power flow (OPF) analysis is often used for monitoring and controlling the operation of the network. The power flow depends, in part, on voltage magnitudes and phase angles. Power flows and voltage levels on the buses are optimized by minimizing an objective function subject to constraints, such as the magnitudes, phases, power transferred, generator capacity, thermal losses, and the like.

Most conventional OPF optimizations:

• • (i) Solve the problem in a centralized manner.
  • (ii) Do not exploit a topology of the network to distribute the optimization over the buses.

Some conventional methods for distributing the optimization problem:

• • (i) Use the Schur complement method, also called iterative substructuring, to distribute the computation of a single step in the optimization problem.
  • (ii) When the optimization problem is distributed by dual decomposition, a conventional subgradient method is used for solving the optimization problem. The subgradient method converges sub-linearly to a solution.

Thus, there remains a need to optimize power flows in electric power networks in an efficient and expedient manner by appropriately distributing the computations.

U.S. Pat. No. 6,625,520 describes a system and method for operating an electric power system that determines optimal power flow and available transfer capability of the electric power system based on the optimal power flow. The system derives data associated with an initial phase angle and maximum electric power value of a generator by determining mechanical output and electrical output of a generator, including a generator phase angle defined by a time function with a constraint condition that the generator phase angle does not exceed a preset value.

SUMMARY OF THE INVENTION

The embodiments of the invention provide methods for optimizing power flows in electric power networks using a decomposition and coordination procedure. The decomposition procedure distributes the optimization problem into a set of smaller disjoint parameterized optimization problems that are independent of each other. The coordination procedure modifies the parameter associated with the individual problems to ensure that a solution of the entire problem is attained.

The methods are based on dualizing coupled constraints to obtain the set of smaller decoupled optimization problems. In one embodiment of the method, the theory of semi-smooth equations is used in the coordination procedure. The semismooth equation theory ensures that superlinear convergence can be theoretically guaranteed in a neighborhood of the solution. Further, the theory allows for using a merit function to ensure global convergence to a solution using initial parameters that are not near the solution.

In one embodiment, the theory of smoothing based methods is used to solve the decomposed problems. A monotonic decrease of the smoothing parameter is used to ensure that superlinear convergence can be theoretically guaranteed in the neighborhood of a solution. Further, the theory allows for using a merit function to ensure global convergence to a solution even when the initial parameters are far from optimal.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 is a schematic of an electric power network and a graph representing the network that can be used by embodiments of the invention;

FIG. 2 is a flow chart of a conventional power optimization method;

FIG. 3 is a flow chart of a power optimization method according to embodiments of the invention; and

FIGS. 4 and 5 are pseudo code for power optimization methods according to some embodiments of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Electrical Power Network Topology and Representative Graph

FIG. 1 shows a topology of an electric power network that can use embodiments of the invention. The network can include AC components, and DC components connected by convertors. The only requirement is that variables and constraints that control the operation of the network are continuously controllable.

The network includes buses 10 locally connected to loads (L) 12 and generators (G) 14. The buses are interconnected by transmission lines 20, also known as branches (B). Some of the transmission lines can be connected to transformers (T) 22.

The generators supply active power (measured in, e.g., Mega Watts (MW)), and reactive power (measured in Mega Volt Ampere Reactive (MVar)). The loads consume the power. The power is defined by voltage magnitude and phase angle.

The parameters for the optimization include, but are not limited to, an admittance matrix based on the branch impedance and bus fixed shunt admittance, and the flow capacity ratings, i.e., the maximal total power flow constrained by thermal ratings.

The topology of the network can be represented by a graph G of nodes (generators and connected loads) 30 connected by edges (transmission lines) 31.

Input

Input to the method optimization includes the following:

• 1) A graph G(N, E) with a set of N nodes connected by a set of E edges (i, j).
• 2) An admittance of the lines y_ij=g_ij+jb_ij∀(i, j)∈E, where g represents conductance of the line, b represents susceptance (imaginary part of the admittance) of the line with j=√{square root over (−1)}.
• 3) Constraints on active power p_i^G,min,P_i^G,maxi∈N that can be produced by the generators, and the reactive power Q_i^G,min,Q_i^G,max∀i∈N that can be produced by the generators.
• 4) Constraints S_ij^max,P_ij^max∀(i, j)∈E on apparent and active power transferred on the lines.
• 5) Limits V_i^mix,V_i^max∀i∈N on voltage magnitudes at the bases.
• 6) Constraints L_ij^max∀(i, j)∈E on thermal losses on the lines.

Output

Output of the method includes complex valued voltages V_i∀i∈N at the buses, and active and reactive power levels P_i^G,Q_i^G∀i∈N of the generators.

The optimization uses a decision function f(P^G,Q^G,V) that depends on active power generation variables P^G=(P₁^G, . . . , P_|N|^G), reactive power generation variables Q^G=(Q₁^G, . . . , Q_|N|^G), and the complex valued voltages V=(V₁, . . . , V_|N|) at the buses.

Optimal Power Flow

The form of the decision function f is quadratic and strictly increasing:

$f\left(P^G,Q^G,V\right)=\sum _{i\in N}\left({c_{2i}\left(P_i^G\right)}^2+c_{li}{P}_i^G+c_{0i}\right),$ where c indicates constants, with c_2i,c_1i≧0∀i∈N.

To model the limits of feasible operation of the network, equality constraints, inequality constraints, and bounds on the decision variables are used. The operation of the electrical network can be modeled by the equality constraintsh_n(P^G,Q^G,V)=0∀n=1, . . . , Ne,where Ne indicates the number of equality constraints.

The constraints on the limits on power generated by generators, limits on voltage magnitudes at the buses, power transferred on the lines and thermal losses ensuring feasible operation are modeled as inequality constraintsg_n(P^G,Q^G,V)≦0∀n=1, . . . , N_i,where N_i indicates the number of inequality constraints.

To determine the voltages at the buses and the powers produced by the generators, the following optimization problem is solved for global optimality:minimize f(P^G,Q^G,V)subject to h_n(P^G,Q^G,V)=0∀n=1, . . . , Neg_n(P^G,Q^G,V)≦0∀n=1, . . . , Ni  (1)where h_n represents equality constraints, and g_n represents equality constraints.

Optimal Power Flow-Constraints

In the preferred embodiment, the equality constraintsh_n(P^G,Q^G,V)=0∀n=1, . . . , Neare represented as

• • Power flows on the lines

$\begin{array}{c}{S}_{\mathrm{ij}}=V_i\left(y_{\mathrm{ij}}\left(V_i-V_j\right)\right)*\\ S_{ji}=V_j\left(y_{\mathrm{ij}}\left(V_j-V_i\right)\right)*\end{array}\}\forall \left(i,j\right)\in E$

• • Power balances at the buses

$\sum _{j\in N\left(i\right)}{S}_{ij}=S_i^G-S_i^D\forall i\in N,$ where S_ij=P_ij+jQ_ij denotes the complex valued power transferred from bus i to bus j, S_ji=P_ji+jQ_ji denotes the complex valued power transferred from bus j to bus i, (V_i)* denotes the complex conjugate of the complex valued variable, S_i^G=P_i^G+jQ_i^G denotes the complex valued power produced by the generators and S_i^D=P_i^D+jQ_i^D denotes the complex valued power demands. The variables representing power flow on the lines are used for convenience.

In the preferred embodiment, the inequality constraintsg_n(P^G,Q^G,V)=0∀n=1, . . . , Niare represented as

• • Limit on apparent power transferred on lines

$\begin{array}{c}S_{\mathrm{ij}}\le S_{\mathrm{ij}}^{\max }\\ S_{\mathrm{ji}}\le S_{\mathrm{ij}}^{\max }\end{array}\}\forall \left(i,j\right)\in E$

• • Limit on active power transferred on lines

$\begin{array}{c}\mathrm{Re}\left(S_{\mathrm{ij}}\right)\le P_{\mathrm{ij}}^{\max }\\ \mathrm{Re}\left(S_{\mathrm{ji}}\right)\le P_{\mathrm{ij}}^{\max }\end{array}\}\forall \left(i,j\right)\in E$

• • Limit on thermal loss on linesRe(S_ij+S_ji)≦L_ij^max∀(i, j)∈E
  • Limit of Power GenerationP_i^G,min≦P_i^G≦P_i^G,max,Q_i^G,min≦Q_i^G≦Q_i^G,max∀i∈N
  • Limit on Voltage MagnitudeV_i^min≦√{square root over (Re(V_i)²+Im(V_i)²)}{square root over (Re(V_i)²+Im(V_i)²)}≦V_i^max∀i∈N.

Conventional Dual Decomposition Based Optimization

FIG. 2 shows a conventional method for solving an optimization problem 201 as defined by equation (2) below:min f₁(x₁)+f₂(x₂)h₁(x₁)=0subject to h₂(x₂)=0A₁x₁+A₂x₂=bx₁,x₂≧0.

Observe that the optimization problem can be decomposed to separate x₁,x₂, but for the equality constraint A₁x₁+A₂x₂=b. These coupling constraints can be removed by dualizing 210 the constraints in the objective function using multipliers ξ as in equation (3) below:min f₁(x₁)+f₂(x₂)+ξ^T(A₁x₁+A₂x₂−b)h₁(x₁)=0subject to h₂(x₂)=0x₁,x₂≧0.

The optimization problem in equation (3) can be decomposed 220 to separate the variables x₁ 221 and x₂ 222 into a set of disjoint parametrized optimization problems, and corresponding constraints, and then solve the set optimization problems x₁,x₂, independently. In the examples used herein the set includes two disjoint problems x₁ and x₂. However, it is understood, that each problem can be further decomposed to a liner granularity as necessary.

Denote by (ξ) 231 and (ξ) 232 the optimal solution to each of the decomposed problems. The correct choice of multipliers ξ is rarely known perfectly. As a consequence, the method iterates 270 until convergence 250 using a convergence test 240, and otherwise updating 260 the multipliers ξ using only (ξ),(ξ), until the constraints A₁(ξ)+A₂(ξ)=b are approximately satisfied.

The updating can be performed using a line search procedure to find a local minimum. The line search finds a descent direction along which the objective function is reduced, and a step size. The descent direction can be computed by various methods, such as gradient descent, Newton's method and the quasi-Newton method.

Alternatively, a trust region can be used to optimizate a subset of the region a (quadratic) model function. If an adequate model of the objective function is found within the trust region, then the region is expanded, otherwise the region is contracted.

Fast Dual Decomposition Based Optimization

FIG. 3 shows an embodiment that accelerates the convergence of dual decomposition based methods. The initial steps are as described above for the conventional method to obtain a set of disjoint optimization problems.

The key difference is as follows. When solving each of the optimization, the invented method obtains the optimal solution (ξ),(ξ) 231-232 as before, but also sensitivities of the solution to the choice of ξ, ∇_ξ(ξ),∇_ξ(ξ) 331-332.

The sensitivity measures the variation in the solutions (ξ),(ξ) to changes in the multiplying parameter ξ.

The updating is also different, in that now the updating 360 uses

The sensitivities can be obtained by a solution of the following linear complementarity problem. In the subproblem for x₁ the sensitivities are obtained as:∇_x1x1L₁(,,;ξ)∇_ξ(ξ)+∇_x1h₁()∇_ξ−∇_ξ=−A₁^Th₁()+∇_x1h₁()^T∇_ξ(ξ)=0(ξ)+∇_ξ(ξ)≧0⊥(ξ)+∇_ξ(ξ)≧0.  (4)where λ and ν are respectively the multipliers corresponding to the inequality constraints and bounds, superscript T is the transpose operator and L₁(x₁,λ₁,ν₁;ξ) is a Lagrangian function defined asL₁(x₁,λ₁,ν₁;ξ)=f₁(x₁)+λ₁^Th₁(x₁)−ν₁^Tx₁+ξ^TA₁x₁.

The sensitivity for x₂ can be obtained in a similar manner. The system in equation (4) can be reduced to a solution of linear equations when the solution satisfies: (ξ)+(ξ)>0, also called a strict complementarity.

Using the sensitivity to the multiplier ξ, a search direction dξ can be determined using the linearization of the dualized constraints as:A₁∇_ξ(ξ)dξ+A₂∇_ξ(ξ)dξ=b−A₁(ξ)−A₂(ξ).  (5)

This constitutes taking a Newton-like direction for the dualized constraints, and this is precisely why the convergence is accelerated. The Newton step is known to converge locally superlinearly when in the neighborhood of the solution. The conventional dual decomposition method does not have this rapid local convergence property.

Further, this approach allows one to define a merit functionΦ(ξ)=∥A₁(ξ)+A₂(ξ)−b∥₂²,to measure the progress of the method towards solving the optimization problem 201 of equation (2), where ∥ ∥₂ represent the Euclidean vector norm.

The merit function influences the choice of the next multipliers through the sufficient decrease requirement where the step length α∈(0,1] is selected to satisfy:Φ(ξ+αdξ)≦Φ(ξ)+βαΦ′(ξ;dξ),  (6)where β>0 is usually a small constant, and Φ′(ξ;dξ) is the directional derivative of the merit function at the point ξ along the direction dξ.The directional derivative is mathematically defined by the following limit:

${\varphi }^{\prime }\left(\xi ;d\xi \right)=\underset{t\to 0}{\lim }\frac{\varphi \left(\xi +td\xi \right)-\varphi \left(\xi \right)}{t}.$

Fast Dual Decomposition Based Optimization for Optimal Power Flow

Consider a partition of the edges in the network into Ng smaller set of edges (E¹, E², . . . , E^Ng) where the sets are disjoint and their union is the set of all edges E. The set of buses in edge set E^k is denoted by N^k. Utilizing this set of smaller networks G(N^k, E^k) the optimal power flow problem can be equivalently formulated as described in the following. Further, denote by K_i the set of networks to which the node i belongs.

The objective function for each sub-network can be posed as:

$f\left(P^{G,k},Q^{G,k},V^k\right)=\sum _{i\in N^k}\left({c_{2i}^k\left(P_i^{G,k}\right)}^2+c_{1i}^k{P}_i^{G,k}+c_{0i}^k\right),$ where P_i^G,k,Q_i^G,k,V_i^k denotes the real power, reactive power and voltage at the node i in node set N^k. The constants are selected as:

$c_{2i}^k=\frac{c_{2i}}{n_i},c_{1i}^k=\frac{c_{1i}}{n_i},c_{0i}^k=\frac{c_{0i}}{n_i},$ where n_i is the number of sub-networks k in which the node i occurs.

The constraints modeling the operation of the electrical network by the equality constraints:h_n^k(P^G,k,Q^G,k,V^k)=0∀n=1, . . . , Ne^k,k=1, . . . , Ngwhere Ne^k indicates the number of equality constraints in the sub-network k.

We model the limits on power generated by generators, limits on voltage magnitudes at the buses, constraints on the power transferred on the lines and thermal losses ensuring feasible operation as inequality constraintsg_n^k(P^G,k,Q^G,k,V^k)≦0∀n=1, . . . , Ni^k,k=1, . . . , Ng,where Ni^k indicates the number of inequality constraints in the sub-network k.

Constraints are also imposed on power generation and voltage magnitudes at the buses.

To determine the voltages at the buses and the powers produced by the generators, the following optimization problem is solved to global optimality:

$\begin{array}{cc}\begin{array}{cc}\mathrm{minimize}& \sum _{k=1}^{\mathrm{Ng}}{f}^k\left(P^{G,k},Q^{G,k},V^k\right)\\ \mathrm{subject}\mathrm{to}& h_n^k\left(P^{G,k},Q^{G,k},V^k\right)=0\forall n=1,\dots ,{\mathrm{Ne}}^k,k=1,\dots ,\mathrm{Ng}\\ & g_n^k\left(P^{G,k},Q^{G,k},V^k\right)\le 0\forall n=1,\dots ,{\mathrm{Ni}}^k,k=1,\dots ,\mathrm{Ng}\\ & \begin{array}{c}{P}_i^{G,k}=P_i^{G,l}\\ Q_i^{G,k}=Q_i^{G,l}\\ \mathrm{Re}\left(V_i^k\right)=\mathrm{Re}\left(V_i^l\right)\\ \mathrm{Im}\left(V_i^k\right)=\mathrm{Im}\left(V_i^l\right)\end{array}\}k=\min \left(K_i\right),l\in N_i,k\ne l,i\in N,\end{array}& \left(7\right)\end{array}$ where the last set of constraints equates the generator power and voltages for nodes that are shared by different sub-networks k. The above formulation is identically equivalent to the optimal power flow formulation.

For ease of exposition, the notation x^k=(P^G,k,Q^G,k,V^k) is used in the following and x_i^k=(P_i^G,k,Q_i^G,k,V_i^G,k) for some i∈N^k. With this notation, the problem in equation (2) can be reformulated as:

$\begin{array}{cc}\mathrm{minimize}& \sum _{k=1}^{\mathrm{Ng}}{f}^k\left(x^k\right)\\ \mathrm{subject}\mathrm{to}& h_n^k\left(x^k\right)=0\forall n=1,\dots ,{\mathrm{Ne}}^k,k=1,\dots ,\mathrm{Ng}\\ & g_n^k\left(x^k\right)\le 0\forall n=1,\dots ,{\mathrm{Ni}}^k,k=1,\dots ,\mathrm{Ng}\\ & x_i^k=x_i^l;k=\min \left(K_i\right),l\in N_i,k\ne l,i\in N.\end{array}$

The above problem equation is decomposed into smaller optimization problems by removing the set of the equality constraints in the optimization and replacing it in the objective function as:

$\begin{array}{cc}\mathrm{minimize}& \sum _{k=1}^{\mathrm{Ng}}{f}^k\left(x^k\right)+\sum _{i\in N}\left(\sum _{k=\min \left(\mathrm{Ni}\right),k\ne l}{\left({\xi }_i^{\mathrm{kl}}\right)}^T\left(x_i^k-x_i^l\right)\right)\\ \mathrm{subject}\mathrm{to}& h_n^k\left(x^k\right)=0\forall n=1,\dots ,{\mathrm{Ne}}^k,k=1,\dots ,\mathrm{Ng}\\ & g_n^k\left(x^k\right)\le 0\forall n=1,\dots ,{\mathrm{Ni}}^k,k=1,\dots ,\mathrm{Ng},\end{array}$ where, ξ_i^kl,k=min(K_i),k≠l are the multipliers for equality constraints that equate the copies of the power generation and voltage variables for nodes that are shared across sub-networks. This procedure is called the dualization of the coupling constraints. This renders the optimization problem decoupled by sub-networks.

The optimization for each sub-network is:

$\begin{array}{cc}\begin{array}{cc}\mathrm{minimize}& \sum _{k=1}^{\mathrm{Ng}}{\tilde{f}}^k\left(x^k;\xi \right)\\ \mathrm{subject}\mathrm{to}& h_n^k\left(x^k\right)=0\forall n=1,\dots ,{\mathrm{Ne}}^k\\ & g_n^k\left(x^k\right)\le 0\forall n=1,\dots ,{\mathrm{Ni}}^k,k=1,\end{array}& \left(8\right)\end{array}$ where the objective function for k=min(K_i) is:

${\tilde{f}}^k\left(x^k;\xi \right)=f^k\left(x^k\right)+\sum _{i\in N}\sum _{l\ne \min \left(K_i\right)}{\left({\xi }_i^{\mathrm{kl}}\right)}^T{x}_i^k$ and for l≠min(K_i) is:

${\tilde{f}}^l\left(x^l;\xi \right)=f^l\left(x^l\right)+\sum _{i\in N}{\left({\xi }_i^{\mathrm{kl}}\right)}^T{x}_i^l$

The decomposition step solves each of the optimization problems for the sub-networks k=l, . . . , Ng for a given choice of the multipliers ξ. Denote by (ξ) the optimal solution to the problem corresponding to the subnetwork k. Further, denote by ∇_ξ(ξ) the sensitivity of the optimal solution to the sub-network k for the given choice of multipliers. The sensitivity can be obtained as follows. The first order optimality conditions for the sub-network k is:∇_xkL^k(x^k,λ^j,k,λ^g,k;ξ)=0h_n^k(x^k)=0∀n=1, . . . , Ne^kg_n^k(x^k)'0⊥λ_n^g,k≧0∀n=1, . . . , Ni^k,  (9)where the Lagrangian function is defined as:

$L^k\left(x^k,{\lambda }^{h,k},{\lambda }^{g,k};\xi \right)={\tilde{f}}^k\left(x^k;\xi \right)+\sum _{n=1}^{{\mathrm{Ne}}^k}{\lambda }_n^{h,k}{h}_n^k\left(x^k\right)+\sum _{n=1}^{{\mathrm{Ni}}^k}{\lambda }_n^{g,k}{g}_n^k\left(x^k\right)$ where the notation λ_n^h,k represent the multipliers for the equality constraints h_n^k; λ_n^g,k represent the multipliers for the inequality constraints g_n^k.

A solution (ξ) to the optimal solution to the problem corresponding to the sub-network k. will necessarily satisfy the first order optimality conditions listed above. The Sensitivity of the optimal solution to the multipliers ξ can be obtained by solving the following linear complementarity problem obtained by differentiating the first order conditions with respect to variables x^k and multipliers (λ^h,k,λ^g,k):

$\begin{array}{cc}{\nabla }_{x^k{x}^k}{L}^k\left({\stackrel{⋒}{x}}^k,{\stackrel{⋒}{\lambda }}^{h,k},{\stackrel{⋒}{\lambda }}^{g,k};\xi \right){\nabla }_{\xi }{\stackrel{⋒}{x}}^k+\sum _{n=1}^{{\mathrm{Ne}}^k}{\nabla }_{x^k}{h}_n^k\left({\stackrel{⋒}{x}}^k\right){\nabla }_{\xi }{\stackrel{⋒}{\lambda }}_n^{h,k}+\sum _{n=1}^{{\mathrm{Ni}}^k}{\nabla }_{x^k}{g}_n^k\left({\stackrel{⋒}{x}}^k\right){\nabla }_{\xi }{\stackrel{⋒}{\lambda }}^{g,k}=-{\nabla }_{x^k\xi }{\tilde{f}}^k\left({\stackrel{⋒}{x}}^k;\xi \right){\nabla }_{x^k}{h_n^k\left({\stackrel{⋒}{x}}^k\right)}^T{\nabla }_{\xi }{\stackrel{⋒}{\lambda }}_n^{g,k}=0\forall n=1,\dots ,{\mathrm{Ne}}^k{g}_n^k\left({\stackrel{⋒}{x}}^k\right)+{\nabla }_{x^k}{g_n^k\left({\stackrel{⋒}{x}}^k\right)}^T{\nabla }_{\xi }{\stackrel{⋒}{x}}^k\le 0\perp {\stackrel{⋒}{\lambda }}_n^{g,k}+{\nabla }_{\xi }{\stackrel{⋒}{\lambda }}^{g,k}\ge 0\forall n=1,\dots ,{\mathrm{Ni}}^k.& \left(10\right)\end{array}$

The sensitivity computations are used to compute a search direction dξ for the multipliers by solving the following equations:(ξ)+∇_ξ(ξ)^Tdξ=(ξ)+∇_ξ(ξ)^Tdξ;k=min(K_i),l∈N_i,k≠l,i∈N.  (11)

Using this search direction a new set of multipliers ξ⁺=ξ+αdξ where 0<α≦1 is selected as described below. Define the function Φ(ξ) termed the merit function as

$\begin{array}{cc}\varphi \left(\xi \right):={\left[\begin{array}{c}{\stackrel{⋒}{x}}_i^k\left(\xi \right)-{\stackrel{⋒}{x}}_i^k\left(\xi \right);k=\min \left(K_i\right),l\in N_i,k\ne l,i=1\\ \dots \\ {\stackrel{⋒}{x}}_i^k\left(\xi \right)-{\stackrel{⋒}{x}}_i^k\left(\xi \right);k=\min \left(K_i\right),l\in N_i,k\ne l,i=N\end{array}\right]}_2^2.& \left(12\right)\end{array}$

The merit function is the Euclidean norm of the residual of the dualized constraints that couple the different sub-networks. This measures the degree to which the original problem in equation (1) has been solved. When Φ(ξ)≈0, the method terminates with (ξ) as the solution for equation (1).

The step for multipliers is selected so that the merit function is decreased as prescribed the following condition:Φ(ξ+αdξ)≦Φ(ξ)+βαΦ′(ξ;dξ)  (13)where Φ′(ξ;dξ) is the directional derivative of the merit function at the point ξ along the direction dξ and is mathematically defined by the following limit:

${\varphi }^{\prime }\left(\xi ;d\xi \right)=\underset{t->0}{\lim }\frac{\varphi \left(\xi +td\xi \right)-\varphi \left(\xi \right)}{t}.$

Using the new multiplier ξ⁺, the decomposed optimization in equation (3), sensitivity computation (4) and multiplier step computation (5) are repeated until the merit function Φ(ξ) is close to zero.

A description of the methodic steps is provided in FIG. 4.

Fast Dual Decomposition with Smoothing for Optimal Power Flow

The fast decomposition approach is based on the computation of sensitivities with respect to the multipliers ξ for the dualized constraints. However, the sensitivity may not exist under some conditions. For instance in optimization problem (3) when (ξ)+(ξ)≯0 does not hold. In such instances only a directional derivative can be obtained.

To rectify this situation consider the modification of the problem solved for each of the sub-networks k as follows. Consider modifying the stationary conditions (9) as,∇_xkL^k(x^k,λ^h,k,λ^g,k;ξ)=0h_n^k(x^k)=0∀n=1, . . . , Ne^kψ(−g_n^k(x^k),λ_n^g,k;τ)=0∀n=1, . . . , Ni^k,  (14)where the function ψ(a,b;τ) is a smoothing function for the complementarity constraints satisfying the property that:

• (i) ψ(a,b;0)=0a≧0⊥b≧0
• (ii) ψ(a,b;τ) is continuously differentiable for all τ>0.

The second property ensures that the for all τ>0 the optimization problems are smooth and continuously differentiable. While the first property ensures that one can recover a solution to equation (9) by solving a sequence smoothed problem when τ→0. There exist several choices for such a function:

• 1. ψ(a,b;τ)=a+b−√{square root over (a²+b²+τ²)}
• 2. ψ(a,b;τ)=a+b−√{square root over ((a−b)²+τ²)}
• 3. ψ(a,b;τ)=ab−τ,a,b≧0

In the third choice the nonnegativity of a and b have to enforced explicitly. This is precisely what interior point methods do.

Using this modification the method where solution to original problem (1) is obtained can be restated by first solving for fixed τ>0 the problem in equation (14) to a certain tolerance and then decreasing τ→0 to obtain a solution to equation (1) in the limit.

Denote by ()(ξ;τ) the solutions to equation (14). In the case of smoothing, the sensitivity of solutions to the multipliers ξ can be computed as the solution of the following linear equations:

$\begin{array}{cc}{\nabla }_{x^k{x}^k}{L}^k\left({\stackrel{⋒}{x}}^k,{\stackrel{⋒}{\lambda }}^{h,k},{\stackrel{⋒}{\lambda }}^{g,k};\xi \right){\nabla }_{\xi }{\stackrel{⋒}{x}}^k+\sum _{n=1}^{{\mathrm{Ne}}^k}{\nabla }_{x^k}{h}_n^k\left({\stackrel{⋒}{x}}^k\right){\nabla }_{\xi }{\stackrel{⋒}{\lambda }}_n^{h,k}+\sum _{n=1}^{{\mathrm{Ni}}^k}{\nabla }_{x^k}{g}_n^k\left({\stackrel{⋒}{x}}^k\right){\nabla }_{\xi }{\stackrel{⋒}{\lambda }}^{g,k}=-{\nabla }_{x^k\xi }{\tilde{f}}^k\left({\stackrel{⋒}{x}}^k;\xi \right){\nabla }_{x^k}{h_n^k\left({\stackrel{⋒}{x}}^k\right)}^T{\nabla }_{\xi }{\stackrel{⋒}{x}}^k=0\forall n=1,\dots ,{\mathrm{Ne}}^k\psi \left(-g_n^k\left({\stackrel{⋒}{x}}^k\right),{\stackrel{⋒}{\lambda }}_n^{g,k};\tau \right)+{\nabla }_{x^k}{\psi \left(-g_n^k\left({\stackrel{⋒}{x}}^k\right),{\stackrel{⋒}{\lambda }}_n^{g,k};\tau \right)}^T{\nabla }_{\xi }{\stackrel{⋒}{x}}^k+{\nabla }_{{\lambda }^{g,k}}\psi \left(-g_n^k\left({\stackrel{⋒}{x}}^k\right),{\stackrel{⋒}{\lambda }}_n^{g,k};\tau \right){\nabla }_{\xi }{\stackrel{⋒}{\lambda }}_n^{g,k}=0\forall n=1,\dots ,{\mathrm{Ni}}^k,& \left(15\right)\end{array}$

Observe these are linear equations as opposed to linear complementarity equations as in equation (10).

The sensitivity computations are used to compute a search direction dξ for the multipliers by solving the following equations:(ξ,τ)+∇_ξ(ξ,τ)^Tdξ=(ξ,τ)+∇_ξ(ξ,τ^Tdξk=min(K_i),l∈N_i,k≠l,i∈N  (16)

Using this search direction a new set of multipliers ξ⁺=ξ+αdξ where 0<α≦1 is selected as described below. Define the function Φ(ξ) termed the merit function as

$\begin{array}{cc}\varphi \left(\xi ;\tau \right):={\left[\begin{array}{c}{\stackrel{⋒}{x}}_i^k\left(\xi ;\tau \right)-{\stackrel{⋒}{x}}_i^l\left(\xi ;\tau \right);k=\min \left(K_i\right),l\in N_i,k\ne l,i=1\\ \dots \\ {\stackrel{⋒}{x}}_i^k\left(\xi ;\tau \right)-{\stackrel{⋒}{x}}_i^l\left(\xi ;\tau \right);k=\min \left(K_i\right),l\in N_i,k\ne l,i=N\end{array}\right]}_2^2.& \left(17\right)\end{array}$

The merit function is the Euclidean norm of the residual of the dualized constraints that couple the different sub-networks. This measures the degree to which the original problem (1) has been solved. When Φ(ξ;τ)≈τ, the iterations for particular smoothing parameter value and decrease the parameter can be terminated. When Φ(ξ;0)≈0, terminate the method terminates with (ξ) as the solution for equation (1).

The step for multipliers is selected so that the merit function decreases as prescribed the following condition:Φ(ξ+αdξ;τ)≦Φ(ξ;τ)+βα∇_ξΦ(ξ;τ)^Tdξ,  (18)where ∇_ξΦ(ξ;τ) is the gradient of the merit function at the given point ξ. A full derivative can be defined, instead of a directional derivative as in equation (13). If the smoothing parameter is selected to be decreased superlinearly, then fast local convergence in the neighborhood of the solution to equation can be obtained (1).

FIGS. 4 and 5 show the pseudo code for the above embodiments. The steps of the method can be performed in one or more processors connected to memory and input/output interfaces as known in the art.

Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.

Claims

1. A method for optimizing a power flow in an electric power network by solving an optimization problem, comprising steps of: decomposing the optimization problem into a set of disjoint parameterized optimization problems, wherein the disjoint parameterized optimization problems in the set are independent of each other, and wherein the decomposition is based on constraints having corresponding multipliers, which are parameters of the disjoint parameterized optimization problems, and wherein the constraints are dualized coupled constraints;solving each disjoint parameterized optimization problem independently to obtain a corresponding solution, and a sensitivity of each solution to changes in the parameter associated with said each disjoint parameterized optimization problem;updating the parameters using the corresponding solutions and the sensitivities; andrepeating the decomposing, the solving, and the updating until reaching a convergence, wherein the steps are performed in a processor.

2. The method of claim 1, wherein variables of the constraints are continuously controllable.

3. The method of claim 1, further comprising: representing the electric power network by a graph of nodes and edges, wherein the nodes represent generators and loads connected buses, and the edges (i, j) between the nodes represent transmission lines.

4. The method of claim 1, wherein an input to the method includes: the graph,transmission lines (ij),an admittance of the transmission lines (ij) yij=gij+jbij,

5. The method of claim 4, wherein the optimization uses a decision function f(PG,QG,V) that depends on active power generation variables PG=(P1G, . . . , P|N|G), reactive power generation variables QG=(Q1G, . . . , Q|N|G), and complex valued voltages V=(V1, . . . , V|N|) at the buses.

6. The method of claim 5, wherein the decision function f is quadratic and strictly increasing:

7. The method of claim 1, wherein limits of feasible operation of the electric power network, equality constraints, inequality constraints, and bounds on decision variables are used.

8. The method of claim 1, wherein inequality constraints are represented as a smoothed equality constraints using a smoothing parameter, and the smoothing parameter is reduced to zero during the repeating of the updating.

9. The method of claim 1, wherein the updating of the parameters are performed using a line search procedure.

10. The method of claim 1, wherein the updating of the parameters are performed using a trust region.

11. The method of claim 1, wherein the updating alternates between updating with and without using the sensitivities.

12. The method of claim 1, wherein the sensitivities measures a variation in the solutions to changes in the multipliers.