Widespread adoption of new grid technologies is continuously changing the face of our modern electricity networks. Increased penetration of renewable energy sources and changing load patterns has lead to higher volatility in power networks [1]. The stochastic nature of renewables and active loads is likely to produce recurring disturbances that will require careful planning and design of power networks with stability in mind [1]. Additionally, the loss of rotational inertia in systems with high percentage of renewables will significantly reduce the capability of power grids to handle such disturbances [2].

While several works have explored the placement of virtual inertia to improve the dynamic performance of power networks [3], it is not well understood how grid topology affects transient stability [4]. Recent works have shown that the grid Laplacian matrix eigenvalues allows us to quantify the effect of grid topology on the power network [4]. Our work is further motivated by the claim that grid robustness against low frequency disturbances is determined by the connectivity of the network [4]. Past studies have looked at designing grid topologies for specific goals such as reduction of transient line losses [1], improvement in feedback control [5], and coherence-based network design [6]. Topologies can also be designed to pursue other objectives such as improving reliability, minimizing investment and operating costs, reducing line losses, and managing network congestion; see [7], [8] and references therein.

This paper investigates the effect of topology on power system dynamics. For a variety of $\mathcal{H}_{2}$-norm based stability metrics, such as reduction of line losses, fast damping of oscillations, and network synchronization, our previous works [9], [10] presented methods that can be used to optimally design power grid topologies. In [10], we presented reformulations of the combinatorial topology design task that allowed us to solve the problem to optimality, albeit with forbiddingly slow run-times. This work significantly improves the computational performance of the previous formulation using the following key contributions discussed in Section IV: i) present methods to derive tight bounds on the continuous optimization variables; ii) generate a-priori valid cuts based on graph-theoretic properties; and iii) showcase an eigenvector-based cut-generation procedure that is able to interject the solver and add constraints on-the-fly. To explore conditions under which a greedy scheme would perform well, we also present new conditions for supermodularity in Section V.

Notation: Column vectors (matrices) are denoted by lower (upper) case letters, and sets by calligraphic symbols. The cardinality of set $\mathcal{X}$ is denoted by $|\mathcal{X}|$ and $\emptyset$ is an empty set. Given a real-valued sequence $\{x_{1},\ldots,x_{N}\}$, $x$ is the $N\times 1$ vector obtained by stacking the entries $x_{i}$ and $\operatorname{dg}(\{x_{i}\})$ is the corresponding diagonal matrix. The operator $(\cdot)^{\top}$ stands for transposition. The $N\times N$ identity matrix is represented by $I_{N}$. The canonical vector $e_{i}$ has a $1$ at the $i$-th entry and is $0$ elsewhere. The time derivative of $\theta$ is denoted by $\dot{\theta}:=\frac{\delta\theta}{\delta t}$. $M\succeq 0$ implies that $M$ is positive semi-definite, and $X_{\mathcal{AB}}$ denotes the matrix obtained from $X$ upon sampling the rows and columns indexed respectively by sets $\mathcal{A}$ and $\mathcal{B}$.

An electric power network with $N+1$ nodes can be represented by a graph ${\cal G}=({\cal V},\cal{E})$, where ${\cal V}:=\{1,2,\ldots,N+1\}$ corresponds to nodes and edges in ${\mathcal{E}}$ correspond to undirected lines. Let node $i=1$ be the reference, and collect the remaining nodes in set $\mathcal{V}_{r}:=\mathcal{V}\backslash\{1\}$. The susceptance of line $(i,j)\in\mathcal{E}$ connecting nodes $i$ and $j\in\cal V$ is denoted by $b_{ij}>0$, while its reactance is denoted by $x_{ij}:=b_{ij}^{-1}$ under a lossless line model. Then, $L$ represents the $(N+1)\times(N+1)$ susceptance Laplacian matrix of graph $\cal G$ and is defined as $L_{ij}:=\sum_{(i,j)\in{\mathcal{E}}}b_{ij},~{}\text{if $i=j$};L_{ij}:=-b_{ij},~{}\text{if $(i,j)\in{\mathcal{E}}$};~{}\text{and}~{}0~{}\text{otherwise}$.

We consider a small-signal disturbance setup where each node $i\in\mathcal{V}$ is associated with a synchronous machine’s rotor angle $\theta_{i}$, frequency $\omega_{i}=\dot{\theta_{i}}$, inertia constant $M_{i}$, and damping coefficient $D_{i}$; see [11]. If a node $i$ hosts an ensemble of devices, the parameters $M_{i}$ and $D_{i}$ represent lumped characterizations of the collective behavior of the hosted devices [3]. Without loss of generality, the quantities $(\theta_{i},\omega_{i})$ will henceforth refer to deviations of nodal phase angles and frequencies from their steady-state values. Using these definitions, the state-space representation of the linearized power grid dynamics can be expressed as [11]

