Strong SOCP Relaxations for the Optimal Power Flow Problem

Strong SOCP Relaxations for the Optimal Power Flow Problem

Burak Kocuk, Santanu S. Dey, X. Andy Sun
Abstract

This paper proposes three strong second order cone programming (SOCP) relaxations for the AC optimal power flow (OPF) problem. These three relaxations are incomparable to each other and two of them are incomparable to the standard SDP relaxation of OPF. Extensive computational experiments show that these relaxations have numerous advantages over existing convex relaxations in the literature: (i) their solution quality is extremely close to that of the SDP relaxations (the best one is within of the SDP relaxation on average for all the IEEE test cases) and consistently outperforms previously proposed convex quadratic relaxations of the OPF problem, (ii) the solutions from the strong SOCP relaxations can be directly used as a warm start in a local solver such as IPOPT to obtain a high quality feasible OPF solution, and (iii) in terms of computation times, the strong SOCP relaxations can be solved an order of magnitude faster than standard SDP relaxations. For example, one of the proposed SOCP relaxations together with IPOPT produces a feasible solution for the largest instance in the IEEE test cases (the 3375-bus system) and also certifies that this solution is within of global optimality, all this computed within seconds on a modest personal computer. Overall, the proposed strong SOCP relaxations provide a practical approach to obtain feasible OPF solutions with extremely good quality within a time framework that is compatible with the real-time operation in the current industry practice.

1 Introduction

Optimal Power Flow (OPF) is a fundamental optimization problem in electrical power systems analysis (Carpentier 1962). There are two challenges in the solution of OPF. First, it is an operational level problem solved every few minutes, hence the computational budget is limited. Second, it is a nonconvex optimization problem on a large-scale power network of thousands of buses, generators, and loads. The importance of the problem and the aforementioned difficulties have produced a rich literature, see e.g. (Momoh et al. 1999a, b, Frank et al. 2012a, b, Cain et al. 2012).

Due to these challenges, the current practice in the electricity industry is to use the so-called DC OPF approximation (FERC 2011). In contrast, the original nonconvex OPF is usually called the AC (alternating current) OPF. DC OPF is a linearization of AC OPF by exploiting some physical properties of the power flows in typical power systems, such as tight bounds on voltage magnitudes at buses and small voltage angle differences between buses. However, such an approximation completely ignores important aspects of power flow physics, such as the reactive power and voltage magnitude. To partially remedy this drawback, the current practice is to solve DC OPF and then to solve a set of power flow equations with the DC OPF solution to compute feasible reactive powers and voltages. However, it is clear such an approach cannot guarantee any optimality of the AC power flow solution obtained. To be concise, we will use OPF to denote AC OPF in the remainder of the paper.

In order to solve the OPF problem, the academic literature has focused on improving nonlinear optimization methods such as the interior point methods (IPM) to compute local optimal solutions, see e.g. (Wu et al. 1994, Torres and Quintana 1998, Jabr et al. 2002, Wang et al. 2007). A well-known implementation of IPM tailored for the OPF problem is MATPOWER (Zimmerman et al. 2011). Although these local methods are effective in solving IEEE test instances, they do not offer any quantification of the quality of the solution.

In the recent years, much research interests have been drawn to the convex relaxation approach. In particular, the second-order cone programming (SOCP) and the semidefinite programming (SDP) relaxations are first applied to the OPF problem in Jabr (2006), and Bai et al. (2008), Bai and Wei (2009) and Lavaei and Low (2012). Among these two approaches, the SDP relaxation and its variations have drawn significant attention due to their strength. Since convex conic programs are polynomially solvable, the SDP relaxation offers an effective way for obtaining global optimal solutions to OPF problems whenever the relaxation is exact. Unfortunately, the exactness of the SDP relaxations can be guaranteed only for a restricted class of problems under some assumptions, e.g. radial networks (Zhang and Tse 2011, Bose et al. 2011, 2012) under load over-satisfaction (Sojoudi and Lavaei 2012) or absence of generation lower bounds (Lavaei et al. 2014), or lossless networks with cyclic graphs (Zhang and Tse 2013). A comprehensive survey can be found in Low (2014a, b). When the SDP relaxation is not exact, it may be difficult to put a physical meaning on the solution.

