A tight compact quadratically constrained convex relaxation of the Optimal Power Flow problem

In this paper, we consider the optimal power flow (OPF) problem which consists in determining the power production at each bus of an electric network by minimizing the production cost. Our contribution is an exact solution algorithm for the OPF problem. It consits in a spatial branch-and-bound algorithm based on a compact quadratically constrained convex relaxation. This compact relaxation is computed by solving the rank relaxation once at the beginning of the algorithm. The key point of this approach is that the lower bound at the root node of the branch-and-bound tree is equal to the rank relaxation value, but is obtained by solving a quadratic convex problem which is much faster than solving a SDP. To construct this compact relaxation, we add only O(n) variables that model the squares of the initial variables, where $n$ is the number of buses in the power system. The relations between the initial and auxiliary variables are therefore non-convex. By relaxing them in our relaxation, we have only O(n) equalities to force by the branch-and-bound algorithm to prove global optimality. Our first experiments show that our new algorithm Compact OPF (COPF) performs better than the methods of the literature we compare it with


Introduction
The Optimal Power Flow (OPF) problem consists in the determination of the power production at different buses of an electric network that minimizes a production cost.The electrical transmission network is modeled by a graph G = (N , E).Each network point belongs to the set N of nodes (i.e the set of buses), and their connections (i.e. the set of transmission lines) are modeled by the set of edges E. We assume that there is an electric demand at each node also called load.We distinguish two classes of nodes: N = N g ∪ N d , where N g is the set of nodes that generates and flows the power (the generator nodes), and N d is the set of nodes that only flows the power (the consuming nodes).The aim of the OPF problem is to satisfy demand of all buses while minimizing the total production costs of the generators such that the solution obeys Ohm's law and and Kirchhoff's law.This problem is naturally formulated with complex variables.Let Y ∈ C n×n , where |N | = n be the admittance matrix, which has component Y i,k = R i,k + jI i,k for each line (i, k) of the network, and R i,i = r i,i − i =k R i,k , I i,i = i i,i − i =k I i,k , where r i,i (resp.i i,i ) is the shunt conductance (resp.susceptance) at bus i, and j 2 = −1.Let p i , q i be the real and reactive power output of the generator node i, and p i , q i the given real and reactive power output of the load node i.We consider here the complex voltage in the rectangular form: i + f 2 i is the voltage magnitude, and we denote by δ(i) the set of adjacent nodes of bus i.With the above notation, the OPF problem can be modeled by the well known rectangular formulation of [27]: −q i = −I i,i (e where C ∈ S + |Ng| is a diagonal and semi-definite matrix, c ∈ R |Ng| is the vector of linear costs of the power injection at each generator node, (v, v) ∈ (R n , R n ) are the bounds on the voltage magnitude, and (p, p, q, q) ∈ (R |Ng| , R |Ng| , R |Ng| , R |Ng| ).This formulation has 2(n + 2|N g |) variables, 2n quadratic equalities (1)-( 4) that enforces the active and reactive power balances at each node, 2n quadratic inequalities (5) that models the voltage magnitude, and 4|N g | box constraints (6)- (7).
The first results of the literature for solving the OPF problem were focused on optimal local solutions, mostly by adapting interior point methods, see, e.g., [30,27,15,29].In the context of global optimization, one requires furthermore to determine lower bounds on the OPF problem.For this, the second-order cone programming (SOCP) and the semidefinite programming (SDP) relaxations were first used (see [13,3,2,19]).The most used SDP relaxation, also named the rank relaxation, leads to very tight lower bounds on the OPF problem.In particular, it was proven in [9] that this relaxation is exact for a restricted class of problems and under some assumptions.In other cases, where it is not exact, it can be strengthen until the optimal solution value, following the ideas of the hierarchy of moment relaxation problems [18,24], that can be applied to any polynomial optimisation problem.This approach was specialized in the context of the OPF problem in [16], and showed its efficiency to solve small-size problems.It was also used in [22] to strengthen the lower bounds for larger problems.Unfortunately, in practice, using interior point methods for solving several large SDP relaxations, which sizes increase at each rank of the hierarchy, is intractable for large networks.Several specialized algorithms that exploit the sparsity of power networks were thus proposed as in [13,14,23,22,20].
More recently, several cheaper computable convex relaxations were introduced for the OPF problem.For instance, linear and quadratic envelopes for trigonometric functions in the polar formulation of the OPF problem are constructed in [6,4,5], and strong SOCP relaxations were introduced in [17].These polynomial bounds can then be used within a spatial branch-and-bound framework to solve the problem to global optimality.
Another exact solution approach called RC-OPF was proposed in [10,7].This method is a specialization of the Quadratic Convex Reformulation approach proposed in [8] for generic quadratic problems.It consits in a branch-and-bound algorithm based on the computation of a tight quadratic convex relaxation the OPF problem that is computed by solving the rank relaxation once at the beginning of the algorithm.The key point of this approach is that the lower bound at the root node of the branch-and-bound tree is equal to the rank relaxation value, but is obtained by solving a quadratic convex problem which is much faster than solving a SDP.This tight quadratic convex relaxation has a quadratic convex objective function and linear constraints.Moreover, it has a size of O((2n) 2 ), since it relies on the introduction of (2n) 2 additional variables that model all the possible products of the original variables.Then a spatial branch-and-bound is performed which aim is to enforce the equality between each additional variable and its corresponding product.Unfortunately, and as illustrated in Section 5, in practice enforcing these (2n) 2 equalities can be very time consuming for the OPF problem.
Our contribution is an exact solution algorithm for the OPF problem that is based on a compact quadratic convex relaxation.This relaxation has only O(n) auxiliary variables and constraints.Moreover, as in method RC-OPF, its optimal value also equals to the optimal value of the well known rank relaxation.We thus get a cheaper computable convex relaxation than is method RC-OPF, that is as tight as the rank relaxation.Finally, to solve (OP F ) to global optimality, we perform a spatial branch-and-bound algorithm based on this new quadratic convex relaxation.A key advantage of our compact relaxation is that we only have to enforce the equality between 2n additional variables and their corresponding products to prove global optimality.Finally, we evaluate our method Compact OPF (COPF) on several instances of the literature and we show that our algorithm outperforms method RC-OPF and the generic global optimization solver Baron [26].
The paper is organized as follows.In Section 2, we describe how we convexify the quadratic constraints of (OP F ).Then, in Section 3, we present a new family of equivalent quadratic convex formulations to (OP F ).In Section 4, we show how we calculate the best quadratic convex relaxation within this family.Finally, in Section 5, we present some computational results.Section 6 draws a conclusion.

A compact convexification of the quadratic constraints
We start by observing that the structure of the formulation (OP F ) is specific.First, only variables p are involved into the objective function.Moreover, the matrix C is diagonal and positive semidefinite, hence the objective function is convex and separable.It follows that the non convexities only come from the quadratic constraints (1)-( 5).More precisely, in the constraints, variables e and f are only involved into quadratic forms, while variables p and q only in linear forms.
Starting from the latter observations, our idea is to build an equivalent problem to (OP F ), where the original constraints are convexified thanks to 2n auxiliary variables z = (z e , z f ) ∈ R 2n that model the squares of the initial variables e and f : Using these new variables it is easy to rewrite Constraints (5) into a convex form by simply linearizing them, and we get: We now handle the equality constraints.For this, the first step is to transform each equality ( 1)-( 4) into two inequalities.We thus introduce for all i ∈ N g , Constraints ( 12)-( 19) obviously equivalent to Constraints (1)-( 4): Then, recall that Inequalities ( 12)-( 19) are only quadratic on variables e and f .To make them convex, we apply the smallest eigenvalue method introduced in [11].Denote by i ) the sub-matrix of the Hessian of the i th Constraints (1) (resp.( 2), ( 3), ( 4)) that corresponds to the quadratic terms involving variables e and f only.Let λ(A 1 i ) be the smallest eigenvalue of matrix A 1 i , and d(λ(A 1 i )) be the diagonal matrix where each diagonal term equals λ(A 1 i ).To rewrite Inequality (12) (resp.( 13)) as a convex function, we add the quadratic quantity −λ( to it.As a consequence the Hessian matrix of the new function is ) that is obviously a positive semidefinite matrix.Moreover, the value of the convexified function remains the same as soon as for all k ∈ N , e 2 k + f 2 k − z e k − z f k = 0, or equivalently when Equalities ( 9)-( 10) are satisfied.We thus obtain the set of convex Constraints ( 20)-( 27): By replacing Constraints (1)-( 5) by Constraints ( 9)-( 11),( 20)-( 27), we then obtain an equivalent problem to (OP F ) where the only non-convexity remains into the equalities ( 9)- (10).Then, we use the McCormick's envelopes [21], to relax the latter equalities and we get a quadratic convex relaxation of (OP F ), that can then be used for the bounding step of a classical spatial branch-and-bound algorithm.Performing such an algorithm is highly dependant on the quality of the bound at the root node.Moreover, we know that for the OPF problem the rank relaxation provides sharp bound.This is why in the rest of the paper, we focus on the computation of quadratic convex relaxation of (OP F ) whose value equals to the optimal value of the rank relaxation.

