236 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 33, NO.

1, JANUARY 2018

Feasible Region of Optimal Power Flow:

Characterization and Applications
Hsiao-Dong Chiang, Fellow, IEEE, and Chu-Yang Jiang

Abstract—The feasible region plays a fundamental role in solv- sults on the geometry of a feasible power injection region for
ing optimal power flow (OPF) problems. In this paper, a math- several network topologies of power grids have been developed
ematical characterization of the feasible region is presented. An (see [10], [12], for example).
equivalence is established between the feasible region of an OPF
problem and the union of regular stable equilibrium manifolds OPF solutions do not exist if the corresponding feasible region
of a quotient gradient system (QGS) that is derived from the set is an empty set. Any solution algorithm proposed for solving
of equality and inequality constraints of the OPF problem. It is OPF problems will not converge if the feasible region does not
further shown that the QGS is completely stable and that each exist. The geometry and characteristics of feasible region indeed
trajectory converges to an equilibrium manifold, making the QGS affect the performance of OPF solution methods. A sufficient
trajectories useful in locating feasible OPF solutions. The theoreti-
cal results developed in this paper have been numerically verified in condition for an OPF problem to have multiple local optimal
several OPF problems. Finally, the notion of a local feasible region solutions is when the feasible region is composed of several
is proposed and discussed. disconnected feasible components. Each feasible component
Index Terms—Optimal power flow, feasible region, characteri- contains at least one local optimal solution while one feasible
zation, visualization, complete stability. component contains the global optimal solution; exceptionally
multiple feasible components contain the global optimal solu-
I. INTRODUCTION tions.
PTIMAL power flow (OPF) is a constrained nonlinear An explicit expression for the feasible region of an OPF prob-
O optimization problem that is complex to solve in mul-
tiple aspects. Many deterministic techniques and intelligent
lem is very challenging to obtain, except for 2-bus OPF problems
[13], [14]. In this paper, a characterization of the feasible region
meta-heuristics have been proposed, including the linear pro- of OPF problems is developed. This mathematical characteriza-
gramming approach, the quadratic programming approach, the tion is computable and can lead to the development of numerical
nonlinear programming approach, the Newton method, the inte- methods for computing feasible solutions. It is shown that the
rior point method, genetic algorithms, particle swarm optimiza- feasible region of a general OPF problem equals the union of
tion, and convex relaxation-based methods (see, for example, regular stable equilibrium manifolds of a non-hyperbolic dy-
[1]– [9]). namical system that is constructed from the set of equality and
In the past, however, few research efforts were directed to- inequality constraints of the OPF problem.
ward investigating another fundamental subject in solving the It is further shown that the non-hyperbolic dynamical system
OPF problem: the feasible region. The feasible region of OPF is completely stable. This complete stability property allows one
problems is defined as the set of control variables for which to locate multiple feasible solutions for the OPF problems via
all the equality and inequality constraints are satisfied. There trajectories of the non-hyperbolic dynamical system, which can
are two kinds of constraints in OPF problems: linear constraints be useful for practical applications.
(box constraints) and nonlinear constraints. Because of the non- In summary, this paper develops and presents:
linear constraints, the feasible region is usually non-convex and 1) A mathematical characterization of the feasible region of
disconnected. Indeed, it was observed in numerical studies that general nonlinear OPF problems utilizing a dynamical
the feasible region is non-convex [4], [10], [11]. Analytical re- system that is constructed from the set of equality and
inequality constraints of the OPF problem.
Manuscript received June 28, 2016; revised October 27, 2016, January 20,
2) A global stability analysis of the constructed non-
2017, and March 17, 2017; accepted March 23, 2017. Date of publication April hyperbolic dynamical system that offers a novel and
27, 2017; date of current version December 20, 2017. This work was supported promising trajectory-based way to find feasible solutions.
in part by the Advanced Research Project Agency-Energy, USA, under Award
#DE-FDA-0001002 and in part by the National Natural Science Foundation
However, some trajectories may converge to degenerate
of China (51337007). Paper no. TPWRS-00930-2016. (Corresponding author: stable equilibrium manifolds, which are not feasible com-
Hsiao-Dong Chiang.) ponents. In this case, an effective method to escape from
H.-D. Chiang is with the School of Electrical and Computer Engineering,
Cornell University, Ithaca, NY 14853, USA (e-mail: hc63@cornell.edu).
degenerate stable equilibrium manifolds and move toward
C.-Y. Jiang is with the School of Electrical Engineering and Automation, regular stable equilibrium manifolds is needed.
Tianjin University, Tianjin 300072, China (e-mail: justpayne@tju.edu.cn). 3) The notion of a local feasible region for large-scale OPF
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
problems for visualization of feasible components and
Digital Object Identifier 10.1109/TPWRS.2017.2692268 local optimal solutions.

0885-8950 © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
CHIANG AND JIANG: FEASIBLE REGION OF OPTIMAL POWER FLOW: CHARACTERIZATION AND APPLICATIONS 237

II. THE OPF FEASIBLE REGION where

We consider an OPF problem subject to the following equality CE = (c1 , . . . , c2N B )T : 2N G +2N B −1 → 2N B ,
and inequality constraint functions:
⎧  CI = (c2N B +1 , . . . , c4N G +4N B +2N L )T : 2N G +2N B −1
⎨PG i − PL i − Vi Vj (Gij cos θij + Bij sin θij ) = 0
∈i
j → 4N G +2N B +2N L .
⎩QG i − QL i − Vi Vj (Gij sin θij − Bij cos θij ) = 0
j ∈i We can use some transformation techniques to transform the
inequalities in (5) into equalities. Here, slack variables are added
i ∈ {1, . . . , NB } (1) to the inequalities as follows:
 