A way to further strengthen the SDP relaxation is to solve a hierarchy of moment relaxation problems (Lasserre 2001, Parrilo 2003). This approach is used in Josz et al. (2015) to globally solve small-size problems, and is also used in Molzahn and Hiskens (2015) to obtain tighter lower bounds for larger problems of 300-bus systems. However due to the NP-hardness of the OPF problem (Lavaei and Low 2012), in general the order of the Lasserre hierarchy required to obtain a global optimal solution can be arbitrarily large. Furthermore, even the global optimal objective function value is achieved, the solution matrices may not be rank one, which poses another challenge in terms of recovering an optimal voltage solution (Lavaei et al. 2014). This indicates the computational difficulty of the SDP relaxation approach to practically solve real-world sized power networks with more than a thousand buses. For such large-scale OPF problems, a straightforward use of IPM to solve the SDP relaxation becomes prohibitively expensive. Interesting works have been done to exploit the sparsity of power networks as in Jabr (2012), Molzahn et al. (2013), Madani et al. (2014b), Molzahn and Hiskens (2015), Madani et al. (2015). The underlying methodology utilizes techniques such as chordal graph extension, tree-width decomposition, and matrix completion, as proposed and developed in Fukuda et al. (2001) and Nakata et al. (2003).

More recently, there is a growing trend to use computationally less demanding relaxations based on linear programming (LP) and SOCP to solve the OPF problem. For instance, linear and quadratic envelopes for trigonometric functions in the polar formulation of the OPF problem are constructed in Coffrin and Van Hentenryck (2014), Hijazi et al. (2013), Coffrin et al. (2015). In Bienstock and Munoz (2014), LP based outer approximations are proposed which are strengthened by incorporating several different types of valid inequalities.

This paper proposes new strong SOCP relaxations of the OPF problem and demonstrates their computational advantages over the SDP relaxations and previously described convex quadratic relaxations for the purpose of practically solving large-scale OPF problems. Our starting point is an alternative formulation for the OPF problem proposed in Expósito and Ramos (1999) and Jabr (2006). In this formulation, the nonconvexities are present in two types of constraints: one type is the surface of a rotated second-order cone, and the other type involves arctangent functions on voltage angles. The SOCP relaxation in Jabr (2006) is obtained by convexifying the first type of constraints to obtain SOCP constraints and completely ignoring the second type constraints. We refer to this relaxation as the classic SOCP relaxation of the OPF problem. We prove that the standard SOCP relaxation of the rectangular formulation of OPF provides the same bounds as the classic SOCP relaxation. Therefore, if we are able to add convex constraints that are implied by the original constraints involving the arctangent function to the classic SOCP relaxation, then this could yield stronger relaxation than the classic SOCP relaxation that may also potentially be incomparable to (i.e. not dominated by nor dominates) the standard SDP relaxations. In this paper, we propose three efficient ways to achieve this goal.

