Virtual-Voltage Partition-Based Approach to Mixed-Integer Optimal Power Flow Problems

Virtual-Voltage Partition-Based Approach to Mixed-Integer Optimal Power Flow Problems

Chin-Yao Chang  Sonia Martínez  Jorge Cortés C.-Y. Chang, S. Martínez, and J. Cortés are with the Department of Mechanical and Aerospace Engineering, UC San Diego. A preliminary version of this work appeared as [1] at the 2017 Allerton Conference on Communications, Control, and Computing.

This paper deals with optimal power flow (OPF) problems with discrete variables that capture binary decisions about network topology configurations and capacitor bank settings. We adopt a semidefinite programming formulation for the OPF problem which, however, remains nonconvex due to the presence of discrete variables and bilinear products between the decision variables. To tackle the latter, we introduce a novel physically-inspired, virtual-voltage approximation that leads to provable lower and upper bounds on the solution of the original problem. To deal with the exponential complexity caused by the discrete variables, we introduce a graph partition-based algorithm which breaks the problem into several parallel mixed-integer subproblems of smaller size. Simulations on the IEEE 30, 118, 300 bus test cases demonstrate the high degree of accuracy and affordable computational requirements of our approach.

I Introduction

Mixed integer programming (MIP) arises in optimal power flow (OPF) problems involving switching control decisions, such as transformer tap settings, transmission line connections, and capacitor bank switching. The optimization of OPF problems with continuous and discrete variables has the potential to yield significant benefits in efficiency and reliability, but the complexity of solving such highly non-convex problems makes achieving this goal difficult. System operators usually select switching control laws offline to handle contingency situations and maintain system stability. In this way, they decouple the decisions about discrete variables from the OPF problem: once discrete controls have been specified, OPF focuses on generation cost or loss minimization. These observations motivate our focus on the co-optimization of discrete decisions and OPF.

Literature review: The works [2, 3] provide surveys on the solution of OPF problems with discrete variables. Many of the methods employed to solve OPF problems have been extended to deal with MIP-OPF problems, e.g., particle swarm optimization [4, 5] and genetic algorithms [6]. Transmission line switching or network topology reconfiguration commonly serves as a corrective mechanism in response to system contingencies see [7, 8] and references therein. In [9, 10], linearized OPF, also known as DCOPF, is used for the fast co-optimization of network topology and OPF. Despite its relative low complexity, DCOPF may lead, especially in congested systems, to poor solutions that can even result in voltage collapse [11, 12]. The work [13] proposes quadratic convex (QC) relaxations for the MIP-OPF problem, which provides more accurate results than DCOPF, while still retaining a fast computation time. Recent studies [13, 14] show that methods based on semidefinite programming (SDP) convex relaxations of ACOPF may lead to better solutions than DCOPF and QC. However, how to handle discrete variables in the context of SDP is challenging and not fully understood. The challenges stem from not only from the integer-valued nature of these variables, but also from the presence of bilinear terms involving the product of discrete decision variables with continuous ones, reflecting the impact on the physical modeling of the line being connected. The paper [14] uses a lift-and-branch-and-bound procedure to deal with the SDP formulation of MIP-OPF problems, but has still exponential complexity in the worst case due to the nature of the branch-and-bound procedure. The work [15] also uses SDP to solve MIP-OPF problems, where bilinear terms associated to line connections are addressed by assuming certain nominal network topology and bilinear terms of other discrete decision variables are dealt with using the McCormick relaxation [16].

Statement of contributions: We consider the co-optimization of transmission line switching and optimal power flow. We introduce the SDP convex relaxation of the OPF problem and add discrete variables that can model transmission line and capacitor bank switching. Our contributions are twofold. First, we introduce a novel way of dealing with the nonconvexity coming from the presence of bilinear terms, which we term virtual-voltage approximation. Our approach is based on introducing virtual voltage variables for the terminal nodes of each switchable line and impose physically-meaningful constraints on them. We show that this approach leads to sound bounds on the optimal power losses and individual flows of the switchable lines, and provides lower and upper bounds on the optimal value of the original problem. Our second contribution deals with the nonconvexity coming from the discrete variables. We build on the virtual-voltage approximation to propose a graph partition-based algorithm that significantly reduces the computational complexity of solving the original problem. This algorithm uses the values of the optimal dual variables from the virtual voltage method to define a weighted network graph, which is then partitioned with a minimum weight edge-cut set. The algorithm breaks the original network into sub-networks so as to minimize the correlation between the solutions to the optimization problem on each sub-network. Finally, the algorithm solves a mixed-integer SDP problem on each sub-network in parallel and combines them to reconstruct the solution of the original problem. We implement the proposed algorithms on IEEE 30, 118, and 300 bus test cases, and compare them with available approaches from the literature to illustrate their superior performance regarding convergence to the optima and computation time.