PGmiin ≤ PG i ≤ PGmiax CE (u, y(u))
i ∈ {1, . . . , NG } (2) H(x) = =0
G i ≤ QG i ≤ QG i
Qm in m ax
CI (u, y(u)) + S 2
Vim in ≤ Vi ≤ Vim ax i ∈ {1, . . . , NB } (3) x = (u, y(u), s(u, y(u))) ∈ 6N G +4N B +2N L −1 ,
 
 f
Sl  ≤ Slm ax S 2 = (s21 , . . . , s24N G +2N B +2N L )T (6)
l ∈ {1, . . . , NL } (4)
|Slt | ≤ Slm ax Hence, the set of constraint functions of the OPF problem can
be represented as the following equalities:
where NB denotes the number of buses, NG denotes the number
of generators, and NL denotes the number of transmission lines. H(x) = 0, x ∈ n (7)
Gij and Bij are equivalent conductance and susceptance of a
where H = (h1 , . . . , hm )T : n → m , n = 6NG + 4NB +
transmission line from bus i to bus j, respectively. And PL i and
2NL − 1, m = 4NG + 4NB + 2NL .
QL i are the real and reactive power loads at bus i.
Definition 1: (The feasible region)
Variables:
The feasible region of general OPF problems (1)–(4) is the
1) θi : voltage angle at bus i, and θij = θi − θj
set of control variables where all the equality and inequality
(The voltage angle of reference bus θr ef is set as a constant) constraints of the problem are satisfied, i.e.,

2) Vi : voltage magnitude at bus i F R = {u ∈ 2N G −1 : H(x) = H(u, y(u), s(u, y(u))) = 0}.
3) PG i : real power of generators at bus i (8)
4) QG i : reactive power of generators at bus i Because of the nonlinear constraints, the feasible region FR
Functional expressions: is usually complex and may be composed of multiple path-
1) Slf : branch flow at the from bus of line l connected feasible components, such as
2) Slt : branch flow at the to bus of line l k
Equations (1) are known as the power flow equations while FR = F Ri , (9)
Equations (2) are the physical limitations imposed on the amount i=1
of real and reactive power that can be generated. The operat- where F Ri is the ith disjoint path-connected feasible component
ing limits imposed on the voltage magnitudes are expressed and k denotes the total number of feasible components.
in Equations (3) while Equations (4) are the thermal limits of We next study the dimension of the feasible region of OPF
branch flows through transmission lines and transformers. problems (1)–(4).
The variables in the OPF problem can be divided into con- Fact 1: (Dimension of the feasible region)
trol variables (u) and state variables (y). The control variables Assume that 0 is a regular value of the function H(·).
include the real power injection and voltage magnitude of the Then, the feasible component F Ri is a (2NG − 1)-dimensional
P-V bus (generator buses are modeled as P-V buses, except one smooth manifold.
generator bus is chosen as the slack bus), the voltage magni- The definitions for regular value and a smooth manifold can
tude of the slack bus, phase shifters, transformer taps, capacitor be found in [21]. Fact 1 can be easily proven using the regularity
shunts, FACTS, and so on. In this paper, we do not consider the condition and the implicit function theorem. Moreover, by the
discrete control variables, rendering the dimension of control Sard theorem, the set of regular values of H(·) is residual and
variables as NG + (NG − 1) = 2NG − 1. All the other vari- dense in m . This fact asserts that each feasible component is
ables are state variables, which can be expressed as functions generally a (2NG − 1)-dimensional smooth manifold. But the
of the control variables via the power flow equations in (1). regularity assumption can be violated in some extreme cases,
Regarding the functional expressions (Slf , Slt ), they can be ex- like the ill-posed cases in [24]. In these cases, the dimension of
pressed as explicit functions of the state variables and control the feasible region is less than 2NG − 1.
variables. We next study a 2-bus system (see Fig. 1). This case is the
For simplicity, Equations (1)–(4) can be rewritten as follows: same as that in [13], except that bus 2 is assumed to be a local
 voltage-controlled bus. Hence, in this 2-bus system, bus 1 is
CE (u, y(u)) = 0
, u ∈ 2N G −1 , y ∈ 2N B (5) the slack bus, and the voltage magnitude of bus 1 (i.e., V1 ) is a
CI (u, y(u)) ≤ 0 control variable. Bus 2 is a locally controlled bus, so V2 is also a
238 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 33, NO. 1, JANUARY 2018

Fig. 1. A simple two-bus system.

Fig. 3. The feasible region shrinks with tightened voltage constraints.

Definition 2: (Equilibrium manifold)