In the following, we summarize the key contributions of the paper.

  1. We theoretically analyze the relative strength of the McCormick (linear programming), SOCP, and SDP relaxations of the rectangular and alternative formulations of the OPF problem. As discussed above, this analysis leads us to consider strengthening the classic SOCP relaxation as a way forward to obtaining strong and tractable convex relaxations.

  2. We propose three efficient methods to strengthen the classic SOCP relaxation.

    1. In the first approach, we begin by reformulating the arctangent constraints as polynomial constraints whose degrees are proportional to the length of the cycles. This yields a bilinear relaxation of the OPF problem in extended space (that is by addition of artificial variables), where the new variables correspond to edges obtained by triangulating cycles. With this reformulation, we use the McCormick relaxation of the proposed bilinear constraints to strengthen the classic SOCP relaxation. The resulting SOCP relaxation is shown to be incomparable to the SDP relaxation.

    2. In the second approach, we construct a polyhedral envelope for the arctangent functions in 3-dimension, which are then incorporated into the classic SOCP relaxation. This SOCP relaxation is also shown to be incomparable to the standard SDP relaxation.

    3. In the third approach, we strengthen the classic SOCP relaxation by dynamically generating valid linear inequalities that separate the SOCP solution from the SDP cone constraints over cycles. We observe that running such a separation oracle a few iterations already produces SOCP relaxation solutions very close to the quality of the full SDP relaxation.

  3. We conduct extensive computational tests on the proposed SOCP relaxations and compare them with existing SDP relaxations (Lavaei et al. 2014) and quadratic relaxations (Coffrin and Van Hentenryck 2014, Coffrin et al. 2015). The computational results can be summarized as follows.

    1. Lower bounds: The lower bounds obtained by the third proposed SOCP relaxation for all MATPOWER test cases from 6-bus to 3375-bus are on average within of the lower bounds of the SDP relaxation. The other two proposed relaxation are also on average within of the SDP relaxation.

    2. Computation time: Overall, the proposed SOCP relaxations can be solved orders of magnitude faster than the SDP relaxations. The computational advantage is even more evident when a feasible solution of the OPF problem is needed. As an example, consider the largest test instance of the IEEE 3375-bus system. Our proposed SOCP relaxation together with IPOPT provides a solution for this instance and also certifies that this solution is within of global optimality, all computed in seconds on a modest personal computer.

    3. Comparison with other convex quadratic relaxation: The proposed SOCP relaxations consistently outperform the existing quadratic relaxation in Coffrin and Van Hentenryck (2014) and Coffrin et al. (2015) on the test instances of typical, congested, and small angle difference conditions.

    4. Non-dominance with standard SDP relaxation: The computation also shows that the proposed SOCP relaxations are neither dominated by nor dominates the standard SDP relaxations.

    5. Robustness: The proposed SOCP relaxations perform consistently well on IEEE test cases with randomly perturbed load profiles.

The paper is organized as follows. Section 2 introduces the standard rectangular formulation and the alternative formulation of the OPF problem. Section 3 compares six different convex relaxations for the OPF problem based on the rectangular and alternative formulations. Section 4 proposes three ways to strengthen the classic SOCP relaxation. Section 5 presents extensive computational experiments. We make concluding remarks in Section 6.

2 Optimal Power Flow Problem

Consider a power network , where denotes the node set, i.e., the set of buses, and denotes the edge set, i.e., the set of transmission lines. Generation units (i.e. electric power generators) are connected to a subset of buses, denoted as . We assume that there is electric demand, also called load, at every bus. The aim of the optimal power flow problem is to satisfy demand at all buses with the minimum total production costs of generators such that the solution obeys the physical laws (e.g., Ohm’s Law and Kirchoff’s Law) and other operational restrictions (e.g., transmission line flow limit constraints).

Let denote the nodal admittance matrix, which has components for each line , and , where (resp. ) is the shunt conductance (resp. susceptance) at bus and . Let (resp. ) be the real and reactive power output of the generator (resp. load) at bus . The complex voltage (also called voltage phasor) at bus can be expressed either in the rectangular form as or in the polar form as , where is the voltage magnitude and is the angle of the complex voltage. In power system analysis, the voltage magnitude is usually normalized against a unit voltage level and is expressed in per unit (p.u.). For example, if the unit voltage is 100kV, then 110kV is expressed as 1.1 p.u.. In transmission systems, the bus voltage magnitudes are usually restricted to be close to the unit voltage level to maintain system stability.

With the above notation, the OPF problem is given in the so-called rectangular formulation:

(1a)
(1b)
(1c)
(1d)
(1e)
(1f)