Ii Preliminaries

This section introduces basic concepts used in the paper111We use the following notation. We denote by , , , and the sets of positive integer, real, positive real, and complex numbers, resp. We denote by the cardinality of . For a complex number , we let and be the complex modulus and angle of , and its real and imaginary parts are and . We let denote the -norm of . Let and be the set of positive semidefinite and -dimensional Hermitian matrices, resp. For , we let and denote its conjugate transpose and trace, resp. We let denote the element of ..

Ii-1 Graph Theory

We review basic notions of graph theory following [17]. A graph is a pair , where is the set of nodes and is the set of edges. A self-loop is an edge that connects a node to itself. The graph is undirected if . A path is a sequence of nodes such that any two consecutive nodes correspond to an edge. The graph is connected if there exists a path between any two nodes. An orientation of an undirected graph is an assignment of exactly one direction to each of its edges. A graph is simple if it does not have self-loops or multiple edges connecting any pair of nodes. Throughout the paper, we limit our discussion to undirected, simple graphs. A vertex-induced subgraph of , written , satisfies and . An edge cut set is a subset of edges which, if removed, disconnects the graph. A weighted graph is a graph where each branch has a weight, . Given the edge weights , the adjacency matrix has if , and otherwise. The degree matrix is a diagonal matrix such that . The Laplacian matrix is . The Laplacian matrix is positive semidefinite, with zero as an eigenvalue and multiplicity equal to the number of connected components in the graph. The Fiedler vector is the eigenvector associated with the second smallest eigenvalue of . An -partition of a connected divides into a number of connected vertex-induced subgraphs, , such that and for all . An -optimal partition of is an -partition of with minimized, where . Spectral graph partitioning partitions a connected graph into two vertex-induced subgraphs, and , where and are the nodes corresponding to the positive and non-positive entries of the Fiedler vector, respectively.

Ii-2 McCormick Relaxation of Bilinear Terms

The McCormick envelopes [16] provide linear relaxations for optimization problems that involve bilinear terms. Consider a bilinear term on the variables , for which there exist upper and lower bounds, . The McCormick relaxation consists of substituting in the optimization problem the term by its surrogate and adding the following McCormick envelopes on ,


Constraints (1) are tight, in the sense that each plane in (1) is tangent to the bilinear-constraint manifold at two boundary lines. The convex polyhedron in the variables encloses the actual bilinear-constraint manifold.

Iii Problem Statement

We begin with the formulation of the OPF problem over an electrical network and its SDP convex relaxation following [18]. Then, we introduce binary variables for the co-optimization of line switching and capacitor banks, leading to the problem formulation of interest in this paper.

Consider an electrical network with generation buses , load buses , and electrical interconnections described by an undirected edge set . Let and denote its cardinality by . We denote the phasor voltage at bus by , where and are the voltage magnitude and phase angle, respectively. For convenience, denotes the collection of voltages at all buses. The active and reactive power injections at bus are given by the power flow equations


where are the active and reactive power demands at bus , and are derived from the admittance matrix Y as follows


Here denotes the canonical basis of . The OPF problem also involves the box constraints


where is the upper bound of the voltage difference between buses , and and are the lower and upper bounds of the voltage magnitude at bus , respectively. All , are defined similarly. The objective function is typically given as a quadratic function of the active power,


where , and . The OPF problem is the minimization over (5) subject to (2) and (III). Such optimization is non-convex due to the quadratic terms on . To address this, one can equivalently define (or and ) as the decision variable (note all the terms in (2), (III) and (5) are quadratic in ). Dropping the rank constraint on makes the OPF problem convex, giving rise to the SDP convex relaxation,

subject to


where are defined so that and .

We are interested in solving the OPF problem with discrete variables that model topology reconfigurations via line switching. Assume that an additional set of edges, , could be added to the network . Topological changes could make the system more robust and reduce generation cost. For each line , we define a binary variable , and we say the line is connected if and disconnected otherwise. If , then the power flow from node to through edge is


where are defined as follows: all entries are prescribed to be zero except the ones defined by222We omit charging susceptance, tap ratio and phase shift of transformers for simplicity.

Here is the admittance of line . Taking (7) into account, the active and reactive power of each node become