A path-connected component of F −1 (0) of system (11) is
called an equilibrium manifold of system (11). In brief, every
point, say x, of an equilibrium manifold, satisfies F(x) = 0.
Fig. 2. The feasible region of the 2-bus system is composed of two path- Definition 3: (Pseudo-hyperbolic equilibrium manifold)
connected 1-dimensional manifolds. An equilibrium manifold Σ of (11) is called pseudo-
hyperbolic if, for each x ∈ Σ, the eigenvalues of DF (x) (the
control variable. The constraint set of this simple 2-bus system Jacobian matrix of F (·) at x) whose corresponding eigenvectors
is expressed as: lie on the normal space Nx (Σ) (the orthogonal complement of
⎧ the tangent space Tx (Σ)) have no zero real parts.
⎪ PG 1 − V1 V2 (G12 cos θ12 + B12 sin θ12 ) = 0
⎪
⎪ Definition 4: (Stable equilibrium manifold)
⎪
⎪ −P
⎪ L 2 − V2 V1 (G21 cos θ21 + B21 sin θ21 ) = 0 For a pseudo-hyperbolic equilibrium manifold Σ, it is an
⎪
⎪
⎪
⎪ QG 1 − V1 V2 (G12 sin θ12 − B12 cos θ12 ) = 0 (asymptotically) stable equilibrium manifold if, for each x ∈ Σ,
⎪
⎪
⎪
⎪ all the eigenvalues of DF (x) whose corresponding eigenvectors
⎪
⎪ −QL 2 − V2 V1 (G21 sin θ21 − B21 cos θ21 ) = 0
⎪
⎪ lie on Nx (Σ) have negative real parts; otherwise, it is termed an
⎨ V1m in ≤ V1 ≤ V1m ax
(10) unstable equilibrium manifold.
⎪
⎪ V2m in ≤ V2 ≤ V2m ax We next review the concept of energy function for the non-
⎪
⎪
⎪
⎪ PGm1in ≤ PG 1 ≤ PGm1ax hyperbolic system (11). A comprehensive treatment of energy
⎪
⎪
⎪
⎪ functions for hyperbolic nonlinear dynamical systems can be
⎪
⎪ QmG 1 ≤ QG 1 ≤ QG 1
in m ax
⎪
⎪ found in [15], [17]. A scalar function E : n →  is termed an
⎪
⎪ 2
≤ (S12
m ax 2
⎪
⎪
S12 )
energy function for non-hyperbolic nonlinear system (11) if the
⎩
S21 ≤ (S21
2 m ax 2
) following three conditions are satisfied [17], [18]:
(E1) F (x) = 0 implies that ∇E = 0.
In Fig. 2, the feasible region of (10) is numerically computed (E2) The function along any non-trivial trajectory Φ(t, x0 )
and displayed on the control variable space (i.e., the V1 × V2 d
is non-increasing, i.e., dt E(Φ(t, x0 )) ≤ 0 and the set {t ∈  :
plane). In this case, the feasible region is composed of two path- Ė(Φ(t, x0 )) = 0} has measure zero.
connected manifolds of dimension 1 as asserted by Fact 1. The (E3) If along any trajectory {E(Φ(t, x0 )) : t ≥ 0} is
feasible region shrinks but remains of dimension 1 as the voltage bounded, then the trajectory {Φ(t, x0 ) : t ≥ 0} is bounded.
bounds are tightened from the default interval [0.950, 1.050] per
unit into [0.951, 1.049] per unit, as shown in Fig. 3.
IV. CHARACTERIZATION OF THE OPF FEASIBLE REGION
III. PRELIMINARIES OF NONLINEAR DYNAMICAL SYSTEMS To characterize the feasible region of the set of nonlinear
In this section, some concepts of general nonlinear dynamical constraint functions (1)–(4), we consider the following nonlin-
systems needed in the subsequent analysis are reviewed. ear dynamical system, which is closely related to the nonlinear
A class of non-hyperbolic dynamical systems is defined by constraint functions (7):
the form: ẋ = QH (x) = −DH(x)T H(x) (12)
ẋ = F (x) = M (x)H(x) (11)
where DH(x) is the Jacobian matrix of H(x).
n ×m
where H :  →  and M :  → 
n m n
, m ≤ n. The so- We term the above system the quotient gradient system
lution of system (11), starting from x0 at t = 0, is called a (QGS), which is a nonlinear, non-hyperbolic dynamical sys-
trajectory, denoted by Φ(·, x0 ) :  → n . tem. We will explore the stable equilibrium manifolds (SEMs)
CHIANG AND JIANG: FEASIBLE REGION OF OPTIMAL POWER FLOW: CHARACTERIZATION AND APPLICATIONS 239

of QGS (12) to characterize the feasible region of an OPF prob-

lem. First, we establish the relationship between the OPF feasi-
ble region and the stable equilibrium manifolds of system (12).
Theorem 2: (The feasible region and SEMs)
Let the feasible region of the constraint set in (7) be composed
of several feasible components. Then a feasible component, say
F Ri of constraint set (7) is a stable equilibrium manifold of
QGS (12), say, Σjs for some j, i.e.,
F Ri = Σjs .
In addition, the (entire) feasible region of constraint set (7)
is contained in the union of the stable equilibrium manifolds of
QGS (12), i.e.,
FR ⊂ Σjs .
j=1

Proof: See the Appendix.

We next show that QGS (12) admits an energy function.
Theorem 3: (Energy function)
The quotient gradient system (12) admits the function E(x) =
1 2
2 H(x) as an energy function.
Proof: See the Appendix.
The existence of an energy function plays an important role
in developing analytical results for the QGS (12). Theorem 3
will be explored to derive the asymptotical behaviors of the
trajectories of QGS (12) in Theorem 6.
Theorem 4: (The feasible region and local minima)
If Σs is a stable equilibrium manifold of QGS (12), then Σs
is either a feasible component of constraint set (7) or a nonzero
Fig. 4. The relationship between the feasible region and regular stable equi-
local minimum point of the energy function E(x) = 12 H(x) 2 . librium manifolds of the corresponding QGS. (a) The feasible region of the
Proof: See the Appendix. 2-bus system. The feasible region is composed of two feasible components.
Theorem 2: and Theorem 4 assert that each feasible compo- (b) Stable equilibrium manifolds of the 2-bus system. There are two regular
stable equilibrium manifolds in the corresponding QGS that correspond to each
nent of constraint set (7) corresponds to a SEM of system (12). feasible component.
Therefore, the feasible region of OPF problems can be mathe-
matically characterized by the SEMs of system (12). However,
the SEMs of system (12) may include the so-called degenerate programming problems in [18]. But the results derived in [18]
SEMs defined below. hold under certain conditions and cannot be directly applied to
Definition 5: (The degenerate SEM) the OPF problem. In the OPF problem, the power flow con-
A stable equilibrium manifold Σs of QGS (12) is said to be straints have their own explicit expressions and each variable
degenerate if H(Σs ) = 0 and DH(Σs )T H(Σs ) = 0. has its own physical upper and lower bounds, all of which have
It is obvious that rank(DH(Σs )) < m if Σs is a degenerate been discussed in the proof of Theorem 3. Additionally, Theo-
SEM. In contrast, the non-degenerate (i.e., regular) SEM satis- rem 3 forms the basis of Theorem 4 and Theorem 5. As a result,
fies H(Σs ) = 0. By combining Theorem 2 and Theorem 4, the the conditions stated in [18] are not required for Theorem 4
following characterization of the feasible region is derived. and Theorem 5. Hence, Theorem 5 gives a sharpened result by
Theorem 5: (Characterization of the feasible region) exploring the structure of OPF problems (1)–(4) and derives a
A path-connected set is a feasible component of constraint necessary and sufficient condition.
set (7) if and only if it is a regular stable equilibrium manifold We next employ the 2-bus system as an example to illus-
of QGS (12), i.e., trate Theorem 5. The OPF constraint set of the 2-bus system
is described in (10). The corresponding QGS is constructed as
FR = Σjr s follows:
j=1

