Improving Clique Decompositions of Semidefinite Relaxations for Optimal Power Flow Problems
Résumé
SSemidefinite Programming (SDP) provides tight lower bounds for Optimal Power Flow problems. However, solving largescale SDP problems requires exploiting sparsity. In this paper, we experiment several clique decomposition algorithms that lead to different reformulations and we show that the resolution is highly sensitive to the clique decomposition procedure. Our main contribution is to demonstrate that minimizing the number of additional edges in the chordal extension is not always appropriate to get a good clique decomposition.
1 Introduction
The Alternating Current Optimal Power Flow (ACOPF) problem is a nonconvex optimization problem whose purpose is to provide an operating point of the power network minimizing generation costs. Solving this problem is essential for RTE, the French transmission system operator, because ACOPF problems are subproblems of decision problems involving millions of euros, e.g. operational planning and grid development. Yet there is still no efficient method to solve ACOPF problems to global optimality.
Many convex relaxations have been proposed to evaluate the quality of feasible solutions computed with local methods [13]. SDP relaxations often provide tight lower bounds [12] that are useful to prove global optimality. Some promising global optimization methods are based on these relaxations: Godard et al. propose an adaptation of the MixedInteger Quadratic Convex Reformulation method to ACOPF problems in [6], Gopalakrishnan et al. present a branchandbound approach using SDP relaxations in [7] and Josz et al. apply the Lasserre hierarchy in [11] to achieve global optimality. All these methods depend on largescale SDP problems being solved efficiently.
Many algorithms are proposed in the literature to solve SDP problems. In this paper, we focus on interiorpoint methods that are the most reliable and accurate algorithms as far as we know. Interiorpoint methods solve efficiently smalltomediumsized SDP problems but do not scale well for largescale problems because of the Hessian equation that implies forming and factorizing a fullydense matrix at each step. Yet, sparse largescale problems can be tackled exploiting sparsity. In this paper, we focus on clique decomposition techniques [5, 15, 20] that perform well for ACOPF problems. Jabr [9] and Molzahn et al. [14] have each proposed a way to use clique decompositions for ACOPF problems and both propose to work on problems formulated in complex variables. However, there are many other possible clique decomposition procedures that provide different reformulations of a SDP problem. Most of the clique decomposition algorithms seek to minimize the number of added edges in the chordal extension. It results in decompositions with many small cliques for ACOPF problems and many linking constraints are required to handle overlaps between cliques. As linking constraints can slow down the resolution, Nakata et al. [15] and Molzahn et al. [14] propose heuristics for merging cliques in order to reduce their number. Yet, merging cliques means adding more edges in the chordal extension. Therefore, classical clique decomposition approaches may not always be the most appropriate since they focus on the number of added edges, regardless of important criteria such as the size of the largest clique or the number of linking constraints.
In this paper, we show that different clique decompositions are not equivalent in terms of resolution by comparing different chordal extension algorithms on RTE [10] and other MATPOWER [22] datasets. Our main contribution is to demonstrate in two ways that computing a chordal extension with a minimum number of added edges is not a relevant choice for ACOPF problems. The first test consists in a comparison between reformulations coming from the SDP problem formulated in complex variables and reformulations coming from the same problem but formulated in real variables. The conclusions of this comparison are that less edges are necessary in the real case but the reformulations are more performant in the complex case. The second test relies on a new clique combination algorithm that improves a given decomposition by adding more edges in the chordal extension. The edges are added in such a way as to minimize the number of linking constraints while keeping small cliques. Computational tests show that this algorithm significantly speeds up the resolution of the largest OPF instances. Both tests highlight the fact that the number of added edges is not always the best criterion to minimize, which proves that the computation of the chordal extension is a question that deserves to be investigated.
This paper is organized as follows. Section 2 describes briefly the OPF problem and the classical rank relaxation. Common clique decomposition techniques are presented in section 3 along with computational comparisons. Section 4 compares reformulations coming from the complex SDP problem and reformulations coming from the real SDP problem. Section 5 details our clique combination algorithm and resulting enhancements. Section 6 concludes the paper.
2 ACOPF problem and rank relaxation
In this section, we present briefly the ACOPF formulation and its classical semidefinite rank relaxation.
2.1 ACOPF formulation
The ACOPF problem is defined on a power transmission network that can be modelled as an oriented graph where represents the electrical buses and the branches (transmission lines, transformers). Let us denote the subset of generator buses. The ACOPF problem is defined as follows:
min  (ACOPF)  
s.t.  
where stands for the real part of the complex number and for the imaginary part. All constant parameters are in bold. and represent respectively the linear cost and the constant cost of the generator bus . represents the load. and are bounds on the voltage magnitude at bus . and (respectively and ) are bounds on the active (respectively reactive) power at generator bus . represents the current limit for branch . The variables are the voltage and the power at each bus with for all buses . Both are complex variables.
is the set of entering branches at bus and the set of exiting branches. stands for the origin of branch , for the destination of branch . The currents and of a branch are linear functions of and depending on physical characteristics. The power on the lines are defined as follows:
(1) 
Note that we use thermal limits modelled with current and linear generation costs to have a formulation closer to what is used in practice at RTE.
2.2 Rank relaxation
Since costs are linear, the ACOPF problem can be written in terms of voltage variables only. The resulting problem is a Quadratically Constrained Quadratic Program (QCQP) with complex variables that can be expressed in a compact way:
(2) 
The classical rank relaxation is constructed by introducing the Hermitian matrix . This constraint is equivalent to and . The rank relaxation is the SDP problem obtained by suppressing the rank constraint:
(3) 
This SDP problem with complex numbers can be converted into a SDP problem with real numbers using the rectangular representation and a symmetric matrix of size .
3 Clique decomposition analysis
In this section, we review the clique decomposition techniques based on matrix completion introduced by Fukuda et al. in [5, 15]. Then we show that the resolution with interiorpoint methods is highly sensitive to the clique decomposition.
3.1 Clique decomposition framework
The clique decomposition relies on the matrix completion theorem [20]. This allows to replace the Positive Semidefinite (PSD) constraint on a big matrix by several PSD constraints on smaller submatrices at the price of adding linking constraints between these submatrices.
Let us define as the aggregate sparsity pattern of a SDP problem, i.e., the matrix of nonzeros entries in the data matrices (objective and constraint matrices). Let us denote as the graph associated to A. The matrix completion theorem can be applied if and only if is chordal. A chordal graph is a graph with no induced cycle of four vertices or more. A clique is a subset of vertices that are all connected together. A maximal clique is a clique that is not included in another clique.
The general framework for clique decomposition contains three steps. The first step consists in computing a chordal extension H of the aggregate sparsity pattern . The second step consists in determining the list of maximal cliques in H, denoted by . These maximal cliques define the submatrices that must be positive semidefinite in the reformulated SDP problem. From the list of maximal cliques , a clique graph W can be defined with nodes corresponding to maximal cliques and weighted edges between each pair of vertices defined by the number of shared nodes between the two corresponding cliques. The third step consists in computing a clique tree U by computing a maximumweight spanning tree for W. This clique tree allows to specify linking constraints between submatrices.
It is clear that the clique decomposition depends on the chordal extension H but it is not completely clear how to precisely define a “good” decomposition from a practical point of view. One reasonable choice is to use decompositions that minimize the number of additional edges. However, finding the minimal chordal extension is NPcomplete [21]. Bergman et al. have recently proposed an exact model with exponentially many constraints to find a minimal extension [2] but it is only usable for small instances. For this reason, several efficient heuristics for adding few edges are available in the literature [8]. The most commonly used is based on a Cholesky factorization of whose rows and columns have been permuted according to a minimum or approximate minimum degree ordering [1]. The ordering has a significant impact on the fillin, that is, on the number of edges added to . Once H is computed, the list of maximal cliques can be computed in linear time [19]. Finally, the most widely used exact algorithm to compute the clique tree is Prim’s algorithm [16] but there are other algorithms and they can lead to different optimal clique trees. However, the decompositions obtained usually do not differ much: the number of linking constraints do not differ significantly from one tree to another because the trees have the same overall weight.
We show in the following subsection that different chordal extension algorithms give different decompositions and that these different decompositions are not equivalent as regards resolution time.
3.2 Numerical comparison of two chordal extension heuristics for SDP relaxations of ACOPF formulated in complex variables
There are several algorithms to compute chordal extensions and the clique decompositions can vary greatly from one algorithm to another, which has an impact on the resolution. In this subsection, we compare clique decompositions coming from two different algorithms: a Cholesky factorization with an AMD ordering and the Minimum Degree (MD) heuristic [3] that is based on a dynamical computation of a minimum degree ordering of the vertices. We have picked these two algorithms because they are among the best heuristics known. Both algorithms are applied on the SDP relaxations of ACOPF formulated in complex variables, therefore providing maximal cliques with complex variables. In order to solve the SDP problems, these cliques are finally converted into cliques with real variables by doubling them.
Table 1 presents some results of this comparison on MATPOWER instances with more than 1000 buses. The tests were carried out on a Processor Intel® Core™ i76820HQ CPU @2.70GHz using the module MathProgComplex.jl [18] with Julia 1.0.3. and the SDP solver MOSEK 8.1.0.72 was used for the resolution. We also used the packages JuMP.jl [4], LightGraphs.jl [17] and MetaGraphs.jl.
Let , and be respectively the number of maximal cliques, the number of linking constraints and MOSEK resolution time obtained with the MD heuristic. Similarly, let , and be respectively the number of maximal cliques, the number of linking constraints and MOSEK resolution time obtained with Cholesky and AMD. Table 1 presents some ratios between these three quantities. All results are presented in percentages.
Instance  

