Novel Approach Towards Global Optimality of Optimal Power Flow Using Quadratic Convex Optimization
Hadrien Godard (CEDRIC - OC, UMA), Sourour Elloumi (CEDRIC), Am\'elie, Lambert (CEDRIC), Jean Maeght, Manuel Ruiz (G-SCOP_OC)

TL;DR
This paper introduces a novel method combining SDP relaxation and branch-and-bound to solve Optimal Power Flow problems to global optimality, demonstrating promising results on various network sizes.
Contribution
It adapts the MIQCR method for OPF, achieving tight bounds and improving global solution accuracy for large-scale power networks.
Findings
Lower bound at root node equals SDP relaxation value.
Method successfully applied to networks with over a thousand buses.
Encouraging initial results on diverse OPF cases.
Abstract
Optimal Power Flow (OPF) can be modeled as a non-convex Quadratically Constrained Quadratic Program (QCQP). Our purpose is to solve OPF to global optimality. To this end, we specialize the Mixed-Integer Quadratic Convex Reformulation method (MIQCR) to (OPF). This is a method in two steps. First, a Semi-Definite Programming (SDP) relaxation of (OPF) is solved. Then the optimal dual variables of this relaxation are used to reformulate OPF into an equivalent new quadratic program, where all the non-convexity is moved to one additional constraint. In the second step, this reformulation is solved within a branch-and-bound algorithm, where at each node a quadratic and convex relaxation of the reformulated problem, obtained by relaxing the non-convex added constraint, is solved. The key point of our approach is that the lower bound at the root node of the branch-and-bound tree is equal to the…
| MIQCR | BARON | |||
|---|---|---|---|---|
| Name | Root gap | Gap-Time | Root gap | Gap-Time |
| WB2 | 2.2% | 1s | 2.2% | 1s |
| LMBM3 | 0% | 1s | 0% | 1s |
| WB3 | 0% | 1s | 0% | 1s |
| WB5 | 16.7% | 23s | 16.8% | 1.92s |
| 6ww | 0% | 1s | 0.2% | 1s |
| 9 | 0% | 1s | 0% | 1.34s |
| 9mod | 0.1% | (0.1%) | 12.4% | 5.83s |
| 14 | 0% | 1s | 100% | (21.1%) |
| 22loop | 0% | 1s | 31.6% | (31.6%) |
| 30 | 0% | 1s | 100% | (100%) |
| 39 | 0% | 2s | 100% | (100%) |
| 39mod1 | 0.1% | (0.1%) | 100% | (100%) |
| 39mod2 | 0.1% | (0.1%) | 100% | (100%) |
| 57 | 0% | 1s | 100% | (100%) |
| 89pegase | 0% | 2s | 72% | (72%) |
| 118 | 0% | 3s | 100% | (100%) |
| 118mod | 0% | 7s | 100% | (100%) |
| 300 | 0% | 10s | 100% | (100%) |
| 300mod | 0.1% | (0.1%) | 100% | (100%) |
| 1354pegase | 0% | 204s | 69% | (69%) |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Novel Approach Towards Global Optimality of Optimal Power Flow Using Quadratic Convex Optimization
Hadrien Godard1,2,3, Sourour Elloumi2,3, Amélie Lambert2, Jean Maeght1 and Manuel Ruiz1 1 R&D Division, RTE (French TSO), 9 rue de la Porte de Buc, 78000, Versailles, France2CEDRIC, CNAM, 292 rue Saint-Martin, 75003, Paris, France3UMA, ENSTA, 828 Boulevard des Maréchaux, 91120, Palaiseau, France
Abstract
Optimal Power Flow (OPF) can be modeled as a non-convex Quadratically Constrained Quadratic Program (QCQP). Our purpose is to solve OPF to global optimality. To this end, we specialize the Mixed-Integer Quadratic Convex Reformulation method (MIQCR) to (OPF).
This is a method in two steps. First, a Semi-Definite Programming (SDP) relaxation of (OPF) is solved. Then the optimal dual variables of this relaxation are used to reformulate OPF into an equivalent new quadratic program, where all the non-convexity is moved to one additional constraint. In the second step, this reformulation is solved within a branch-and-bound algorithm, where at each node a quadratic and convex relaxation of the reformulated problem, obtained by relaxing the non-convex added constraint, is solved. The key point of our approach is that the lower bound at the root node of the branch-and-bound tree is equal to the SDP relaxation value.
We test this method on several OPF cases, from two-bus networks to more-than-a-thousand-buses networks from the MATPOWER repository. Our first results are very encouraging.
I Introduction
The Optimal Power Flow (OPF) problem deals with determining power production at different nodes of an electric network where a production cost is minimized. If we model this network with an AC framework, the optimization problem is quadratic and non-convex [16].
(OPF) can be solved with general purpose solvers such as Baron [21]. A branch-and-bound algorithm specialized to (OPF) has been introduced in [7] using the SDP rank relaxation [15] as a lower bound provider. More recently [13] introduces an SOCP-based branch-and-bound algorithm.
Our goal is to design a branch-and-bound algorithm that closes the gap between lower and upper bounds. The lower bound is obtained with the rank relaxation and the upper bound is obtained with a feasible point computed, for instance, by an interior point method. It has been already observed that this gap is quite small for (OPF) instances as these lower and upper bounds are very sharp in general [9].
To this end, we work on a specialization of MIQCR [2, 3, 6] which is a method designed to solve non-convex and mixed integer quadratic programs to global optimality. MIQCR works in two steps: first a semi-definite relaxation is used to reformulate the problem into another quadratic program. Then this reformulated problem is solved within a branch-and-bound framework where at each node a quadratic and convex program is solved to get a local lower bound. A key advantage of this method is that it requires the solution of only one SDP relaxation as a preprocessing step. Then, the strength of this SDP lower bound is captured onto quadratic and convex programming. However, the solution of this large SDP is the bottleneck of MIQCR when handling larger instances. The contribution of this paper is a specialization of MIQCR to (OPF) where we prove that we can reach the strengthened lower bound of the MIQCR method by solving a smaller semidefinite relaxation than in the original approach.
In this paper, we adapt method MIQCR to the OPF problem. In Section 2, we recall the formulation of the OPF problem as a quadratic program. In Section 3 we introduce the semi-definite relaxation used in our algorithm, known as the rank (or Shor) relaxation of the OPF problem [15]. We also give in this section a new proof of the strong duality of the rank relaxation. Then, in Section 4 and in Section 5, we present our contribution: the specialization of MIQCR to solve the OPF problem. In Section 4, we prove that solving a smaller SDP relaxation than in the original MIQCR approach is sufficient to reformulate the OPF problem and to attain the sharp rank relaxation bound at the root node of the branch-and-bound tree. Indeed, solving the dual of the rank relaxation is equivalent to finding the best quadratic reformulation of an OPF problem among the MIQCR family of reformulations. Moreover, in Section 5, we present a branch-and-bound framework to solve the reformulated problem. Finally, in Section 6, we illustrate our method by computational experiments on small and medium-sized instances of OPF problems. Section 7 draws a conclusion.
Notations
- •
i is the complex number whose real part is null and imaginary one equals one.
- •
is the set of symmetric matrices on .
- •
is the zero matrix of size .
- •
is the identity matrix of size .
- •
is the smallest eigenvalue of a symmetric matrix .
- •
is the optimal value of an optimization problem .
- •
is the canonic scalar product between matrices and .
- •
is the cardinal of a set .
II A quadratic formulation of the OPF problem
Driving power flows from producers to consumers in an electric network constitutes the OPF problem. Usually the amount of consumed electric power is known at each node of the network. On the contrary the production is unknown and OPF deals with determining its value. The goal is to minimize electric power production costs under the constraint that the demand is satisfied at each node, and that active and reactive powers are box constrained for each production unit. Because real electric transmission networks work with alternative current, one must consider voltage at each node in order to compute where the power flows through the network. Engineering limits such as bounds on voltage magnitude are also considered as constraints. See the bus injection model in [16] for more precisions on this topic.
We only consider a cost linearly linked to the power produced at each unit. Note that a linear cost on the power production handles the case where one tries to minimize losses on the network. The following program is a classic model for the (OPF) problem:
[TABLE]
Where is the set of network nodes and their number. , respectively , is the set of production, respectively consumption, nodes and . Variables are the complex voltages of network buses. Variables are the generated active powers at production buses. Variables are the generated reactive powers at production buses. is the vector of linear costs, where is the production cost at node . is the complex admittance matrix at node . and are lower and upper bounds on active generated powers. and are lower and upper bounds on reactive generated powers. and are lower and upper bounds on voltage magnitudes.
The objective function (1) is linear relatively to active produced powers . Constraints (2) (resp. (3)) are power balance equations at production (resp. consumption) nodes. Constraints (4), (resp. (5)), are bounds on active (resp. reactive), produced powers. Constraints (6) are bounds on voltage magnitudes.
Substituting voltages to other variables, from equations (2), one can see that (OPF) can be modeled as a pure quadratic program, without linear terms, in the real and imaginary parts of voltages.
To simplify the notations, we rewrite the above problem in a more abstract style as the following quadratically constrained quadratic program with variables :
[TABLE]
Variables are the real and imaginary parts of the voltage at each network node. . is the set of constraints indices. At each node there are two constraints that bound the voltage magnitude, and four other constraints that model the complex power balance. Hence, the number of constraints is . and .
The objective function and the constraints are convex if and only if matrices and are positive semi-definite (PSD) which is not the case for .
In this formulation with continuous variables there are no natural lower and upper bounds on variables . However to perform a spatial branch-and-bound algorithm, initial lower and upper bounds ( and ) on each variable are needed.
To obtain such bounds one can use the fact that the modulus of each complex voltage is upper-bounded. Those upper-bounds are also at the heart of the proofs of Propositions 1 and 2.
Suppose that is the real part of the complex voltage at node , and its imaginary part. By Constraints (6), we have:
[TABLE]
It follows that:
[TABLE]
Below, we take:
[TABLE]
III The rank relaxation of OPF
In this section we recall the rank relaxation of , that we call :
[TABLE]
The dual of is:
[TABLE]
where is the dual variable associated with constraint .
Strong duality for feasible OPFs has been already proved (see for instance [8]). In this paper, we propose another proof based on the fact that the modulus of each complex voltage is bounded. Parts of this proof will be used later to demonstrate Proposition 2.
Proposition 1
If is feasible, there is no duality gap between and . In other words, strong duality holds for the rank relaxation of feasible OPF problems.
Proof:
We first prove that if is feasible then is feasible too. Indeed from each feasible solution to one can build a feasible solution to .
Let us now prove that is strictly feasible, i.e. finding satisfying .
In the following we assume that the elements of are integers from 1 to , and that the first elements of are the indices of voltage magnitude upper-bound constraints (7) on the network nodes. For from to : all entries of are zeros except the -th and -th entries of the diagonal, which are equal to 1. It follows that .
For : take . Consider the matrix
.
For from to :
take .
Then:
[TABLE]
As has positive entries and , is strictly feasible.
To sum up, is feasible, is strictly feasible, so strong duality holds. ∎
We have shown that strong duality holds for the rank relaxation of , using the fact that the modulus of each complex voltage is bounded. This is a special case of ball constraints on an optimization problem with complex variables, see [8] to get more details on this subject.
relaxation of is known to give a sharp lower bound. We want to point out that we study transmission networks which are not tree networks and are highly meshed. Moreover we do not have a sufficient number of phase-shifter-transformers on network to use results from [20]. Thus rank relaxation does not necessarily lead to an optimal solution.
In the next sections, we present a branch-and-bound algorithm that starts with this sharp lower bound which is based at each node on a quadratic and convex relaxation of .
IV An equivalent formulation to
The first step of MIQCR consists in reformulating a QCQP like into an equivalent quadratic problem that has a quadratic and convex objective function, linear constraints, and additional variables that are meant to satisfy quadratic constraint .
Let be a positive semi-definite matrix. We reformulate as:
[TABLE]
Observe that the reformulated objective and constraints functions have the same value as the original ones for a same if Constraints (10) are satisfied. Moreover, the new objective function is convex since matrix is positive semi-definite. The constraints are now linear, and thus convex. This is why, this reformulation is called a quadratic and convex reformulation.
In only Constraint (10) is non convex. In a way, all the non-convexity has been moved into this constraint.
To solve to global optimality MIQCR uses a branch-and-bound algorithm where, at each node, we relax Constraint (10) and add the linear McCormick inequalities (11)-(14) [17] to tighten the relaxation.
Therefore, at the first node of the branch-and-bound, the quadratic and convex relaxation is solved.
[TABLE]
where are defined as in (8) and (9) and
E=\Big{\{}(i,j)\in\{1,\ldots,2n\}^{2}\,:\,i\leq j\Big{\}}.
Notice that in practice, the full variables matrix is not considered. Coefficient is considered if and only if the entry is non zero in a matrix among .
Every positive semi-definite matrix gives a different reformulation. In the particular case where , the objective function is linear. It is the linearization of . For example the Baron solver [21] relies on this linearization within a branch-and-bound framework.
We are now interested in finding a ”best” matrix , i.e. a matrix which gives the largest lower bound at the root node of the branch-and-bound tree. That is to say, a matrix which maximizes the value of . More formally we want to solve:
[TABLE]
It is proved in [3] that a best matrix can be computed from optimal dual variables of a semi-definite relaxation of which is but with additional constraints and variables to raise the McCormick’s inequalities.
In this paper, we characterize , which does not contain the McCormick inequalities, as the relaxation that can be used to compute a best positive semi-definite matrix. The fact that the McCormick inequalities are redundant in the relaxation is a significant result and is the main difference between Proposition 2 and the result in [3]. As a consequence, to reformulate , we have to solve a problem with a smaller size than the one solved in the original method MIQCR. Moreover, we prove that the optimal value of this ”best” quadratic and convex relaxation is equal to the optimal value of .
Proposition 2
Let be an optimal solution to , take:
[TABLE]
If and are defined by (8) and (9) we have:
[TABLE]
Proof:
For each PSD matrix , let us introduce the optimization problem:
[TABLE]
is the relaxation of where inequalities (11)-(14) have been dropped. Thus, it is a relaxation of .
The proof is divided in two parts.
- •
First we prove that:
[TABLE]
Let us rewrite the dual of the SDP relaxation by introducing a slack matrix :
[TABLE]
For a given , the optimization problem in is a linear program. We replace it by its LP-dual and obtain the equivalent problem:
[TABLE]
Now we can observe that as S is positive semi-definite, one can add to the objective function together with variables . Indeed will be equal to 0 in any optimal solution. Therefore:
[TABLE]
We have proved:
[TABLE]
- •
Now, in the second part, we prove that:
[TABLE]
From the first part, and by Proposition 1, it follows that , and, as is a relaxation of :
[TABLE]
Let us now prove that:
[TABLE]
Let be a solution to . Let us show that is a feasible solution to with a lower objective value.
The inequalities are trivially satisfied. Let us now prove that (11)-(14) are satisfied. Which amounts to prove:
[TABLE]
When , as , then for all . Moreover, from (7):
[TABLE]
And therefore:
[TABLE]
When , as . From (17), it follows that . Then (15) and (16) are satisfied.
Let us now compare the objective solution values of and :
[TABLE]
as and are PSD.
Therefore has a lower solution value than .
To sum up for any PSD matrix , has a lower optimal solution value than .
We can now conclude that and that maximizes the value of . ∎
Remark 1
In the proof above, we demonstrate that any solution of ”satisfies” the McCormick inequalities. This is because bounds and were obtained with constraints (7) which are in . See Section II.
Remark 2
The optimal matrix is not unique, other matrices may give a root node relaxation with the same value.
Proposition 2 ensures that the lower bound obtained at the root node of our branch-and-bound framework is equal to the rank relaxation bound. For many test cases, this bound seems to be very sharp [9]. To solve and compute matrix, one can use the solver introduced in [9, 18].
Those results allow us to build the following algorithm to solve the OPF problem to global optimality:
Solve the rank relaxation and deduce optimal dual variables .
Define the PSD matrix .
Solve within a branch-and-bound algorithm.
V Solution within a branch-and-bound algorithm
In the previous section we showed how to build an ”optimal” reformulation of in the sense that it maximizes the lower bound at the root node of our branch-and-bound tree. In this section we describe the second step of the MIQCR method: the solution within a branch-and-bound algorithm. This algorithm is used to solve and hence .
Let us recall that a branch-and-bound is an enumeration tree used to solve an optimization problem. Each node of the tree represents a sub problem of the original one. There are multiple ways to divide the original problem into sub problems. The classic way is to divide each variable interval into different subintervals, those branch-and-bound algorithms are called spatial branch-and-bound. We choose to implement this type of branch-and-bound, that is why we need bounds and on variables . To sum up, at each node we modify values of and to build the sub problem.
A branch-and-bound implementation is defined by:
- •
Actions performed at each node of the tree,
- •
Next node selection strategy.
V-A How to deal with a node ?
At each node we solve where bounds and are different. This change modifies the relaxation value since and are involved in the McCormick inequalities (11)-(14). We recall that this relaxation is convex and quadratic and that it gives a lower bound of the node subproblem.
Next step depends on the result of the node relaxation:
- •
If the relaxation is infeasible: the branch is pruned.
- •
If the lower bound from relaxation is greater than the best current upper bound: the branch is pruned.
- •
If the solution from the relaxation satisfies constraint (10): is then a solution of , the branch is pruned and the upper bound is potentially updated.
- •
Else: two new nodes are built as children of the current node.
About branching: To build two children nodes from a parent node, a branching variable and a branching value are chosen.
Variable selection strategy: is chosen among variables that violate the most Equality (10) (for the Euclidian norm).
Interval division: is chosen between the middle of the current interval of variable (denoted ) and , the value of in the node relaxation solution. Let be a parameter between 0 and 1:
[TABLE]
Lower bounds on the local optimal solution value are computed at each node, however we need to find upper bounds and we cannot rely on relaxations results to do so. That is why every three nodes, a heuristic, consisting in finding a point satisfying first order optimality conditions of , is launched.
V-B How to select the following node ?
At each node is attached a potential which is the lower bound found at its parent node. At this node the relaxation value cannot be lower than this potential.
Our next node selection strategy is designed as a ”best-first” strategy. Indeed we want to see the global lower bound of the tree (which is the minimum among the potentials of the leafs in the tree) increase. To do so, we select the node with the lowest potential as the following node to handle.
Branch-and-bound algorithm is terminated when the relative difference between the lowest upper bound (the best feasible solution) and the global lower bound is less than an -value.
VI Numerical experiments
To illustrate our method, we present results on a batch of 20 OPF instances involving networks from 2 to 1354 nodes. All these instances come from the MATPOWER repository [23], except the ”WB” and ”LMBM” instances that come from [4].
In our experiments (as in [9]) we do not consider any constraints on current magnitudes, and we only consider the linear part on the active power cost.
For each instance, we launch our implementation of MIQCR along with the Baron solver [21], keeping default options.
For MIQCR, the rank relaxation is solved with the Mosek solver [19] and a chordal decomposition to exploit the problem sparsity. The quadratic and convex relaxation at each branch-and-bound node is solved with the Xpress [22] solver. Heuristic performed at each three nodes is based on the solution of the first order optimality conditions of with the Knitro solver [5, 12]. To perform computation we use the following parameters:
[TABLE]
Hence a 0% gap means that the gap is under the parameter value. The choice of came after some numerical experiments along with the reading of [1].
Table 1 presents numerical results where each line refers to one OPF test case. In the first column, each instance name contains the number of nodes in the associated electric network. We recall that in the number of real variables is twice the number of nodes. For each solver, the Root gap column gives the relative gap between the best known solution for the instance and the lower bound found at the root node of the branch-and-bound tree, i.e., if the best known solution is denoted by and the lower bound by , the gap equals . The Gap-Time column gives the execution time if global optimality is reached within five minutes of computation, if this is not the case, this column gives the relative final gap between the best known solution and the final lower bound.
On small test cases (under 10 nodes) the gap is closed within the branch-and-bound tree by both solvers, except for 9mod instance where MIQCR fails. We observe that the root gap of MIQCR is very tight. This reflects the quality of the SDP rank relaxation lower bound as already observed by several authors [9, 15]. On the contrary, the root gap obtained with the Baron solver [21] is often very large for instances with more than 10 nodes. Moreover its root gap is not improved during the five minutes of branch-and-bound computation.
Notice that a feasible solution is always found for each test case by MIQCR, and may be improved within the branch-and-bound. For instance, for case300 first upper bound found equals 475783 and it is improved to 475462.2 during branch-and-bound tree iterations.
In [10], authors successfully use the Lasserre hierarchy [14] to also solve smallest instances of (OPF) problems to global optimality. Although their method can be extended to deal with some well-conditioned larger instances [8, 11], it fails with largest generic ones.
VII Conclusion
In this paper we show how to adapt the MIQCR method to the OPF problem. We prove that the well-known rank relaxation is sufficient to build an ”optimal” quadratic reformulation of . This result can be extended to each quadratic problem with complex variables whose magnitudes are upper-bounded.
When solving this optimal reformulation within a branch-and-bound framework, we can close the gap between rank relaxation and known feasible solutions.
First numerical results are encouraging. Future work consists in a specialization of the branch-and-bound framework in order to close the gap by raising the lower bound in larger instances of the OPF problem. To do so we will focus on bound tightening techniques.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] P. Belotti, J. Lee, L. Liberti, F. Margot and A. Wächter. Branching and bounds tightening techniques for non-convex MINLP. Optimization Methods & Software , 24(4-5):597–634, 2009.
- 2[2] A. Billionnet, S. Elloumi, and A. Lambert. Extending the QCR method to general mixed-integer programs. Mathematical Programming , 131(1):381–401, 2012.
- 3[3] A. Billionnet, S. Elloumi, and A. Lambert. Exact quadratic convex reformulations of mixed-integer quadratically constrained problems. Mathematical Programming , 158(1):235–266, 2016.
- 4[4] W. A. Bukhsh, A. Grothey, K. I. M. Mc Kinnon and P. A. Trodden. Local Solutions of the Optimal Power Flow Problem. IEEE Transactions on Power Systems , 28(4):4780–4788, 2013.
- 5[5] R. H. Byrd, J. Nocedal, and R. A. Waltz. Knitro: An integrated package for nonlinear optimization. Springer Verlag , 2006, pp. 35–59.
- 6[6] S. Elloumi and A. Lambert. Global solution of non-convex quadratically constrained quadratic programs. Optimization Methods and Software , 34(1):98–114, 2019. https://doi.org/10.1080/10556788.2017.1350675 . · doi ↗
- 7[7] A. Gopalakrishnan, A. U. Raghunathan, D. Nikovski and L. T. Biegler. Global optimization of Optimal Power Flow using a branch & bound algorithm. Communication, Control, and Computing (Allerton) , 609–616, 2012.
- 8[8] C. Josz. Application of polynomial optimization to electricity transmission networks. Ph D. Thesis , 2016. http://arxiv.org/abs/1608.03871 .