where Σjr s
is a regular stable equilibrium manifold of QGS (12). ẋ = −DH2 − bu s (x)T H2 − bu s (x). (13)
Proof: See the Appendix.
Theorem 5 asserts that the feasible region of general OPF The feasible region of (10) and the SEMs of system (13)
problems (1)–(4) equals the union of regular stable equilibrium are numerically computed and compared in Fig. 4. The feasible
manifolds of QGS system (12). We point out that the QGS has region in Fig. 4(a) is computed using the interior point method
been developed to obtain feasible points of general nonlinear with a constant objective function value. Indeed, the feasible
240 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 33, NO. 1, JANUARY 2018

region in this case equals the union of the two SEMs in Fig. 4(b), stability region of the corresponding SEM, which is a
as asserted by Theorem 5. smooth manifold. This is because the vector field of the
QGS (12) is at least continuously differentiable [15].
Hence, the proposed method is not sensitive to initial
V. APPLICATIONS conditions in the sense that the initial points converge to
The proposed QGS nonlinear system (12) can provide novel different points lying on the same SEM as long as the
computational mechanisms to locate a path-connected feasible initial points lie inside the same stability region.
component. A conceptual numerical algorithm will be intro- ii) Since the proposed QGS-based trajectory method does
duced in this section. To shed some light on this mechanism, not require matrix inversion, the ill-conditioning of a Ja-
we next show that the QGS possesses a desirable numerical cobian matrix along a trajectory of QGS (12) is avoided.
property: each trajectory converges to an equilibrium manifold. iii) The QGS-based method provides a novel and promising
Theorem 6: (Complete stability) trajectory-based way to compute a feasible solution, if it
QGS (12) is completely stable, i.e., each trajectory of system exists. This method can be combined with other methods,
(12) converges to an equilibrium manifold. such as the interior point method or other optimization-
Proof: See the Appendix. based methods, to obtain feasible solutions and local op-
Theorem 5 and Theorem 6 lead to the following (conceptual) timal solutions.
QGS-based method for computing multiple feasible solutions
of general OPF problems (1)–(4). VI. NUMERICAL STUDY
Given: The nonlinear constraint set of OPF problems (1)–(4)
In this section, the feasible regions of several test systems and
and an initial guess.
the stable equilibrium manifolds of the corresponding QGS will
Step 1: Construct the corresponding QGS (12).
be simulated to illustrate Theorem 2, Theorem 4, and Theorem
Step 2: Integrate QGS (12), starting from the initial guess
5. In addition, the notion of a local feasible region will be
until it converges to a point lying on one stable equilibrium
proposed and employed to visualize the feasible region of high-
manifold, say x.
dimensional OPF problems.
Step 3: If |H(x)| ≤ ε is a sufficiently small positive scalar,
then x is a feasible solution and go to Step 4; otherwise, try an-
A. The Feasible Region and SEMs
other point that lies outside the stability region of the degenerate
SEM containing x and go to Step 2. The feasible regions of OPF problems can be mathematically
Step 4: If a sufficient number of feasible solutions are ob- characterized by the regular stable equilibrium manifolds of
tained, then stop; otherwise, try another initial point lying in a system (12), as asserted by Theorem 5. To visualize the feasible
neighborhood of x and go to Step 2. region, we sample a variety of initial points to obtain the feasible
Remarks: solutions lying on the stable equilibrium manifolds. To compute
1) Steps 1 through 3 are built on Theorem 6 and on the di- the feasible region and the stable equilibrium manifolds in a
mension property of the stability region and the stability straightforward way, we employ the Latin hypercube sampling
boundary [15]. As asserted by Theorem 6, the state space technique to uniformly sample a large number of points, say
of system (12) is composed of the closure of the stability 10000 points, within the bounds of each variable of the test
regions of SEMs of the system. Therefore, only trajec- systems. We then integrate system (12), starting from these
tories on the stability boundaries converge to unstable initial points, to obtain points lying on the stable equilibrium
equilibrium manifolds. Since the dimension of stability manifolds. As asserted by Theorem 6, the state space of system
regions is n while the stability boundaries are of n−1 (12) is composed of the closure of the stability regions of stable
[15], Theorem 6 implies that almost every trajectory of equilibrium manifolds of the system. Hence, these trajectories
QGS converges to a point lying on one SEM. will converge to either regular or degenerate SEMs. Since the
2) Theorem 6 suggests that computing a feasible point of the Latin hypercube sampling technique is applied to uniformly
general OPF problem can be implemented by integrating sample points in the search space, we believe that, for small test
QGS (12) to a point lying on one regular stable equilibrium systems, all the SEMs in the search space may be located as
manifold. This point is a feasible solution, as assured by long as the sample size is sufficiently large, although this claim
Theorem 5. cannot be guaranteed.
3) The global convergence property of this trajectory-based We consider a 9-bus system composed of three generators
method is assured by Theorem 6. However, some trajec- and nine branches [13]. By Fact 1, the feasible region of the
tories may converge to degenerate SEMs. Step 3 partially 9-bus system is 5-dimensional. The feasible region and stable
deals with this situation. An effective method to system- equilibrium manifolds on the PG projection plane are shown in
atically escape from a degenerate SEM and move toward Fig. 5, where the feasible region is obtained in the same way
a regular SEM is needed. as that obtained in Fig. 4(a). The numerical results shown in
Some distinguishing features of the QGS-based method in- Fig. 5 suggest that this 9-bus test system has three disconnected
clude the following: feasible components containing four OPF solutions with the
i) The convergence region of a feasible component of local optima and global optimum marked with ‘×’ and ‘◦’,
an OPF problem under the QGS-based method is the respectively. We have the following observations:
CHIANG AND JIANG: FEASIBLE REGION OF OPTIMAL POWER FLOW: CHARACTERIZATION AND APPLICATIONS 241