case1888rte  0.17%  1.7%  6.1% 
case3012wp  0.14%  3.3%  15.1% 
case6468rte  0.10%  0.31%  6.8% 
case9241pegase  0.02%  4.9%  62% 
case13659pegase  0.04%  4.7%  48% 
This table is not exhaustive but some conclusions can be drawn from it. First, clique decompositions are different from one algorithm to another. In particular, there are large differences in resolution time even if the decompositions have about the same number of cliques. Theses differences can be explained by the differences in the number of linking constraints, especially for the biggest instances. Another factor can be the size of the biggest clique in the decomposition. It is also possible to generate decompositions which differ both in the number of cliques and in the number of linking constraints so differences in resolution time can be much more impressive. In particular, in preliminary tests we experienced that there exists “bad” clique decompositions by computing several “bad” orderings for Cholesky (e.g. random ordering, maximum degree ordering). In this case, a “bad” clique decomposition means a decomposition with many linking constraints, which can lead to memory issues with interiorpoint methods. This confirms the idea that the computation of the chordal extension is a point worth studying.
The results shown in this section are a good motivation to continue exploring the impact of the clique reformulation on the resolution. The next section focus on the impact of computing the clique decomposition on SDP relaxations of ACOPF problems formulated in complex variables or on SDP relaxations of ACOPF problems formulated in real variables.
4 Comparison of clique decompositions computed in the complex and in the real framework
SDP relaxations of ACOPF problems are naturally expressed with complex numbers but in practice, in order to be solved by the available solvers, they are rewritted in terms of only real variables. The process of rewriting the model with real variables leads to a different (and bigger) aggregate sparsity pattern A. Such a A could be itself decomposed in a way that would not be possible with the original formulation with the complex variables. Therefore, there are two possibilities for the clique decomposition: either applying the procedure on the complex SDP formulation and converting the complex cliques to real cliques or directly applying the clique decomposition procedure on the real SDP formulation. However, computing the chordal extension on the complex or the real problem is not equivalent. We first show it theoretically on a small MATPOWER [22] example, LMBM3.
Let us define (respectively ) as the graph associated to the aggregate sparsity pattern A of the complex (respectively real) SDP relaxation. For any ACOPF instance as defined in section 2, is the network graph . Let us denote (respectively ) the chordal extension computed from (respectively ) with a given algorithm. The chordal extension can be converted to real numbers in the same way as the complex graph can be converted to the real graph . Let us denote this conversion of . is also a chordal extension of . The objective of this section is to compare and .
For LMBM3, is a clique of size 3, denoted by and is the graph with the 6 following nodes and the 12 following edges:
Since is complete and thus chordal, regardless of the chordal extension algorithm used. Clearly, there is one maximal clique with 3 nodes and no linking constraints. When converted to real numbers, we get , the clique of size 6, which accounts to adding 3 edges to . This chordal extension gives one maximal clique with 6 real nodes and no linking constraints.
On the other hand, is not chordal because there are several induced subgraphs of size 4, e.g. . It suffices to add the two edges and to to get a minimal chordal extension
This small example proves that it is not equivalent to compute the clique decomposition on the complex problem or on the real problem from the theoretical point of view: more edges are added in the complex case to get a chordal extension. This idea is confirmed by computational tests that show that there are on average 50% more edges added when the chordal extension is computed in the complex case rather than in the real case for MATPOWER instances with more than 1000 buses. However, adding less edges in the chordal extension is not necessarily advantageous from the numerical resolution point of view. Indeed, adding less edges in the chordal extension can mean having more cliques and more linking constraints as shown in Table 2. This table presents comparisons between clique decomposition coming from the complex case (computed on ) and clique decomposition coming from the real case (computed on ). More precisely, both clique decompositions are based on a chordal extension obtained by a Cholesly factorization with an Approximate Minimum Degree (AMD) [1] ordering. Let (respectively ) be the list of maximal cliques obtained from the chordal extension (respectively ). The second column represents the ratio for the number of cliques defined as . Similarly, the third column represents the ratio for the number of linking constraints, i.e., with (respectively ) the number of linking constraints when the clique decomposition is done for the complex (respectively real) problem. Finally, the fourth column represents the ratio for MOSEK resolution time, i.e., with (respectively ) the resolution time with the clique decomposition coming from the complex (respectively real) problem. All results are presented in percentages and the tests were carried out in the same manner as in subsection 3.2.
Instance  