where ${M}:=\operatorname{dg}(\{M_{i}\})$ and ${D}:=\operatorname{dg}(\{D_{i}\})$ are the $(N+1)\times(N+1)$ matrices collecting the inertia and damping constants across all nodes; the $(N+1)$-long vectors ${\theta}$, ${\omega}$, and ${u}$ stack respectively the phase angles, frequencies, and power disturbances at each node. The subsequent analysis relies on the ensuing assumption, which has been justified in several existing works [10, 1, 12].

Given the state-space model in (1), our goal is to design power network topologies that minimize the expected steady-state value of a generalized control objective combining angle deviations and frequency excursions [9], [3]

for given non-negative weights $\{w_{ij}\}$ and $\{s_{i}\}$. The weights $w_{ij}$ induce a connected weighted graph $\mathcal{G}_{w}$ that may be different from $\cal G$. Let $W$ be the Laplacian matrix of graph $\mathcal{G}_{w}$ and define $S:=\operatorname{dg}(\{s_{i}\})$. Then, it is easy to see that $f(t)=\|y(t)\|_{2}^{2}$, where

The importance of the generalized control objective in (2) is that by varying $(W,S)$, one can capture different stability metrics and study them under a unified framework [9]. Note that the network coherence metric considered in the numerical tests corresponds to the case of $W=I_{N+1}-\frac{1}{N+1}1_{N+1}1_{N+1}^{\top}$ and $S=0$ [9].

The expected steady-state value of $f(t)$ can be interpreted as the squared ${\cal H}_{2}$-norm of the linear time-invariant (LTI) system described by (1) and (3). This system will be compactly denoted by $H:=(A,B,C)$. Leveraging this link, the generalized control objective can be expressed as

where $Q\succeq 0$ is the observability Gramian matrix of the LTI system $H$, and can be computed as the solution to the Lyapunov equation [13, Ch. 5]

The ensuing sections select network topologies that minimize the stability metric of (4).

Consider a connected graph $\hat{\mathcal{G}}=(\mathcal{V},\hat{\mathcal{E}})$, where $\hat{\mathcal{E}}$ is the set of all candidate lines. The goal is to find a subset of lines ${\mathcal{E}}\subseteq\hat{\mathcal{E}}$, so that the resultant power network minimizes the stability objective in (4). Because adding lines can be costly and utilities have limited budgets, we further impose the constraint that $|\mathcal{E}|\leq K$ with $K\geq N$. By setting $K=N$, a radial topology can be enforced. The topology design task can be now posed as [10]

Under Assumption 1, the objective of (6) can be shown to be proportional to $\operatorname{Tr}(\tilde{W}\tilde{L}^{-1})$, where $\tilde{W}$ and $\tilde{L}$ are the $N\times N$ matrices obtained after removing the first row and column from $W$ and $L$, respectively; see [10, Lemma 1]. Problem (6) can be then rewritten as [10]

where the rank constraint ensures that $\mathcal{E}$ is connected.

To express the optimization in (7) over $\hat{\mathcal{E}}$ in a more convenient form, let us associate each line $\ell\in\hat{\mathcal{E}}$ with a binary variable $z_{\ell}$, which is $z_{\ell}=1$ if line $\ell$ is selected (that is $\ell\in\mathcal{E}$); and $z_{\ell}=0$, otherwise. If we collect variables $\{z_{\ell}\}_{\ell\in\hat{\mathcal{E}}}$ in vector $z$, then $z$ must lie in the set $\mathcal{Z}:=\left\{z:z^{\top}{1}_{|\hat{\mathcal{E}}|}\leq K,\ z\in\{0,1\}^{|\hat{\mathcal{E}}|}\right\}.$ Based on the line selection vector $z$, the reduced susceptance Laplacian of $\hat{\mathcal{G}}$ can be expressed as

Here, each $a_{ij}\in\{0,\pm 1\}^{1\times N}$ is the row vector of the reduced branch-bus incidence matrix corresponding to line $(i,j)\in\hat{\mathcal{E}}$ [10]. Given the explicit form of $\tilde{L}(z)$, it is not hard to see that the objective in (7) is a monotone function that is minimized when all lines in $\hat{\mathcal{E}}$ are selected. Utilizing (8), problem (7) can be equivalently written as an MINLP [10]

Constraint (9b) enforces $X=\tilde{L}^{-1}(z)$, and thus, matrix $\tilde{L}(z)$ to be full-rank at optimality. Problem (9) is non-convex due to the binary nature of $z$ and the bilinear constraints in (9b). To handle the latter, we adopt McCormick linearization [14], which is briefly reviewed next; see also [10] for details.

Constraint (9b) involves bilinear terms $z_{\ell}X_{ij}$ for all $\ell\in\hat{\mathcal{E}}$ and $i,j\in\mathcal{V}_{r}$. For every term, introduce an auxiliary variable