Fig. 6. Local feasible components and local optimal solutions of the 118-bus
Fig. 5. The feasible region and stable equilibrium manifolds of the 9-bus sys-
system. (a) Local feasible components containing the three local OPF solutions
tem. The feasible region is composed of three feasible components: the local
on the P G 6 9 × P G 7 0 × P G 7 2 projection plane. (b) Local feasible compo-
OPF solutions are located in two feasible components, and the global OPF solu-
nents containing the three local OPF solutions on the P G 6 9 × P G 7 0 projection
tion is located in the other feasible component. (a) Feasible components, three
plane.
local OPF solutions, and one global OPF solution on the P G 1 × P G 2 × P G 3
projection plane. (b) Stable equilibrium manifolds on the P G 1 × P G 2 × P G 3
projection plane. There are three regular stable equilibrium manifolds and one
degenerate stable equilibrium manifold in the state space of the corresponding feasible region and then projecting the feasible region onto a 2-
QGS. or 3-dimensional subspace) is employed to visualize the feasi-
ble region. However, the projection method is not suitable for
1) Each feasible component corresponds to a SEM. visualizing the feasible region of large OPF problems because
2) There is one degenerate SEM (the red point) whose en- of the two main problems:
ergy function value is non-zero, and thus, not a feasible 1) The search space is large
component. 2) Even if the original feasible components are disconnected,
3) The energy function achieves its local minimum at the the projection may be connected.
degenerate SEM and is consistent with Theorem 4. It is indeed very challenging to completely visualize the en-
While Fig. 5 suggests that the feasible region of this 9-bus tire high-dimensional feasible region of an OPF problem. To
system is composed of three disconnected components, there is this end, the notion of a local feasible region is proposed. We
no rigorous proof for this suggestion. In one respect, it is possible argue that, for high-dimensional OPF problems, one is usually
(though unlikely) that the feasible region in Fig. 5 is connected interested in local feasible regions surrounding some important
via very thin regions that are not detected. On the other hand, points such as the current operating point or certain local optimal
the number of disconnected components in the projection can solutions.
only give a lower bound for the number in the original space. A 118-bus OPF test system from [13] is examined. We adopt
The above numerical studies are in agreement with the three-point strategy used in [20] to compute the local feasible
Theorem 2, Theorem 4, and Theorem 5, and they confirm that region:
the union of regular stable equilibrium manifolds of the QGS Step 1: Choose 3 points (a, b, c) of interest in the search
(12) indeed equals the feasible region of OPF problems. space and define the following 2-dimensional plane in high-
dimensional search space:

B. The Local Feasible Region y = a + (b − a) s + (c − a) t,