Building a compact family of equivalent quadratic formulations
We now present a compact family of quadratic reformulations to (OP F ).For simplicity, we start by rewriting the initial equality constraints ( 1)-( 4) by using the notation y = (p, q) ∈ R 2|Ng| and x = (e, f ) ∈ R 2n as follows: where C = (N g , N g , N d , N d ), with |C| = 2n, and ∀ r ∈ C, A r ∈ S 2n is the Hessian matrix of constraint r (i.e.matrices A 1 i , A 2 i , A 3 i , and A 4 i ), a r ∈ R 2|Ng| is the vector of linear coefficients of constraint r, and b ∈ R 2n , where coefficient b r is the the right-hand side of constraint r.
Let (φ, γ) ∈ (R 2n , R 2n ) be two vector parameters, we build the following parameterized function: where h(y) is the initial objective function, and we recall that for all i ∈ N , x i = e i , x i+n = f i , z i = z e i , and z i+n = z f i .Observe that there exist parameters (φ, γ) such that h φ,γ is a convex function.Indeed, as mentioned above, function h(y) is convex and separable.Now, the two additional terms are linear in y and z, and pure quadratic in x.By taking ∀ r ∈ C, φr = 0, and ∀i ∈ N , γi any non negative value, the associated function hφ ,γ (x, y, z) is obviously convex.
By denoting λ r = λ min (A r ) and λ r = λ min (−A r ), we are now able to build a family of equivalent formulations to (OP F ): where Constraints ( 28)-( 29) are Constraints ( 20)-( 27) written in a compact form.It is easy to see that problem (OP F φ,γ ) is equivalent to problem (OP F ), in the sense that any optimal solution from one is an optimal solution from the other.Indeed, h φ,γ (x, y, z) = h(y) when Constraints ( 1)-( 4), (9), and (10) are satisfied.Moreover, we already observed that we can choose parameters such that the objective function h φ,γ is a convex function.Finally, the only constraints that remain non-convex are Constraints ( 9) and (10).To derive a quadratic convex relaxation of (OP F φ,γ ), we relax them by use of the McCormick's upper envelopes (see [21]), keeping the quadratic convex lower envelope.However, for this, we need upper and lower bounds on each variable x i .Some trivial initial bounds can easily be deduced from Constraints (5), i. e. − √ By denoting with ( , u) ∈ (R 2n , R 2n ) these bounds (i.e.= (− , we replace Constraints ( 9) and ( 10) by the following set of convex inequalities: We get (OP F φ,γ ), a family of quadratic convex relaxations to (OP F ): Problem (OP F φ,γ ) is a compact quadratic convex relaxation to (OP F ), since we only add O(n) variables and constraints to the original formulation.