Here, the objective function is typically linear or convex quadratic in the real power output of generator . Constraints (1b) and (1c) correspond to the conservation of active and reactive power flows at each bus, respectively. denotes the set of neighbor buses of bus . Constraint (1d) restricts voltage magnitude at each bus. As noted above, and are both close to 1 p.u. at each bus . Constraints (1e) and (1f), respectively, limit the active and reactive power output of each generator to respect its physical capability.

Note that the rectangular formulation (1) is a nonconvex quadratic optimization problem. However, quite importantly, we can observe that all the nonlinearity and nonconvexity comes from one of the following three forms: (1) , (2) , (3) . To capture this nonlinearity, we define new variables , and for each bus and each transmission line as , , . These new variables satisfy the relation . With a change of variables, we can introduce an alternative formulation of the OPF problem as follows:

(2a)
(2b)
(2c)
(2d)
(2e)
(2f)

This formulation was first introduced in Expósito and Ramos (1999) and Jabr (2006). It is an exact formulation for OPF on a radial network in the sense that the optimal power output of (2) is also optimal for (1) and one can always recover the voltage phase angles ’s by solving the following system of linear equations with the optimal solution :

(3)

which then provide an optimal voltage phasor solution to (1) (see e.g. Zhang and Tse (2011)). Here, in order to cover the entire range of , we use the function111, which takes value in , rather than as is the case of the regular arctangent function. Unfortunately, for meshed networks, the above formulation (2) can be a strict relaxation of the OPF problem. The reason is that, given an optimal solution for all edges of (2), it does not guarantee that sums to zero over all cycles. In other words, the optimal solution of (2) may not be feasible for the original OPF problem (1). This issue can be fixed by directly imposing (3) as a constraint (Jabr 2007, 2008). Thus (2) together with (3) is a valid formulation for OPF in mesh networks. Note that the constraints involving the function are nonconvex.

Line Flow Constraints: Typically, the OPF problem also involves the so-called transmission constraints on transmission lines. In the literature, different types of line flow limits are suggested. We list a few of them below (Madani et al. 2015, Coffrin et al. 2015) and present how they can be expressed in the space of variables:

  1. Real power flow on line :

  2. Squared voltage difference magnitude on line :

  3. Squared current magnitude on line :

  4. Squared apparent power flow on line :

  5. Angle difference on line : or , see equation (55) for details

We omit such constraints for the brevity of discussion. However, last two constraints are included in our computational experiments whenever explicit bounds are given.

3 Comparison of Convex Relaxations

In this section, we first present six different convex relaxations of the OPF problem. In particular, we consider the McCormick, SOCP, and SDP relaxations of both the rectangular formulation (1) and the alternative formulation (2). Then, we analyze their relative strength by comparing their feasible regions. This comparison is an important motivator for the approach we take in the rest of the paper to generate strong SOCP relaxations. We discuss this in Section 3.3.

3.1 Standard Convex Relaxations

3.1.1 McCormick Relaxation of Rectangular Formulation ().

As shown in McCormick (1976), the convex hull of the set is given by

which we denote as . We use this result to construct McCormick envelopes for the quadratic terms in the rectangular formulation (1). In particular, let us first define the following new variables for each edge : , , , and for each bus : , . Consider the following set of constraints:

(4a)
(4b)
(4c)
(4d)
(4e)
(4f)

where the McCormick envelops in (4e)-(4f) are constructed using the bounds given in (4d). Using these McCormick envelopes, we obtain a convex relaxation of the rectangular formulation (1) with the feasible region denoted as

(5)

Note that this feasible region is a polyhedron. If the objective function is linear, then we have a linear programming relaxation of the OPF problem.

3.1.2 McCormick Relaxation of Alternative Formulation ().

Using the similar technique on the alternative formulation (2), we define new variables , , for each edge , and consider the following set of constraints:

(6a)
(6b)
(6c)

where the McCormick envelops in (6c) are constructed using the bounds given in (6b) and (2d). Denote the feasible region of the corresponding convex relaxation as

(7)

Again, is a polyhedron.

3.1.3 SDP Relaxations of Rectangular Formulation ().