In the above small test systems, a simple projection method
(i.e., random initial points in the search space for finding the where a, b, c and y ∈ n , and s and t are scalars.
242 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 33, NO. 1, JANUARY 2018

Step 2: Starting from some points lying on the 2-dimensional the fact that certifying the feasibility of an OPF problem has been
plane defined in Step 1, integrate the corresponding trajectories proven to be strongly NP-hard [22]. In the proposed QGS-based
until they converge to some SEMs of the QGS. method, the main challenge is that even if (7) is feasible, the
Step 3: Select 2 or 3 variables from the control variables of trajectories may converge to degenerate SEMs. Therefore, an
the OPF problem. effective method to escape from degenerate SEMs and move
Step 4: Obtain a local feasible region by projecting these toward regular SEMs is needed.
points on the SEMs, obtained in Step 2, onto the projection The analytical results developed in this paper can be fur-
plane defined by the variables selected in Step 3. ther explored. Additional development of the visualization tech-
Remarks: niques for displaying the local feasible regions is needed. When
1) Step 1 defines a 2-dimensional cutting plane that contains the QGS-based method converges to a degenerate SEM, the
the 3 points of interest, while Step 2 obtains a local feasible problem of how to recover a feasible solution by changing cer-
region near the cutting plane. tain constraint parameters via exploring the degenerate SEM
2) Step 3 and Step 4 project the obtained local feasible region requires further investigation.
onto a 2- or 3-dimensional subspace for better visualiza-
tion. APPENDIX
For the 118-bus test system, the three known local optimal
A. Proof of Theorem 2: Proof: First, we will prove that
solutions (marked with ‘×’ in Fig. 6) are chosen as the points
each feasible region component F Ri of constraint set (7) is a
of interest and the PG 69 × PG 70 × PG 72 plane is chosen as
pseudo-hyperbolic equilibrium manifold of QGS (12).
the projection plane because, on this PG subspace, the sum of
Let x̄ be a feasible point on F Ri , i.e., x̄ ∈ F Ri ,H(x̄) = 0.
the Euclidean distances between pairs of the three local opti-
Then the Jacobian of system (12) at x̄ is:
mal solutions is the maximum. The local feasible components
containing the three local optimal solutions are shown in Fig. 6. DQH (x̄) = −D(DH(x̄)T )H(x̄) − DH(x̄)T DH(x̄)
We observe from Fig. 6 that the three known local optimal solu-
= −DH(x̄)T DH(x̄)
tions are located on two disconnected feasible components, but
this observation needs further examination. We believe that this DQH (x̄) is a real symmetric matrix, so the eigenvalues of
notion of local feasible regions may offer valuable insights and DQH (x̄) that have zero real parts can only be zeros.
information regarding the entire feasible region. Hopefully, this Assume that rank(DH(x̄)) = r. Then it is obvious that
work will spur more research into improving this notion.
dim(DH(x̄)⊥ ) = n − r
where DH(x̄)⊥ := {u : DH(x̄)u = 0}.
VII. CONCLUSION
Let the basis vectors that form DH(x̄)⊥ be
We have developed a mathematical characterization of the
feasible region of nonlinear optimal power flow problems. This vi , i = 1, 2, . . . , n − r.
characterization is developed through a QGS derived from the Then, for each vi ,
set of equality and inequality constraints of the underlying non-
linear optimal power flow problem. An equivalence between DH(x̄)vi = 0 ⇒ −DH(x̄)T DH(x̄)vi = 0, i = 1, 2, . . . , n−r.
the entire feasible region of an OPF problem and the union of Thus, vi is an eigenvector associated with the zero eigen-
regular stable equilibrium manifolds of the QGS has been es- value. Since rank(DH(x̄)T DH(x̄)) = rank(DH(x̄)) = r,
tablished. It is further shown that the QGS is completely stable the number of zero eigenvalues of DQH (x̄) is n − r. Hence, ev-
and that each trajectory converges to an equilibrium manifold, ery eigenvector associated with the zero eigenvalue of DQH (x̄)
making the QGS-based method useful in locating OPF feasible can be represented by vi , i = 1, 2, . . . , n − r.
solutions with the following features: Moreover, for any vector u ∈ DH(x̄)⊥ , u is on the tangent
1) It is not sensitive to initial conditions in the sense that space Tx̄ (F Ri ) [21], that is, all the eigenvectors associated
the initial points converge to different points lying on the with the zero eigenvalues are on the tangent space Tx̄ (F Ri ).
same SEM as long as the initial points lie inside the same Therefore, according to Definition 2 and Definition 3, F Ri is a
stability region. pseudo-hyperbolic equilibrium manifold of QGS (12).
2) It does not require matrix inversion, hence possible ill- Next, we show that F Ri is a pseudo-hyperbolic equilibrium
conditioning of the Jacobian matrix along a trajectory is manifold and a stable equilibrium manifold of system (12):
avoided. 1) For any vector u ∈ / DH(x̄)⊥ ,
3) The QGS-based method provides a novel and promising
trajectory-based way to compute feasible solutions. This uT DQH (x̄)u = −(DH(x̄)u)T DH(x̄)u
method can be combined with other methods, such as the = − DH(x̄)u 2
< 0.
interior point method or other optimization-based meth-
ods, to obtain feasible solutions as well as local optimal 2) For any vector u ∈ DH(x̄)⊥ , we have uT DQH (x̄)u = 0.
solutions. So, the Jacobian of system (12) at x̄ is negative semi-definite.
In future work, we plan to develop effective numerical QGS- Moreover, F Ri is a pseudo-hyperbolic equilibrium manifold.
based methods for fast computation of feasible solutions, despite Therefore, according to Definition 3 and Definition 4, F Ri is a
CHIANG AND JIANG: FEASIBLE REGION OF OPTIMAL POWER FLOW: CHARACTERIZATION AND APPLICATIONS 243

stable equilibrium manifold of system (12). This completes the is constant and SVV2 h is nonnegative. Similarly, (H7’) becomes
proof. unbounded (positive infinity) if V goes to negative infinity. In
B. Proof of Theorem 3: Proof: We will show that function other words, E(x) always becomes unbounded, regardless of
E(x) satisfies the three conditions of energy functions (E1), whether V goes to positive or negative infinity.
(E2), and (E3). Case 2.2: If V is bounded but one or more variables in
First, voltage angle constraints are added to constraint set (7): {SVV h , SVV l } go to infinity, then (H6’) becomes unbounded as
SVV h goes to infinity, since V is bounded and V m ax is constant.
θi 2 ≤ (θim ax )2 i ∈ {1, . . . , NB }
Similarly, (H7’) becomes unbounded as SVV l goes to infinity.
where θim ax is usually chosen as 2π, which does not affect the Therefore, E(x) becomes unbounded if either SVV h or SVV l
final solutions of constraint set (7) and will help avoid periodic goes to infinity.
solutions. To summarize, the value of E(x) becomes unbounded if one
By adding slack variables to each inequality constraint, the or more variables in {V, SVV h , SVV l } go to infinity. By anal-
constraint functions of OPF problems can be defined: ogy, the cases for one or more variables in {PG , SVP h , SVP l }
⎛ ⎞ ⎛  ⎞ or {QG , SVQ h , SVQ l } that go to infinity can be similarly ana-
PG − PL − P (θ, V ) H1
⎜ ⎟ ⎜  ⎟ lyzed. The above discussion leads to the conclusion that E(x)
⎜ QG − QL − Q(θ, V ) ⎟ ⎜H2 ⎟ becomes unbounded as well if x → ∞. Hence, condition (E3)
⎜ ⎟ ⎜ ⎟
⎜ θ2 − (θm ax )2 + SV 2 ⎟ ⎜H3 ⎟ is verified.
⎜ θ ⎟ ⎜ ⎟
⎜S 2 − (S m ax )2 + SV 2 ⎟ ⎜  ⎟ Above all, E(x) satisfies conditions (E1) – (E3) and is an
⎜ f l Sf ⎟ ⎜H4 ⎟
⎜ 2 ⎟ ⎜ ⎟ energy function for system (12). This completes the proof.
⎜ St − (Slm ax )2 + SVS2t ⎟ ⎜H5 ⎟
⎜ ⎟ ⎜ ⎟ C. Proof of Theorem 4: Proof: Let Σs be a stable equi-
⎜ ⎟ ⎜ ⎟
H(x) = ⎜ V − V m ax + SVV2 h ⎟ = ⎜H6 ⎟ = 0. librium manifold of (12). We first show that E(x) = constant
⎜ ⎟ ⎜ ⎟
⎜ V m in − V + SVV2 l ⎟ ⎜H7 ⎟ for all x ∈ Σs . Fix x0 ∈ Σs and let x1 ∈ Σs too. Then, by
⎜ ⎟ ⎜ ⎟
⎜ P − P m ax + SV 2 ⎟ ⎜  ⎟ the connectedness of equilibrium manifold Σs , there exists
⎜ G G Ph ⎟ ⎜H8 ⎟
⎜ m in ⎟ ⎜ ⎟ a path c :  → Σs such that c(0) = x0 and c(1) = x1 . Since
⎜ PG − PG + SVP2 l ⎟ ⎜H9 ⎟
⎜ ⎟ ⎜ ⎟ ∇E(c(s)) = DH(c(s))T H(c(s)) = 0 for all s, we have
⎜ Q − Qm ax + SV 2 ⎟ ⎜ ⎟
⎝ G G Qh ⎠ ⎝ H10 ⎠
QG − QG + SVQ l
m in 2
H11  d d
E(c(s)) = ∇E(c(s))T c(s) = 0.
ds ds
By differentiating E(x), one obtains:
Hence, E(c(s)) = constant for all s, which implies that
dE(x)  
= ∇E(x), ẋ = DH(x)T · H(x), −DH(x)T · H(x) E(x0 ) = E(c(0)) = E(c(1)) = E(x1 ). Therefore, E(x) = a
dt constant for all x ∈ Σs . Next, since Σs is stable, there exists
 2 an ε > 0 such that Φ(t, z) → Σs for each z ∈ Bε (Σs ). Since
 