Computing a sharp quadratic convex relaxation
We are now interested in the best parameters (φ * , γ * ) that maximize the optimal value of (OP F φ,γ ) while making convex the parameterized function h φ,γ .We prove that these best parameters can be deduced from the dual optimal solution of the rank relaxation of (OP F ).For simplicity, we denote by γ = (γ, γ) ∈ R 2n , and with this notation, we can rewrite function h φ,γ as follows: where d(γ ) is the diagonal matrix of dimension 2n, where the i th -diagonal coefficient equals γ i .We formally pose the problem we aim to solve as follows: (P ) : max v(OP F φ,γ ) : where v(OP F φ,γ ) is the optimal value of problem (OP F φ,γ ).
We state in Theorem 1 that the optimal value of (P ) equals the optimal value of the following the so-called rank relaxation of (OP F ): (SDP ) where Constraints (32) are the compact form of the following set of constraints: The optimal value of (P ) equals the optimal value of (SDP ). Proof.
To prove that v(P ) ≤ v(SDP ), we show that for any feasible solution ( φ, γ ) to (OP F φ,γ ), we have v(OP F φ,γ ) ≤ v(SDP ), which in turn implies that v(P ) ≤ v(SDP ) since the right hand side is constant.Let (ȳ, Ȳ , X) be a feasible solution of (SDP ), and build the solution (x = 0, y = ȳ, z = D( X)), where D(X) is the vector composed of the diagonal terms of matrix X.We show that: i) (x, y, z) is feasible for (OP F φ,γ ), and i,i) its objective value is less or equal than v(SDP ).Since (OP F φ,γ ) is a minimization problem, (OP F φ,γ ) ≤ v(SDP ) follows.
By constraints (32), an by definition of γ , we have: Indeed, the first term is non positive since C i,i ≥ 0, and by Constraint (37), we have that ∀i, ȳ2 i − Ȳi,i ≤ 0. This is also the case for the last term, by Constraint (38), and since r∈C φr A r + d(γ ) 0 by feasibility of (P ).
Denote by ∆ the difference between the objective values, and we below prove that ∆ ≥ 0.
Hence, for computing the best convex relaxation of the OPF problem, we need to solve once problem (SDP ).Then, to solve (OP F ) to global optimality, we perform a spatial branch-and-bound algorithm based on the bound given by problem (OP F φ * ,λ * ).As proven in Theorem 1, the bound at the root node of the tree equals the optimal value of (SDP ).
Another advantage of our compact relaxation, in comparison to the complete linearization of the constraints of method RC-OPF, is that we only have to enforce 2n equalities to prove global optimality in our branch-and-bound.We sum up our global optimization algorithm Compact OPF (COPF) in Algorithm 1.