where . Given , we use (P1)- to refer to the OPF problem solved with the network topology with extra lines as determined by .

We are interested in solving what we call (P2), which is the optimization (P1) with constraints (6a) and (6b) replaced by (7) and (8). In addition to optimizing the power flows, this corresponds to also selecting the optimal set of lines to connect among the ones in . The problem (P2) is non-convex for two reasons: the binary variables and the bilinear products of and . The first problem can be addressed using existing integer-programming solvers [19, 15]. The McCormick relaxation described in Section II-2 is the standard way to deal with the second problem. In this paper, we instead provide alternative routes to address each of these problems for the optimization (P2).

Remark III.1.

(OPF with capacitor bank switching). The formulation above focuses on the co-optimization of OPF and line switching, but it can also accommodate capacitor-bank switching. This generalization entails an equation similar to (8) to adjust the reactive power injection as follows


where , and can take multiple positive integer values. The latter can be equivalently represented via binary variables. The last term in (9) entails non-convex bilinear products of decision variables, as described above.

Iv Virtual-Voltage Approximation of Bilinear Terms

We introduce here a novel way to deal with the bilinear terms in (P2) which we term virtual-voltage approximation. We start by noting that every binary variable only appears in the bilinear products in (8) together with another continuous variable . If we convexify the binary variables by having them take values in , then we can interpret each bilinear term corresponding to as a line power flow from to , with the magnitude bounded by what indicates. Following this reasoning, if the direction of power flow of every line was known, then the bilinear term would no longer be an issue. For example, if we knew that and , then we could define new variables, and , replacing and in (8), respectively, and impose


This would eliminate the bilinear terms and the only remaining non-convexity would be that the physical feasible solution should satisfy and . In general, however, the direction of power flow of the lines is not known a priori and, hence, the trivial convex constraints (10) for the relaxation are no longer valid.

Iv-a Convex Relaxation Via Virtual Voltages

Our idea to approximate each bilinear term builds on defining a virtual voltage for the terminal nodes of the line and impose constraints on them to make sure they have physical sense. We make this precise next. Let be an arbitrary orientation of . To define the virtual voltages, and in keeping with the SDP approach, for each we introduce a two-by-two positive semidefinite matrix . This matrix encodes physically meaningful voltages at the terminal nodes if its rank is one, namely, , with and corresponding to the voltages of nodes  and , respectively. For convenience, we introduce and impose the following constraints on 


Constraints (11a) and (11b) ensure that the voltage magnitudes of and derived from are no bigger than the ones from . Constraint (11c) ensures that the voltage difference between nodes and computed from is less than the corresponding difference from . Therefore, if the matrix has rank one, constraints (11) ensure that we obtain physically meaningful and feasible voltage values.

Let be the principal sub-matrix of by only keeping the rows and columns associated with nodes  and . We define similarly. We replace and in (8) by and , respectively. We now have all the elements necessary to convexify (P2) as follows

subject to (6c)-(6e), (11), and

Remark IV.1.

(OPF with capacitor bank switching – continued). The formulation (P3) can also be employed to relax the bilinear terms involved in capacitor banks, cf. Remark III.1. Let be the binary control of a capacitor bank. We replace in (9) by a new variable and impose the following constraints on 

where is a two-by-two real-valued positive semidefinite matrix. In this way, if , and when the inequalities become equalities.

Each optimal solution of (P3) has a dominant eigenvalue, much larger than the other one. To formally state the result, let denote the principal sub-matrix of the optimum of (P3) obtained by removing from all columns and rows except the ones corresponding to and . We use the spectral decomposition to rewrite  as

where , , and is the condition number of . Lemma IV.2 establishes a useful lower bound on .

Lemma IV.2.

(Lower bound on the condition number). For all , the optima of (P3) has


By constraint (11c) and (6e), we have

Lemma IV.2 follows by rearranging the inequality above. ∎

For all practical purposes, the result of Lemma IV.2 implies that the matrix specifies well-defined virtual voltages at the terminal nodes, as we explain next.

Remark IV.3.

(Optimal solutions have well-defined virtual voltages). Using Lemma IV.2, we justify that the optimal solution has a dominant eigenvalue as follows. The denominator of (13) is always non-negative due to (11c). The order of the denominator of (13) is at most as in most test cases. On the other hand, when the virtual voltage satisfies or , then the numerator is lower bounded by a scalar close to one, as . As a consequence, the fraction in (13) is usually bigger than . Our simulations on IEEE 118 and 300 bus text cases confirm that is at least .