To apply SDP relaxation to the rectangular formulation (1), define a hermitian matrix , i.e., , where is the conjugate transpose of . Consider the following set of constraints:

(8a)
(8b)
(8c)
(8d)

where and are the real and imaginary parts of the complex number , respectively. Let denote the vector of voltage phasors with the -th entry for each bus . If , then rank and rank are both equal to 2, and (8a)-(8c) exactly recovers (1b)-(1d). By ignoring the rank constraints, we come to the standard SDP relaxation of the rectangular formulation (1), whose feasible region is defined as

(9)

This SDP relaxation is in the complex domain. There is also an SDP relaxation in the real domain by defining . The associated constraints are

(10a)
(10b)
(10c)
(10d)

where and . If , i.e. rank()=1, then (10a)-(10c) exactly recovers (1b)-(1d). We denote the feasible region of this SDP relaxation in the real domain as

(11)

The SDP relaxation in the real domain is first proposed in Bai et al. (2008), Bai and Wei (2009), and Lavaei and Low (2012). The SDP relaxation in the complex domain is formulated in Bose et al. (2011) and Zhang and Tse (2013) and is widely used in the literature now for its notational simplicity. These two SDP relaxations produce the same bound, since the solution of one can be used to derive a solution to the other with the same objective function value. See e.g., Section 3.3 in Taylor (2015) for a formal proof. Henceforth, we refer to the SDP relaxation as that is , where the second equality (with some abuse of notation) is meant to imply that the two relaxations have the same projection in the space of the variables.

3.1.4 SOCP Relaxation of Rectangular Formulation ().

We can further apply SOCP relaxation to the SDP constraint (8d) by only imposing the following constraints on all the submatrices of for each line ,

(12)

This is a relaxation of (8d), because (8d) requires all principal submatrices of be SDP (see e.g., Horn and Johnson (2013)). Note that (12) has a hermitian matrix, i.e., . Using the Sylvester criterion, (12) is equivalent to and . The first inequality can be written as , which is an SOCP constraint in the real domain. Thus, we have an SOCP relaxation of the rectangular formulation with the feasible region defined as

(13)

This SOCP relaxation is also proposed in the literature, see e.g., Madani et al. (2013). In Low (2014a), this relaxation is proven to be equivalent to the SOCP relaxation of DistFlow model proposed in Baran and Wu (1989a, b).

3.1.5 SOCP Relaxation of Alternative Formulation ().

The nonconvex coupling constraints (2f) in the alternative formulation (2) can be relaxed as follows,

(14)

It is easy to see that (14) can be rewritten as , which represents a rotated SOCP cone in . Note that the SOCP cone (14) is the convex hull of (2f). Using (14), we have an SOCP relaxation of the alternative formulation (2). Denote the feasible region of this SOCP relaxation as

(15)

This is the classic SOCP relaxation first proposed in Jabr (2006).

3.1.6 SDP Relaxation of Alternative Formulation ().

We can also apply SDP relaxation to the nonconvex quadratic constraints (2f) in the alternative formulation (2). In particular, define a vector , of which the first components are indexed by the set of branches , and the last components are indexed by the set of buses . Essentially, represents . Then define a real matrix variable to approximate and consider the following set of constraints:

(16a)
(16b)
(16c)
(16d)

where and . Constraints (16a) and (16b) are the usual SDP relaxation of (2f), and constraints (16c) and (16d) are used to properly upper bound the diagonal elements of , where constraint (16d) follows from applying McCormick envelopes on the squared terms . In particular, if is restricted to be between , then the McCormick upper envelope for the diagonal element is given as Denote the feasible region of this SDP relaxation of the alternative formulation (2) as

(17)

Note that this SDP relaxation of the alternative formulation is derived using standard techniques, but to the best of our knowledge, it has not been discussed in the literature.

3.2 Comparison of Relaxations

The main result of this section is Theorem 3.1, which presents the relative strength of the various convex relaxations introduced above. In order to compare relaxations, they must be in the same variable space. Also the objective function depends only on the value of the real powers. Therefore, we study the feasible regions of the above convex relaxations projected to the space of real and reactive powers, i.e. the space. We use ‘’ to denote this projection. For example, is the projection of to the space.