Some experimental results
In this section, we illustrate on some experiments the behavior of the algorithm COPF for exact solution of instances of the PG-lib library [25].We compare COPF with the nonlinear solver Baron [26], and with the initial quadratic convex reformulation approach RC-OPF [10,7].Our experiments were carried out on a server with 2 CPU Intel Xeon each of them having 12 cores and 2 threads of 2.5 GHz and 4 * 16 GB of RAM using a Linux operating system.We set the time limit to 3 hours for all methods.For the solver Baron, we use the multi-threading version of Cplex 12.9 [12] with up to 64 threads.For methods COPF and RC-OPF, we used the semidefinite solver Mosek [1] for solving semidefinite programs.At each node of the spatial branch-and-bound, we used the solver Mosek for solving the QCQP of method COPF, and the solver Cplex 12.9 for solving QP of method RC-OPF.For computing feasible local solutions, we use the local solver Ipopt [28].
For our experiences, we considered medium-sized data of power networks having 3 to 300 buses, and we took the formulation of the OPF described by Constraints (1)- (7).We report in Table 1 the characteristics of each instance: its Name, and the number of Buses, Generators, and Lines of the considered power network.We also indicate in the column Opt the best solution found by COPF and RC-OPF, within 3 hours of computing time.The column |(y, x)| specifies the dimension of the variable vector (y, x) in (OPF).
RC-OPF does not increase the lower bound, despite a large number of explored nodes.This is not the case for COPF, which, with a much smaller number of nodes, slightly increases the lower bound over the course of the branch-and-bound, even for the largest instances.This is because the number of auxilliary variables and relaxed equalities z i = x 2 i is strongly reduced in COPF.Let us finally note that determining a feasible solution is also difficult in practice, at least with a generic interior point algorithm.Moreover, the CPU time for this step is not constant since for example it goes from 400 to 1380 seconds for the instance pglib opf case200 activ.

Conclusion
We consider the OPF problem that determines the power production at each bus of an electric network minimizing a production cost.In this paper, we introduce a global optimisation algorithm that is based on a new quadratically constrained quadratic relaxation.This relaxation is compact in the sense that it has only O(n) auxiliary variables and constraints, where n is the number of buses of the network.We moreover prove that our quadratic relaxation has the same optimal value as the rank relaxation.Finally, to solve (OP F ) to global optimality, we perform a spatial branch-and-bound algorithm based on our new quadratic convex relaxation.Another advantage of our approach is that to prove global optimality, we have a reduced number of non-convex equalities to force in spatial the branch-and-bound.We report computational results on instances of the literature.These results show that this new approach is more efficient than the standard solver Baron.
Moreover, for the most difficult instances, and with a basic implementation, it is able to improve the lower bound unlike other methods in the literature.A future work consists in using OBBT techniques to further improve the behaviour of COPF on the most difficult instances.