and let the entries $X_{ij}$ lie within bounds $[\underline{X}_{ij},\overline{X}_{ij}]$. Since $z_{\ell}$ is binary, the following inequalities hold true

One can replace the bilinear terms in (9b) by $y_{\ell ij}$’s, drop equation (10), and enforce (11) as additional constraints for all $\ell\in\hat{\mathcal{E}}$ and $i,j\in\mathcal{V}_{r}$ to get an MILP reformulation of (9). It is not hard to show that this reformulation is exact, i.e., $y_{\ell ij}=z_{\ell}X_{ij}$ because $z$ is binary [14].

Through the aforementioned process, problem (9) is reformulated to an MILP over variables $\{X_{ij}\}$, $\{z_{\ell}\}$, and $\{y_{\ell ij}\}$, and can thus be handled by modern MILP solvers such as Gurobi or CPLEX. Nonetheless, MILPs with McCormick linearization can be forbiddingly complex to solve if the bounds $(\underline{X}_{ij},\overline{X}_{ij})$ on each $X_{ij}$ are arbitrarily wide.

This section develops graph theoretic arguments to provide easy-to-compute non-trivial bounds on the entries of $X$ and derive valid cuts for two topology design tasks: i) adding lines to an existing connected network; and ii) designing a new network afresh.

For a given finite set $\mathcal{S}$, a function $f:2^{\mathcal{S}}\rightarrow\mathbb{R}$ is said to be supermodular if for all subsets $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{S}$ and all $\mathcal{C}\in\mathcal{S}\backslash\mathcal{B}$, it holds that $f(\mathcal{A}\cup\mathcal{C})-f(\mathcal{A})\leq f(\mathcal{B}\cup\mathcal{C})-f(\mathcal{B})$ [20]. In other words, the returns due to selection of $\mathcal{C}$ are non-diminishing where adding elements to the larger set $\mathcal{B}$ gives larger gains. Minimizing a supermodular decreasing function is NP-hard. However, it is known that a greedy scheme that iteratively minimizes the objective $f(A\cup\{s\})$ for $s\in\mathcal{S}\backslash\mathcal{A}$ is at least $1-1/e\simeq 63\%$ close to the optimal cost [20].

In general, the edge addition problem is not supermodular and hence strong theoretical guarantees are not permissible. The following theorem lays a restrictive condition under which the set function $f(\hat{\mathcal{E}})=\text{Tr}(\tilde{W}\tilde{L}^{-1})$ is a supermodular decreasing function of edge addition.

To the best of our knowledge, this is the first time explicit conditions for supermodularity have been provided for the topology design problem in power grids. Since these conditions are restrictive and do not hold in general, a greedy scheme is unlikely to perform well.

All tests were carried out on a 2.3 GHz Intel Dual-Core i5 laptop with 16GB RAM. The MILP formulations were solved in Julia/JuMP v0.6 [21] using Gurobi v8.1.1. We first tested the performance of the Tightened MILP formulation for augmenting existing networks. For this design task, the IEEE $39$-bus system was used as the pre-existing connected network [22]. Set $\hat{\mathcal{E}}$ consisted of edges in this base network and an additional $22$ randomly placed lines. From these lines, we solved (9) for $K=\{5,6,7,8\}$. To satisfy Assumption 1, we assumed $M_{i}=10^{-4}$ on all nodes that did not host generators, and $D_{i}=c=0.025$ for all $i\in\mathcal{V}$. Table I compares the computation time of the MILP formulation in [10] with the Tightened MILP formulation discussed in this paper. On average, the run-time required to find the optimal solution to the Tightened MILP decreased by $61\%$.

We next considered the radial topology design problem with $\hat{\mathcal{E}}$ composed of all edges in the IEEE $39$-bus network. To generate the eigenvector-based cuts we utilized a threshold value of $\gamma=-0.95$ and set the sparsity parameter $k$ to one. Compared to the MILP formulation presented in [10] that required $8,034$ sec. to reach the optimal solution, the Tightened MILP formulation with (18d)-(18e) was solved in $7,361$ sec. Augmenting the previous formulation with the eigenvector-based cuts reduced the run-time further to $4620$ sec. The overall reduction in computation time was $40\%$ - a non-trivial improvement in comparison to state-of-the-art methods. Similarly, we considered the optimal design of a meshed network with $39$ edges, where $\hat{\mathcal{E}}$ was composed of all edges in the $39$-bus network. The optimal solution in this case was found in $9$ hours. Longer run-times for the latter task can be attributed to looser bounds on the entries of $X$.

We have presented tight mathematical formulations and numerous valid cutting planes that can substantially speed up the computational time required to find an optimal topology design (radial or meshed) for enhanced power system stability.