= −DH(x)T · H(x) ≤ 0 E(Φ(t, z)) is an energy function of t for each z ∈ Bε (Σs ), we
have E(z) ≥ E(Σs ), which implies that Σs is locally optimal
Therefore, condition (E1) holds. Clearly, dEdt(x) = 0 if and for the following minimization problem:
only if x is an equilibrium point. Thus, condition (E2) is also
satisfied. It remains to be shown that condition (E3) is valid. To min E(x).
this end, we will prove the inverse negative proposition of (E3),
i.e., If x → ∞ then E(x) = 12 H(x) 2 → ∞. Moreover, if the global minimal value (i.e., E(x) = 0) is
Observe that if x → ∞, then at least one variable in set achieved, then Σs is a feasible component of constraint set (7).
{θ, V, PG , QG , SV } goes to infinity. In this regard, we discuss This completes the proof.
the following cases. D. Proof of Theorem 5: Proof: We will show the necessity
Case 1: If one or more variables in {θ, SVθ , SVS f , SVS t } go and sufficiency. The sufficiency follows from Theorem 4 and the
to infinity, then one can easily check that either (H3’), (H4’), definition of non-degenerate stable equilibrium. Assume Σr s
or (H5’) must become unbounded because θm ax and Slm ax are is a regular stable equilibrium manifold of QGS (12). Then
fixed constants, while the terms θ2 , SVθ2 , SVS2f and SVS2t are H(Σr s ) = 0, i.e., Σr s is a feasible component of (7).
all nonnegative. This shows that E(x) becomes unbounded if Next, we show the necessity. Assume that F Rj is a feasi-
one or more variables in {θ, SVθ , SVS f , SVS t } go to infinity. ble component of constraint set (7). According to Theorem 2,
Case 2: If one or more variables in the complement F Rj is a stable equilibrium manifold of (12). Moreover, if
{V, SVV h , SVV l , PG , SVP h , SVP l , QG , SVQ h , SVQ l } go to in- H(F Rj ) = 0, then F Rj is a non-degenerate stable equilibrium
finity, the same conclusion can be derived as follows. To fix the manifold. This completes the proof.
ideas, one can start by considering the case in which one or more E. Proof of Theorem 6: Proof: The conditions (E1) and
variables in {V, SVV h , SVV l } go to infinity. Indeed, if one or (E2) of the energy functions ensure that every trajectory of
more variables in {V, SVV h , SVV l } go to infinity, the discussion system (12) either converges to an equilibrium manifold or
can be further sub-divided below for the sake of clarity. diverges. We next show that it is not possible for a trajectory of
Case 2.1: If V goes to positive infinity, one observes that the QGS to diverge. This is because of the following fact, which
(H6’) must become unbounded (positive infinity), since V m ax can be found in the proof of Theorem 3:
244 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 33, NO. 1, JANUARY 2018