case1888rte  51%  24%  17% 
case3012wp  39%  15%  4% 
case6468rte  40%  19%  58% 
case9241pegase  48%  18%  28% 
case13659pegase  59%  20%  23% 
This table shows that the clique decomposition procedure on gives approximatly 50% more cliques and 20% more linking constraints than on , which slows down the resolution (the computing time increases by at least 4% more and up to 58% more). Therefore, even if less edges are added to the graph when working directly in real numbers, it seems better to work in complex numbers. This comparison demonstrates that computing a chordal extension with a minimum number of added edges is not necessarily profitable. The next section presents a clique combination algorithm that leads to the same conclusion.
5 A new clique combination approach
Several heuristics do already exist to compute chordal extensions but to the best of our knowledge, they all try to minimize the number of added edges. However, nothing proves that a chordal extension with few added edges gives a decomposition easier to solve. In contrast, minimizing the size of the biggest clique or the number of linking constraints could be more interesting because these quantities have an impact on the resolution. To support this claim, we propose a new clique combination algorithm that allows to reduce the number of linking constraints in a given decomposition while keeping small matrices. This algorithm consists in merging some cliques that are adjacent in the clique tree by solving the problem (IP) below. There is one binay variable per edge in the clique tree. If this variable equals 1, it means that the cliques at the ends of the edge will be merged. We allow at most one combination per clique and combinations that result in a clique of size greater than a certain size are forbidden. The objective is to maximize the number of linking constraints that will be eliminated thanks to the combinations. More precisely, the model is the following:
(IP) 
with