Theorem 3.1.

Let , , , , , and be the McCormick relaxation of the rectangular formulation (1), the McCormick relaxation of the alternate formulation (2), the SDP relaxation of the rectangular formulation (9), the SOCP relaxation of the rectangular formulation (13), the SOCP relaxation of the alternative formulation (15), and the SDP relaxation of the alternative formulation (17), respectively. Then:

  1. The following relationships hold between the feasible regions of the convex relaxations when projected onto the space:

    (18)
  2. All the inclusions in (18) can be strict.

In Section 3.2.1 we present the proof of part (1) of Theorem 3.1, and in Section 3.2.2 we present examples to verify part (2) of the theorem.

3.2.1 Pairwise comparison of relaxations

The proof of part (1) of Theorem 3.1 is divided into Propositions 3.1-3.4.

Proposition 3.1.

.

Proof.

In order to prove this result, we want to show that for any given , we can find such that is in . For this purpose, set for , , =0, , for each , and and for each . By this construction, we have , , and . Therefore, (2b)-(2c) implies that the constructed satisfy (4a)-(4b); (2d) implies (4c); (4d) is trivially satisfied since . Now to see the McCormick envelopes (4e)-(4f) are satisfied, consider . Using the bounds (4d) on , can be written as

Since for all , the above inequalities reduce to , which is implied by and the bounds (6b). Similar reasoning can be applied to verify that the other McCormick envelopes in (4e)-(4f) are satisfied. Finally, it is straightforward to see that . Therefore, the constructed is in . ∎

In fact, the above argument only relies on constraints (2b)-(2d) and (6b) in . This suggests that may be strictly contained in , which is indeed the case shown in Section 3.2.2.

Proposition 3.2.

.

Proof.

It suffices to prove that . For this purpose, we want to show that for each , there exists such that . In particular, set , , and for each . Note that and satisfy constraints (6c) by the definition of McCormick envelopes. So, it remains to verify if satisfies . This involves verifying:

(19)

The first inequality holds because and . The second inequality holds since Thus, the result follows. ∎

The next result is straightforward, however we present it for completeness.

Proposition 3.3.

.

Proof.

Note that , and forms a bijection between the variables in and the variable in . Given this bijection, (1b)-(1d) are equivalent to (8a)-(8c) and (2e) is equivalent to . As mentioned in Section 3.1.4, the hermitian condition in (12) is equivalent to and , which is exactly (14). ∎

Proposition 3.4.

.

Proof.

It suffices to show that . For this, we can show that given any , there exists a such that . In particular, we construct , where , and is a diagonal matrix with the first entries on the diagonal equal to for each and all other entries equal to zero. By this construction, (16a) is satisfied. Also , i.e. (16b) is satisfied. Note that the first entries on the diagonal of is equal to . Therefore, the bounds in (16c) are equivalent to , , where the first one follows from (2d) and the second one follows from (2d) and (14). Finally, constraint (16d) is the McCormick envelope of with the same bounds in (2d), therefore it is also satisfied by the solution of . ∎

3.2.2 Strictness of Inclusions.

Now, we prove by examples that all the inclusions in Theorem 3.1 can be strict.

  1. : Consider the 9-bus radial network with quadratic objective function from Kocuk et al. (2015). We compare the strength of these three types of relaxations in Table 1, which shows that the inclusions can be strict.

    case objective
    9-bus quadratic 89.46 69.13 52.69
    Table 1: Percentage optimality gap for , , and with respect to global optimality found by BARON.

    Here, the percentage optimality gap is calculated as , where is the global optimal objective function value found by BARON (Tawarmalani and Sahinidis 2005) and is the optimal objective cost of a particular relaxation.

  2. : Consider a 2-bus system. Since there is only one transmission line, has only one SOCP constraint

    (20)

    while the SDP relaxation has