If x → ∞, then the energy function of (12),E(x) = [18] J. Lee and H. D. Chiang, “Theory of stability regions for a class of nonhy-
1
H(x) 2 , goes to infinity. perbolic dynamical systems and its application to constraint satisfaction
2 problems,” IEEE Trans. Circuits Syst. I, vol. 49, no. 2, pp. 196–209,
Since the energy function decreases along any non-trivial tra- Feb. 2002.
jectory, it follows that every trajectory of system (12) converges [19] W. A. Bukhsh, A. Grothey, K. McKinnon, and P. A. Trodden. Test Case
to an equilibrium manifold. This completes the proof. Archive of Optimal Power Flow (OPF) Problems With Local Optima.
[Online]. Available: http://www.maths.ed.ac.uk/optenergy/LocalOpt/
[20] W. Suampun and H. D. Chiang, “Critical evaluation of methods for es-
ACKNOWLEDGMENT timating stability boundary for transient stability analysis in power sys-
tems,” in Proc. IEEE Power Energy Soc. General Meeting, 2010, pp. 1–8.
The authors would like to thank the reviewers for their de- [21] J. Milnor, Topology From the Differentiable Viewpoint. Charlottesville,
tailed and insightful comments, which greatly improved the VA, USA: Univ. Virginia Press, 1965.
[22] D. Bienstock and A. Verma, “Strong NP-hardness of AC power flows
presentation of this paper. feasibility,” arXiv preprint arXiv: 1512.07315, 2015.
[23] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MAT-
REFERENCES POWER: Steady-state operations, planning and analysis tools for power
systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1,
[1] J. Carpentier, “Contribution to the economic dispatch problem,” Bull. de pp. 12–19, Feb. 2011.
la Société Française des Electriciens, vol. 3, no. 8, pp. 431–447, 1962. [24] K. C. Almeida and A. Kocholik, “Solving Ill-posed optimal power flow
[2] M. Huneault and F. D. Galiana, “A survey of the optimal power flow problems via Fritz-John optimality conditions,” IEEE Trans. Power Syst.,
literature,” IEEE Trans. Power Syst., vol. 6, no. 2, pp. 762–770, May 1991. vol. 31, no. 6, pp. 4913–4922, Nov. 2016.
[3] J. A. Momoh, M. E. El-Hawary, and R. Adapa, “A review of selected
optimal power flow literature to 1993. Part I & II,” IEEE Trans. Power
Syst., vol. 14, no. 1, pp. 96–111, Feb. 1999.
[4] B. C. Lesieutre and I. A. Hiskens, “Convexity of the set of feasible injec-
tions and revenue adequacy in FTR markets,” IEEE Trans. Power Syst.,
vol. 20, no. 4, pp. 1790–1798, Nov. 2005.
[5] K. S. Pandya and S. K. Joshi, “A survey of optimal power flow methods,”
J. Theor. Appl. Inf. Technol., vol. 4, no. 5, pp. 450–458, 2008. Hsiao-Dong Chiang (M’87–SM’91–F’97) received
[6] J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow problem,” the Ph.D. degree in electrical engineering and com-
IEEE Trans. Power Syst., vol. 27, no. 1, pp. 92–107, Feb. 2012. puter sciences from the University of California at
[7] D. K. Molzahn, J. T. Holzer, B. C. Lesieutre, and C. L. DeMarco, “Imple- Berkeley, Berkeley, CA, USA. Since 1998, he has
mentation of a large-scale optimal power flow solver based on semidefinite been a Professor with the School of Electrical and
programming,” IEEE Trans. Power Syst., vol. 28, no. 4, pp. 3987–3998, Computer Engineering, Cornell University, Ithaca,
Nov. 2013. NY, USA. His current research interests include non-
[8] S. Low, “Convex relaxation of optimal power flow: Parts I & II,” IEEE linear system theory, nonlinear computation, nonlin-
Trans. Control Netw. Syst., vol. 1, no. 1, pp. 15–27, Mar. 2014. ear optimization, and their practical applications. He
[9] D. Molzahn and I. Hiskens, “Sparsity-exploiting moment-based relax- holds 17 U.S. and overseas patents and several con-
ations of the optimal power flow problem,” IEEE Trans. Power Syst., sultant positions. He is the author of two books: Direct
vol. 30, no. 6, pp. 3168–3180, Nov. 2015. Methods for Stability Analysis of Electric Power Systems: Theoretical Founda-
[10] J. Lavaei, D. Tse, and B. Zhang, “Geometry of power flows and optimiza- tions, BCU Methodologies, and Applications (Hoboken, NJ, USA: Wiley, 2011)
tion in distribution networks,” IEEE Trans. Power Syst., vol. 29, no. 2, and (with Luis F. C. Alberto) Stability Regions of Nonlinear Dynamical Sys-
pp. 572–583, Mar. 2014. tems: Theory, Estimation, and Applications (Cambridge, U.K.: Cambridge Univ.
[11] A. Schecter and R. P. O’Neill. (2013). Exploration of the ACOPF Feasible Press, 2015). He has served as an Associate Editor for several IEEE transactions
Region for the Standard IEEE Test Set. FERC Staff Technical Paper. [On- and journals and is the founder of Bigwood Systems, Inc., Ithaca.
line]. Available: http://www.ferc.gov/industries/electric/indus-act/market-
planning/opf-papers/acopf-6-test-problem-properties.pdf
[12] B. Zhang and D. Tse, “Geometry of injection regions of power networks,”
IEEE Trans. Power Syst., vol. 28, no. 2, pp. 788–797, May 2013.
[13] W. A. Bukhsh, A. Grothey, K. McKinnon, and P. A. Trodden, “Local
solutions of the optimal power flow problem,” IEEE Trans. Power Syst.,
vol. 28, no. 4, pp. 4780–4788, Nov. 2013.
[14] D. Alves and G. da Costa, “An analytical solution to the optimal power
flow,” IEEE Power Eng. Rev., vol. 22, no. 3, pp. 49–51, Mar. 2002. Chu-Yang Jiang received the B.Sc. degree in elec-
[15] H. D. Chiang and L. F. C. Alberto, Stability Regions of Nonlinear Dynam- trical engineering and automation from Tianjin Uni-
ical Systems: Theory, Estimation, and Applications. Cambridge, U.K.: versity, Tianjin, China, in 2011, where he is currently
Cambridge Univ. Press, 2015. working toward the Ph.D. degree at the School of
[16] H. D. Chiang, M. W. Hirsch, and F. F. Wu, “Stability region of nonlinear Electrical Engineering and Automation. His areas of
autonomous dynamical systems,” IEEE Trans. Autom. Control, vol. 33, interest include optimal power flow problems, nonlin-
no. 1, pp. 16–27, Jan. 1988. ear dynamic systems, and their application to electric
[17] H. D. Chiang, Direct Methods for Stability Analysis of Electric Power power systems.
Systems: Theoretical Foundations, BCU Methodologies, and Applications.
Hoboken, NJ, USA: Wiley, 2011.