the clique tree computed for a given decomposition ;

the number of common nodes for ;

the size of the merged clique for ;

a parameter that defines the maximal size of merged cliques;

the set of edges that have as extremity.
Number of times (IP) is applied  0  1  2  3  4  5  6 
Total resolution time (seconds)  3397  2316  1807  1691  1533  1710  1826 
This algorithm has been tested with the solver Xpress for and the results are presented in Table 3. This table represents the evolution of the total resolution time depending on the number of times the combination algorithme (IP) is used. This table shows that applying this algorithm at least once improves significatively MOSEK resolution time and the best time is obtained when the algorithm is applied 4 times. Therefore, it is worth combining cliques, i.e., adding more edges in the chordal extension, to improve the performance of a given decomposition for ACOPF problems.
6 Conclusion and future work
Sections 4 and 5 show that, for ACOPF problems, there is no real reason to build chordal extensions with as few added edges as possible since adding more edges can lead to better performance. These results encourage us to change the approach for computing chordal extensions and in the future we will focus on the formulation of an adapted heuristic that computes performant chordal extensions for ACOPF problems.
Footnotes
 There are other possibilities for , e.g., if the edges and are added, but the results will be the same.
References
 Patrick R Amestoy, Timothy A Davis, and Iain S Duff. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications, 17(4):886–905, 1996.
 David Bergman, Carlos H Cardonha, Andre A Cire, and Arvind U Raghunathan. On the minimum chordal completion polytope. Operations Research, 67(2):532–547, 2019.
 Anne Berry, Pinar Heggernes, and Genevieve Simonet. The minimum degree heuristic and the minimal triangulation process. In International Workshop on GraphTheoretic Concepts in Computer Science, pages 58–70. Springer, 2003.
 Iain Dunning, Joey Huchette, and Miles Lubin. Jump: A modeling language for mathematical optimization. SIAM Review, 59(2):295–320, 2017.
 Mituhiro Fukuda, Masakazu Kojima, Kazuo Murota, and Kazuhide Nakata. Exploiting sparsity in semidefinite programming via matrix completion i: General framework. SIAM Journal on Optimization, 11(3):647–674, 2001.
 H. Godard, S. Elloumi, A. Lambert, J. Maeght, and M. Ruiz. Novel approach towards global optimality of optimal power flow using quadratic convex optimization. In 2019 6th International Conference on Control, Decision and Information Technologies (CoDIT), pages 1227–1232, April 2019.
 Ajit Gopalakrishnan, Arvind U Raghunathan, Daniel Nikovski, and Lorenz T Biegler. Global optimization of optimal power flow using a branch & bound algorithm. In 2012 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 609–616. IEEE, 2012.
 Pinar Heggernes. Minimal triangulations of graphs: A survey. Discrete Mathematics, 306(3):297–317, 2006.
 Rabih A Jabr. Exploiting sparsity in sdp relaxations of the opf problem. IEEE Transactions on Power Systems, 27(2):1138–1139, 2012.
 Cédric Josz, Stéphane Fliscounakis, Jean Maeght, and Patrick Panciatici. Ac power flow data in matpower and qcqp format: itesla, rte snapshots, and pegase. arXiv preprint arXiv:1603.01533, 2016.
 Cédric Josz, Jean Maeght, Patrick Panciatici, and Jean Charles Gilbert. Application of the momentsos approach to global optimization of the opf problem. IEEE Transactions on Power Systems, 30(1):463–470, 2015.
 Javad Lavaei and Steven H. Low. Zero duality gap in optimal power flow problem. IEEE Transactions on Power Systems, 27(1):92–107, February 2012.
 D. K. Molzahn and I. A. Hiskens. A Survey of Relaxations and Approximations of the Power Flow Equations. Foundations and Trends in Electric Energy Systems, 4(12):1–221, February 2019.
 Daniel K Molzahn, Jesse T Holzer, Bernard C Lesieutre, and Christopher L DeMarco. Implementation of a largescale optimal power flow solver based on semidefinite programming. IEEE Transactions on Power Systems, 28(4):3987–3998, 2013.
 Kazuhide Nakata, Katsuki Fujisawa, Mituhiro Fukuda, Masakazu Kojima, and Kazuo Murota. Exploiting sparsity in semidefinite programming via matrix completion ii: Implementation and numerical results. Mathematical Programming, 95(2):303–327, 2003.
 Robert Clay Prim. Shortest connection networks and some generalizations. The Bell System Technical Journal, 36(6):1389–1401, 1957.
 James Fairbanks Seth Bromberger and other contributors. Juliagraphs/lightgraphs.jl: an optimized graphs package for the julia programming language, 2017.
 J. Sliwak, M. Ruiz, M. F. Anjos, L. Létocart, and E. Traversi. A julia module for polynomial optimization with complex variables applied to optimal power flow. In 2019 IEEE Milan PowerTech, pages 1–6, June 2019.
 Robert E Tarjan and Mihalis Yannakakis. Simple lineartime algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs. SIAM Journal on computing, 13(3):566–579, 1984.
 Lieven Vandenberghe, Martin S Andersen, et al. Chordal graphs and semidefinite optimization. Foundations and Trends in Optimization, 1(4):241–433, 2015.
 Mihalis Yannakakis. Computing the minimum fillin is npcomplete. SIAM Journal on Algebraic Discrete Methods, 2(1):77–79, 1981.
 Ray Daniel Zimmerman, Carlos Edmundo MurilloSánchez, Robert John Thomas, et al. Matpower: Steadystate operations, planning, and analysis tools for power systems research and education. IEEE Transactions on power systems, 26(1):12–19, 2011.