Iv-B Physical Properties of the Convex Relaxation

The active and reactive power flows in (P3) on a switchable line are determined by  according to


The next result shows that the optimal power losses on each edge are bounded by the ones computed from .

Lemma IV.4.

(Bounds on the sums of line active and reactive powers). If the line charging susceptance is zero for all , then the following inequalities hold


If the line charging susceptance is zero, then and take the following form

Since both and are non-negative, the result of Lemma IV.4 follows from (11c) and the equalities (14). ∎

We next seek to upper bound the individual flows and . Our next result shows that, under certain conditions for (P3), stronger properties hold on the active power retrieved from the optimal solution and .

Proposition IV.5.

(Bounds on directional power flow). Let for each . Assume is purely inductive and has zero charging susceptance, and


Then the following inequalities hold


If is purely inductive and has zero charging susceptance, then

Note that only the off-diagonal entries of and are non-zero, making and . If , (17) follows as . It is then enough to show that if , . We show it by contradiction. If , then and , where , and is the angle difference between and . Using (11a), we define such that , and rewrite the inequality as . Then,


Rewriting (11c) as a function of ,


where we use in the inequality. Using (18), the RHS of (19) is less than


The derivative (20) with respect to is


The first two elements summed up to a non-positive value due to (16). We then conclude that (21) is non-positive with given and fixed. Therefore, (20) is non-positive for all because it is zero when and is non-increasing. But (20) is strictly larger than the RHS of (19), contradicting (19). ∎

Condition (16) holds for most existing power systems [20]. An analogous result holds by restricting  instead.

Proposition IV.6.

(Bounds on directional power flow. II). If is purely inductive and has zero charging susceptance, and is sufficiently small such that (16) holds, then (17) follows.

The proof of Proposition IV.6 is analogous to that of Proposition IV.5 and therefore we omit it. Similar bounds as (17) follow for reactive power if the sum of the cosine terms in the bracket of (19) is non-negative (however, in general, this is not the case). In addition, more involved, inequalities as (17) hold for the general impedance case, but we do not pursue them here. Propositions IV.5 and IV.6 show that, when the diagonal entries of are at the boundary points of their constraints, (P3) eliminates the bilinear terms on the active line power flow of (P2) in the same way as (10). The difference between the relaxations is that there is no need in (P3) to know the direction of line power flow a priori, as opposed to (10).

Iv-C Reconstructed Solution to the MIP-OPF Problem

We note that the ratio of the voltage magnitudes derived from and provides an approximation of the discrete variables in (8) as


Note that because of (11). If we round the entries of to the closest number in , we obtain a candidate solution to (P2). The following result, whose proof is straightforward, states the relationship between (P2) and (P3) based on the rounded solution .

Proposition IV.7.

(Properties of the reconstructed solution). The optimal values of (P1)-, (P2), and (P3) satisfy . Moreover, if , then the optimal solution of (P1)-, , combined with , is an optimal solution of (P2).

Note that even if , does not necessary equal . The reason is that (22) computes  from the diagonal of , and hence we can not conclude any equality for the off-diagonal elements of and . Hence, even if , the optimal solution of (P3) does not necessarily lie in the feasible region of (P2).

Remark IV.8.

(Comparison with the McCormick relaxation). We explain how we implement the McCormick relaxation on the problem (P2) for comparison purposes. For each , we define new variables to substitute the bilinear terms , respectively. Then, we impose constraints of the form (1) on the new variables based on and upper and lower bounds of active/reactive line power flow, , . If these bounds are far from the actual optimal line power flows, this can significantly affect the quality of the solution obtained by the McCormick relaxation, a point that we illustrate later in our simulations, along with rationale for how to select them. In contrast, the proposed relaxation (P3) is not sensitive to those line power bounds, as the virtual voltages are bounded by the power computed from . Additionally, the variables and in the McCormick relaxation are loosely tied to the decision variable , whereas (P3) introduces constraints (11a)-(11c) enforcing a stronger physical connection between the virtual voltages and .

V Partition-Based MIP-OPF Algorithm

In general, the solution to the convex relaxation (P3) of the MIP-OPF problem introduced in Section IV provides a lower bound on its optimal value, cf. Proposition IV.7. Here, we discuss how to refine the reconstructed solution to better approximate the solution of the original optimization problem. One approach consists of using the branch-and-bound algorithm [21] and relying on (P3) to generate the required branch lower bounds in its execution. However, this approach can easily become intractable as the number of switchable lines grows because of the large number of cases where (P3) must be executed. To address the exponential complexity, we propose an alternative method, described as follows:

[Informal description]: The algorithm partitions the network graph into smaller components to break the original problem into subproblems of smaller size. To do this, we view the disconnection of a line as a perturbation on the constraint (6c). This is natural, in the sense that the disconnection changes the active and reactive power injections of the terminal nodes, which in turn may cause the constraints to be violated. Hence, we partition the graph so that the correlation between the solutions to the optimization problem on each of the resulting subgraphs is minimized. The idea is that, if the optimal solution in each subgraph minimally violates the constraints that connect them to other graphs’ solutions, then, when put together, the reconstructed solution to (P2) would be of better quality.

Algorithm 1 formally presents the partition-based MIP-OPF algorithm.

1:Compute the optimal solution of (P3)
2:Construct graph reduction (Section V-1)
3:Assign adjacency matrix to (Section V-2)
4:Compute cut set to partition into subgraphs (Section V-2)
5:Solve integer optimization problem (P4) on each subgraph to find (Section V-3)
6:Solve (P1)- (Section V-4)
Algorithm 1 Partition-Based MIP-OPF Algorithm

Next, we describe in detail each of its steps.

V-1 Graph reduction

This is a step prior to graph partitioning which is motivated by the following observation. The graph partitioning should not result in nodes connected by a switchable line belonging to different subgraphs. This is because, if that were the case, then solving the OPF associated with each subgraph cannot capture how the switch in that specific line affects the optimal value of the original MIP-OPF problem. To address this, we ‘hide’ the nodes that are connected by  to the partitioning algorithm that finds the edge cut  so as to ensure . Let and let be the set of nodes that are connected to node through a line in . All nodes in are clustered as one representative node and all the edges connected to one of are considered being connected to the representative node. This results in a graph , where is the collection of representative nodes. Notice that and is a strict subset of if there is such that a path connecting nodes and exists in the graph . Figure 1 illustrates the construction of  and has as one edge of is dropped in the process of graph reduction.

Fig. 1: Graph reduction. Nodes connected by  are collapsed into one node. The dash lines denote edges in ; the solid lines denote the edges in .
Remark V.1.

(Networks where all lines are switchable). The graph reduction step described above only makes sense when not all lines are switchable because otherwise, it results in a graph with a single node. Some works [22, 13], however, consider scenarios where all lines are switchable. For such scenarios, one can either skip the graph reduction step, or pre-select a set of lines that should remain non-switchable. We come back to this point in Section VI later.

V-2 Graph partitioning

Our next step is to find an edge cut set of the graph . In order to minimally affect the optimal value , the graph partitioning is based on the optimal dual variables of (P3). The optimum dual variables measure how the optimal value of (P3) changes with respect to the corresponding constraint. Formally, by taking the derivative of the Lagrangian of (P3), we have the following for ,

and for ,

With this interpretation, we define edge weights as follows

Given the adjacency matrix , we do an -optimal partition on , which gives with . Since all the removed edges are in , we can use the same cut for the partition of : with . Such partition ensures . The intuition is that the cut minimally perturbs because it select edges with minimal weight for the weighted graph .

Finding the optimal cut set is NP-hard. There are algorithms [23, 24] that can approximate it in a few seconds for graphs of the order of a thousand nodes. However, they do not guarantee that the resulting subgraphs are connected. To ensure this property, we resort to spectral graph partitioning.

Theorem V.2.

(Fiedler’s theorem of connectivity of spectral graph partitions). The two subgraphs resulting from spectral graph partitioning on a connected graph are also connected.

The proof is available in [25, Corollary 2.9]. To derive a -partition, one can implement spectral graph partitioning recursively  times. Since we aim for subgraphs with similar size, each iteration applies spectral graph partitioning on the subgraph with the largest number of nodes. The most computationally expensive step in this process is the eigenvector computation, which only takes linear time, or . This low complexity is reflected on the negligible computational time in our simulations. Even though this recursive spectral partitioning does not in general lead to an -optimal partition, the sum of the weights of the edge cut set, , can be upper bounded via Cheeger inequalities [26, 27].

V-3 Integer optimization on subgraphs

Given a -partition , we define an optimization problem associated with each subgraph. This problem is a variant of (P2) that is convenient for the reconstruction of the solution of (P2) over the original . For subgraph , let be its set of edges, the decision variable, the set of switchable lines, and the set of nodes in that connects to at least one node of another subgraph. Each subgraph solves

subject to