Designing Power Grid Topologies for Minimizing Network Disturbances: An Exact MILP Formulation

TL;DR

This paper presents an exact MILP formulation for designing power grid topologies that minimize network disturbances, improving robustness and stability through optimized line placement and topology configuration.

Contribution

It introduces a novel MILP-based approach for power grid topology design that effectively minimizes disturbances and enhances stability, addressing scalability issues of previous methods.

Findings

Optimal topologies reduce network disturbances in IEEE test cases

The MILP approach outperforms traditional design methods in stability metrics

Graphical properties improve computational efficiency of the optimization

Abstract

The dynamic response of power grids to small transient events or persistent stochastic disturbances influences their stable operation. This paper studies the effect of topology on the linear time-invariant dynamics of power networks. For a variety of stability metrics, a unified framework based on the $H_2$-norm of the system is presented. The proposed framework assesses the robustness of power grids to small disturbances and is used to study the optimal placement of new lines on existing networks as well as the design of radial (tree) and meshed (loopy) topologies for new networks. Although the design task can be posed as a mixed-integer semidefinite program (MI-SDP), its performance does not scale well with network size. Using McCormick relaxation, the topology design problem can be reformulated as a mixed-integer linear program (MILP). To improve the computation time, graphical properties are exploited to provide tighter bounds on the continuous optimization variables. Numerical tests on the IEEE 39-bus feeder demonstrate the efficacy of the optimal topology in minimizing disturbances.

Figures

Tables

Table 1. TABLE I: Cost of designed grids using the IEEE 39 39 39 -bus network.

# of lines	Optimal Topology	Suboptimal Topology
$2$	$0.570$	$0.585$
$3$	$0.564$	$0.582$
$4$	$0.557$	$0.576$
$5$	$0.552$	$0.570$

Equations

$$L_{ij}:=\{\begin{array}{ll}-b_{ij}& ,\text{if (i,j)∈E}\\ {\sum }_{(i,j)\in \mathcal{E}}{b}_{ij}& ,\text{if j=i}\\ 0& ,\text{otherwise}.\end{array}$$

$$M_i{\dot{\omega }}_i+D_i{\omega }_i=P_i^m-P_i^e+u_i$$

$$P_i^e\simeq \sum _{(i,j)\in \mathcal{E}}{b}_{ij}({\theta }_i-{\theta }_j).$$

$$\left[\begin{array}{c}\dot{\theta }\\ \dot{\omega }\end{array}\right]=\underset{A:=}{\underbrace{\left[\begin{array}{cc}0& I_N\\ -M^{-1}L& -M^{-1}D\end{array}\right]}}\left[\begin{array}{c}\theta \\ \omega \end{array}\right]+\underset{B:=}{\underbrace{\left[\begin{array}{c}0\\ M^{-1}\end{array}\right]}}u$$

$$f_c(t):=\sum _{i=1}^N{\left({\theta }_i(t)-\frac{1}{N}\sum _{j=1}^N{\theta }_j(t)\right)}^2.$$

$$f(t):=\sum _{\forall i\ne j}{w}_{ij}({\theta }_i(t)-{\theta }_j(t){)}^2+\sum _{i\in \mathcal{V}}{s}_i{\omega }_i^2(t)$$

$$y(t):=\underset{C:=}{\underbrace{\left[\begin{array}{cc}{W}^{1/2}& 0\\ 0& S^{1/2}\end{array}\right]}}\left[\begin{array}{c}\theta (t)\\ \omega (t)\end{array}\right].$$

$$\Vert H{\Vert }_{{\mathcal{H}}_2}^2:=\underset{t\to \infty }{\lim }\mathbb{E}\left[\Vert y(t){\Vert }_2^2\right].$$

$$\Vert H{\Vert }_{{\mathcal{H}}_2}^2:=\sum _{i=1}^N{\int }_0^{\infty }\Vert y_i(t){\Vert }_2^{2}dt$$

$$\Vert H{\Vert }_{{\mathcal{H}}_2}^2=\mathrm{Tr}(B^{\top }QB)$$

$$A^{\top }Q+QA=-C^{\top }C.$$

$$Q=\left[\begin{array}{cc}{Q}_1& Q_0\\ Q_0^{\top }& Q_2\end{array}\right].$$

$$\arg \underset{\mathcal{E}\in \hat{\mathcal{E}}}{\min }$$

$$s.to$$

$$Q\text{ satisfies}\text{eqref}observability$$

$$\mathcal{E}\text{ is connected}.$$

$$\mathrm{Tr}(B^{\top }QB)=\frac{\mathrm{Tr}(W{L}^{+})+\mathrm{Tr}(S{M}^{-1})}{2d}.$$

$$\arg \underset{\mathcal{E}\in \hat{\mathcal{E}}}{\min }$$

$$s.to$$

$$\mathrm{rank}(L)=N$$

$$\mathrm{Tr}(W{L}^{+})=\mathrm{Tr}(\tilde{W}{\tilde{L}}^{-1}).$$

$$W=R\tilde{W}{R}^{\top }\text{and}L=R\tilde{L}{R}^{\top }$$

$$R^{\top }R=I_N+1_N{1}_N^{\top }=J^{-1}.$$

$$L^{+}=RJ{\tilde{L}}^{-1}J{R}^{\top }.$$

$$\mathrm{Tr}(W{L}^{+})=\mathrm{Tr}\left(R\tilde{W}{R}^{\top }RJ{\tilde{L}}^{-1}J{R}^{\top }\right)=\mathrm{Tr}(\tilde{W}{\tilde{L}}^{-1})$$

$$\mathcal{Z}:=\left\{z:z^{\top }{1}_{|\hat{\mathcal{E}}|}\le K,z\in \{0,1{\}}^{|\hat{\mathcal{E}}|}\right\}.$$

$$\tilde{L}(z)=\sum _{(i,j)\in \hat{\mathcal{E}}}{z}_{ij}{b}_{ij}{a}_{ij}{a}_{ij}^{\top }$$

$$[a_{ij}{]}_n:=\{\begin{array}{ll}+1& ,\text{ if }n=i\\ -1& ,\text{ if }n=j\\ 0& ,\text{ otherwise.}\end{array}$$

$$\arg \underset{X,z\in \mathcal{Z}}{\min }$$

$$s.to$$

$$(X^{\ast },z^{\ast })\in \arg \underset{X,z\in \mathcal{Z}}{\min }$$

$$s.to$$

$$y_{mij}=z_m{X}_{ij}.$$

$$z_m(X_{ij}-{\underline{X}}_{ij})\ge 0$$

Full text

I Introduction

The electric power system is continually changing. It is expected that the grid of the future will have higher variability due to renewables, changing load patterns, and distributed energy sources [1]. This paradigm shift will pose an enormous challenge for design and stable operation of power networks. The inherent uncertainty associated with renewable energy sources and active loads is likely to produce more frequent and higher amplitude disturbances [1]. In addition, owing to the lower aggregate inertia of systems with high penetration of renewables, the capability of power networks to handle such disturbances may be significantly reduced [2].

Thus, improving the dynamic performance of the power grid is of importance and has received greater attention from academia and industry. Efforts in this direction include development of payment structures and novel markets, as well as analysis of techniques to incentivize load-side participation [3], [4], [5], [6]. The benefits of load-side controllers has motivated a series of recent works to understand how different system parameters and controller designs impact the transient response of the network [7], [8]. Recent work in [9] has also explored the optimal placement of virtual inertia to improve stability.

Compared to the previously mentioned system parameters, the effect of network topology on transient stability is less well understood. Without detailed simulations, it is usually hard to infer how a change in network topology influences the overall grid behavior and performance. Recent work in [8] shows that the impact of network topology on the power system can be quantified through the network Laplacian matrix eigenvalues. In addition, grid robustness against low frequency disturbances is mostly determined by network connectivity [8], further motivating this study. Past studies in the power and control systems communities have also looked at designing network topologies for specific goals using system theoretic tools. Such goals include reduction of transient line losses [1], improvement in feedback control [10], [11], coherence based network design [12] and augmentation [13]. Semidefinite programming (SDP) based tools have also been utilized to design and augment network topologies for dynamic control [14], [12], [15].

While the primary focus here is topology design, we recognize that there is line of related work dealing with learning network topologies and line parameters [16], [17], [18]. Schemes that rely on passive data have been used in [18] and [16] for learning radial topologies. Different from them, the work in [17] actively probes the grid to recover radial topologies and verify line statuses.

In this paper we are interested in studying the effect of topology on the power grid dynamics. For a variety of objective functions, such as line loss reduction, fast damping of oscillations, and network coherence, our previous work [19] presented a unified framework to study topology design based on the $\mathcal{H}_{2}$-norm. In [19], the focus was on topology reconfiguration rather then topology design. Further, the work in [19] developed suboptimal algorithms, albeit with guarantees on optimality gap, to tackle the combinatorial design problems involved. Here, we present reformulations of the topology design task that allow us to solve the problem to optimality.

Our contributions are as follows. First, we provide a comprehensive modeling and analysis framework for the topology design problem to optimize a $\mathcal{H}_{2}$ norm based performance metric subject to budget constraints in Section II. Second, in Section III we show that although the topology design task is inherently non-convex, it is possible to exactly reformulate the problem in tractable form using McCormick relaxation (or linearization). This can then be used with off-the-shelf solvers to determine the optimal solution. Further, we show that exploiting graph-theoretic properties to tighten bounds on the continuous optimization variables yields significant improvements in computation time. Section V discusses numerical tests based on the IEEE $39-$bus test case followed by conclusions and future directions in Section VI.

Notation: Column vectors (matrices) are denoted by lower- (upper-) case letters and sets by calligraphic symbols, unless noted otherwise. The cardinality of set $\mathcal{X}$ is denoted by $|\mathcal{X}|$. Given a real-valued sequence $\{x_{1},x_{2},\ldots,x_{N}\}$, $x\in\mathbb{R}^{N}$ is the vector obtained by stacking the scalars $x_{i}$ and $\operatorname{dg}(\{x_{i}\})$ is the corresponding diagonal matrix. The operator $(\cdot)^{\top}$ stands for transposition. The $N$-dimensional all ones vector is denoted by $1_{N}$; $I_{N}$ is the $N\times N$ identity matrix; and $e_{i}$ is the canonical vector with a $1$ at the $i$-th entry and zero everywhere else. The notation $\dot{\theta_{i}}$ denotes its time derivative $\frac{\delta\theta_{i}}{\delta t}$.

II Grid Modeling and Network Coherence

This section presents a linearized power system model, defines network coherence metrics, and reviews their connection to studying the stability of a linear dynamical system.

II-A Dynamic Power System Model

A power network having $N+1$ buses can be modeled as a connected graph ${\cal G}=({\cal V},\cal{E})$, whose nodes ${\cal V}:=\{0,1,\ldots,N\}$ correspond to buses, and edges ${\mathcal{E}}\subseteq\cal V\times\cal V$ to undirected lines. Bus $i=0$ is selected as the reference; all other buses constitute set $\mathcal{V}_{r}:=\mathcal{V}\backslash\{0\}$. Let $b_{ij}>0$ be the susceptance of line $(i,j)\in\mathcal{E}$ connecting nodes $i$ and $j\in\cal V$. Then, the susceptance Laplacian matrix $L\in\mathbb{R}^{(N+1)\times(N+1)}$ associated with the grid graph $\cal G$ can be defined as

[TABLE]

Each node $i\in\mathcal{V}$ is associated with a phase angle $\theta_{i}$, frequency $\omega_{i}=\dot{\theta_{i}}$, inertia constant $M_{i}$, and damping coefficient $D_{i}$; see [20]. Since bus $i$ may host an ensemble of devices such as synchronous machines, renewable or energy storage sources, frequency-dependent or actively controlled frequency-responsive loads, the parameters $M_{i}$ and $D_{i}$ characterize their aggregate behavior [9].

Focusing on small-signal stability, the quantities $(\theta_{i},\omega_{i})$ will henceforth refer to the deviations of nodal voltage angles and frequencies from their nominal values. The grid dynamics at bus $i$ can then be described by the linearized swing equation [20]

[TABLE]

where $P^{m}_{i}$ denotes the mechanical power input; $P_{i}^{e}$ is the electric power flowing from bus $i$ to the grid; and $u_{i}$ is the power consumed at bus $i$. Again, the aforementioned quantities refer to the deviations from their scheduled values. Under the linearized DC model, the power flowing from bus $i$ to the grid can be approximated as [20]

[TABLE]

Combining (1) and (2), the state-space representation of the power grid is

[TABLE]

where ${M}:=\operatorname{dg}(\{M_{i}\})$ and ${D}:=\operatorname{dg}(\{D_{i}\})$ are diagonal matrices containing the inertia and damping coefficients; the states ${\omega}\in\mathbb{R}^{N+1}$ and ${\theta}\in\mathbb{R}^{N+1}$ are accordingly the stacked vectors of nodal frequencies and angles; and ${u}\in\mathbb{R}^{N+1}$ is the vector of local power disturbances. The subsequent analysis relies on the ensuing assumption.

Assumption 1**.**

The inertia coefficients are strictly positive and damping coefficients are identical for all buses, that is $M_{i}>0$ and $D_{i}=d$ for all $i\in\cal V$.

The non-zero inertia assumption is not necessary, but simplifies our presentation. If $M_{i}=0$ for a bus $i\in\cal V$, then a Kron-reduced network containing only nodes with inertia can be obtained. The second assumption of constant damping has been previously used in [1], [10]. When it is not satisfied, the stability metric defined in the next section does not enjoy a closed-form expression, and only bounds can be derived. In our future work, we plan to extend our approach to the case with variable $D$.

II-B Generalized Network Coherence Metrics

Given the state-space model in (3), our goal is to design network topologies or augment existing ones to minimize the voltage angle deviations caused by load disturbances. These angle deviations are formally captured by the metric of network coherence [21]. The latter is defined as $\lim_{t\rightarrow\infty}\mathbb{E}[f_{c}(t)]$ where $f_{c}(t)$ is the steady-state deviation of the angles from their grid-average

[TABLE]

In other words, network coherence quantifies how tightly bus angles drift together. Larger variances in angle deviations reflect a more disordered network [1].

Instead of network coherence, the system operator may be interested in minimizing the expected steady-state value of a generalized function combining both voltage angle and frequency excursions [19], [9]

[TABLE]

where $w_{ij}$ and $s_{i}$ are given non-negative scalars. The weights $w_{ij}$ induce a connected weighted graph $\mathcal{G}_{w}$ that is not necessarily identical to $\cal G$; see Fig. 1. Let $W$ be the Laplacian matrix of graph $\mathcal{G}_{w}$ and $S:=\operatorname{dg}(\{s_{i}\})$. Then, it is not hard to verify that $f(t)=\|y(t)\|_{2}^{2}$ where

[TABLE]

Being a Laplacian matrix, matrix $W$ is positive semidefinite, and so its matrix square root $W^{1/2}$ is well-defined. The matrix $S^{1/2}$ is well-defined too. The importance of the generalized coherence metric of (5) is that for different choices of $(W,S)$, one can represent different grid performance metrics and study them in a unified manner [19]: For instance, if the system operator is only interested in minimizing frequency excursions, one can set $W=0$ and $S=I_{N+1}$. Similarly, if the goal is to reduce transient line losses, one can choose $W=L$ and $S=0$. Lastly, the network coherence metric of (4) corresponds to the case of $W=I_{N+1}-\frac{1}{N+1}1_{N+1}1_{N+1}^{\top}$ and $S=0$. Figure 1 shows the physical power network and its associated coherence graph. Note that the choice of $(W,S)$ for network coherence penalizes global deviations, whereas that for line loss reduction penalizes local deviations.

II-C Relation to Stability Analysis

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 (3) and (6). This system will be compactly denoted by $H:=(A,B,C)$. Leveraging this link, the generalized network coherence is amenable to a closed-form expression [19].

The ${\cal H}_{2}$ norm is widely used as a stability performance metric and has several interpretations [22]: For unit-variance stochastic white noise $u(t)$, the ${\cal H}_{2}$-norm of an LTI system equals the steady-state output variance [22, Ch. 5], [1], [9]

[TABLE]

For unit-impulse disturbances $u_{i}(t)=e_{i}\delta(t)$ for $i=1,\ldots,N$, the ${\cal H}_{2}$-norm can be equivalently written as

[TABLE]

where $y_{i}(t)$ is the system output corresponding to disturbance vector $u_{i}(t)$. Instead of evaluating the expectation in (7) or the time integral in (8), the $\mathcal{H}_{2}$-norm for system $H$ can be expressed as [22, Ch. 5]

[TABLE]

where $Q:=\int_{0}^{\infty}e^{A^{\top}t}C^{\top}Ce^{At}dt$ is the observability Gramian matrix of the LTI system $H$. In fact, the matrix $Q$ can be computed as the solution to the linear Lyapunov equation [22, Ch. 5]

[TABLE]

Matrix $Q$ is known to be symmetric positive semidefinite, so it can be partitioned as

[TABLE]

Based on the reformulation in (9), we next design topologies attaining minimum generalized network coherence.

III Optimal Grid Topology Design

Among other criteria, the topology of a grid can be designed to minimize the generalized network coherence metric of (5). Given a graph $\hat{\mathcal{G}}=(\mathcal{V},\hat{\mathcal{E}})$, where $\hat{\mathcal{E}}$ is the set of candidate lines weighted by their susceptances, the goal is to find the subset ${\mathcal{E}}\subseteq\hat{\mathcal{E}}$ of cardinality $|\mathcal{E}|\leq K$ with $K\geq N$ attaining the minimum generalized network coherence. Given the equivalences of the previous section, this task can be posed as the optimization problem

[TABLE]

The constraint in (12b) reflects the budget on the number of edges. For $K=N$, the problem in (12) yields a tree topology, which is important for typically radially operated distribution grids. Interestingly, leveraging the problem structure under Assumption 1, it can be shown that the objective of (12) simplifies as [9], [19], [1]

[TABLE]

If Assumption 1 is not met, i.e., the damping coefficients are not identical $(D\neq dI_{N})$, then it may not be possible to find a closed-form expression for the objective in (12). If the damping coefficients are known to lie within $[d_{\min},d_{\max}]$, then the objective of (12) can be bounded as $\frac{\operatorname{Tr}(WL^{+})+\operatorname{Tr}(SM^{-1})}{2d_{\max}}\leq\operatorname{Tr}(B^{\top}QB)\leq\frac{\operatorname{Tr}(WL^{+})+\operatorname{Tr}(SM^{-1})}{2d_{\min}}$; see [1], [9]. Therefore, as the range $[d_{\min},d_{\max}]$ becomes narrower, minimizing the numerator of (13) approaches the optimal solution to (12). From Assumption 1, we consider $d_{\min}=d_{\max}$ here.

The second summand in the right-hand side of (13) is independent of the grid topology, and can thus be ignored. Problem (12) can then be reformulated as

[TABLE]

where the rank constraint ensures that the graph induced by $\mathcal{E}$ is connected. Problem (14) could be tackled with brute-force algorithms over all the possible topologies of budget $K$ and below; though that would incur exponential complexity.

The objective in (14) can be written in terms of the inverse of the reduced Laplacian matrix of $\mathcal{G}$ as explained next. The claim has appeared in [1, Lemma 3.2], albeit for the restricted case where $W$ and $L$ have the same structure.

Lemma 1**.**

If $\tilde{W}$ and $\tilde{L}$ are the $N\times N$ matrices obtained after removing the first row and column from $W$ and $L$, respectively, then

[TABLE]

Proof.

The Laplacian matrices can be described in terms of their reduced counterparts as

[TABLE]

where $R:=[1_{N}~{}-I_{N}]^{\top}$. This is because a Laplacian matrix satisfies $L1_{N+1}=0_{N+1}$. If we define matrix $J:=I_{N}-\frac{1}{N+1}1_{N}1_{N}^{\top}$, then it is not hard to see that

[TABLE]

Then the pseudoinverse of $L$ can be expressed as

[TABLE]

The latter can be shown by simply verifying that $LL^{+}L=L$ and $L^{+}LL^{+}=L^{+}$. From (15) and (17), we get that

[TABLE]

where the last equality follows from (16) and the cyclic property of the trace. ∎

To express the optimization in (14) over $\hat{\mathcal{E}}$ in a more convenient form, let us introduce a binary variable $z_{m}$ for every line $m\in\hat{\mathcal{E}}$. This variable is $z_{m}=1$ if line $m$ is selected, i.e., $m\in\mathcal{E}$; and $z_{m}=0$, otherwise. If we stack variables $\{z_{m}\}_{m\in\hat{\mathcal{E}}}$ in vector $z$, then $z$ has to lie in the set

[TABLE]

Based on the line selection vector $z$, the reduced susceptance Laplacian of $\mathcal{G}$ can be expressed as

[TABLE]

where each vector $a_{ij}$ corresponds to line $(i,j)\in\hat{\mathcal{E}}$, and its $n$-th entry is defined as

[TABLE]

Given Lemma 1 and (19), the optimization in (14) can be equivalently written as the mixed-integer semidefinite program (MI-SDP)

[TABLE]

To see this, the constraint in (20b) is equivalent to $X\succeq 0$ and $X\succeq\tilde{L}^{-1}(z)$; see [23, Sec. A.5.5]. Since $\tilde{W}\succ 0$, the optimal $X$ can be shown to be $X=\tilde{L}^{-1}(z)$. In fact, constraint (20b) waives the possibility of the optimal $\tilde{L}$ being singular, and thus, ensures the connectedness of the graph.

To overcome the computational complexity of the MI-SDP in (20), one may relax the binary variables to box constraints as $z\in[0,1]^{|\hat{\mathcal{E}}|}$ to get an ordinary SDP, which can be handled by off-the-shelf solvers for moderately-sized networks. Being a relaxation, the SDP solution provides a lower bound on the cost of (20). If the obtained solution of the SDP turns out to be binary, then this $z$ minimizes the MI-SDP in (20) as well. Otherwise, (randomized) rounding heuristics can be adopted to acquire a feasible $z$.

Aiming at an exact solver, we will next show how the MI-SDP of (20) can be equivalently formulated as an MILP. To this end, we first rewrite (20) as

[TABLE]

Note that for constraint (21b) to hold, $\tilde{L}(z^{*})$ must be non-singular and $X^{*}=\tilde{L}^{-1}(z^{*})$. Although its cost is linear, problem (21) is non-convex due to the bilinear constraints in (21b) and because vector $z$ is binary. To handle the former, we adopt the McCormick relaxation technique [24], which is briefly reviewed next.

Constraint (21b) involves the bilinear terms $z_{m}X_{ij}$ for $m\in\hat{\mathcal{E}}$ and $i,j\in\mathcal{V}_{r}$. For each such term, introduce a new variable

[TABLE]

Suppose that the entries $X^{*}_{ij}$ are known to lie within the range $[\underline{X}_{ij},\overline{X}_{ij}]$. Since $z_{m}\in[0,1]$, the ensuing inequalities hold trivially [24]

[TABLE]

Substituting $z_{m}X_{ij}$ by $y_{mij}$ in (23) provides

[TABLE]

One can replace the bilinear terms in (21b) by $y_{mij}$’s and enforce (22) and (24) as additional constraints for all $m\in\hat{\mathcal{E}}$ and $i,j\in\mathcal{V}_{r}$. In that case, the constraints (24) are apparently redundant. However, one may simplify the problem by dropping (22) to get an MILP reformulation of (21). Interestingly, this reformulation is exact due to the binary nature of $z$. To see this, if $z_{m}^{*}=1$ for some $m$, then (24b) and (24d) imply $y_{mij}^{*}=X^{*}_{ij}$ for all $i,j\in\mathcal{V}_{r}$. Otherwise, if $z_{m}^{*}=0$, then (24a) and (24c) imply $y_{mij}^{*}=0$ for all $i,j\in\mathcal{V}_{r}$.

Through the aforementioned process, problem (21) has been converted to an MILP over the variables $\{X_{ij}\}$, $\{z_{m}\}$, and $\{y_{mij}\}$. MILPs can be handled by state-of-the-art solvers such as Gurobi [25], and are known to scale better than MI-SDPs. Recall that the MILP reformulation of (21) requires the bounds $(\underline{X}_{ij},\overline{X}_{ij})$ on each $(i,j)$-th entry of $X^{*}$. Lacking prior information on $X^{*}$, one could select an arbitrarily wide range $[\underline{X}_{ij},\overline{X}_{ij}]$. However, this could slow down the MILP solver significantly. In the other extreme, if $X^{*}$ is known, that is $\underline{X}_{ij}=\overline{X}_{ij}$ for all $i,j\in\mathcal{V}_{r}$, then the binary vector $z$ can be recovered by simply solving the linear equations in (21b). To improve the run time of the involved MILP, we next derive tighter, non-trivial bounds on the entries of $X_{ij}^{*}$.

IV Bounding the Continuous Variables

Depending on the problem structure, different bounds can be derived on $X^{*}_{ij}$s. This section considers two classes of topology design tasks. In the first task, some lines are already energized and the operator would like to augment a connected network by additional lines to further improve its stability. In the second task, a network topology is designed from scratch with the additional requirement of a radial grid.

IV-A Augmenting Existing Networks

Given the graph $\hat{\mathcal{G}}=(\mathcal{V},\hat{\mathcal{E}})$, this problem setup considers a pre-existing connected network described by $\mathcal{G}_{e}=(\mathcal{V},\mathcal{E}_{e})$, and the goal is to energize additional lines from $\hat{\mathcal{E}}\setminus\mathcal{E}_{e}$ to improve its stability. In essence, this corresponds to the problem in (21) with the entries of $z$ corresponding to the lines in $\mathcal{E}_{e}$ being set to one. Then, based on (19), the reduced Laplacian matrix of the existing network is obviously

[TABLE]

Under this setup, the entries of $X^{*}$ minimizing (21) under the additional constraints $z_{\ell}=1$ for all $\ell\in\mathcal{E}_{e}$, can be bounded as follows.

Lemma 2**.**

The entries of $X^{*}$ are bounded by

[TABLE]

Proof.

Observe that $\tilde{L}_{e}$ can be written as $\tilde{L}(z)$ by fixing the entries of $z$ corresponding to the lines in $\mathcal{E}_{e}$ to $1$. Since the same entries will remain $1$ in $z^{*}$, it readily follows that $\tilde{L}(z^{*})\succeq\tilde{L}_{e}\succ 0$, where $\tilde{L}_{e}$ is non-singular since the existing network is already connected. Then, it follows that $X^{*}\preceq\tilde{L}_{e}^{-1}$ and $v^{\top}(X^{*}-\tilde{L}_{e}^{-1})v<{0}$ for any $v\neq{0}$. Setting $v=e_{i}$, the diagonal entries of $X^{*}$ can be bounded as $X^{*}_{ii}\leq[\tilde{L}_{e}^{-1}]_{ii}$ for all $i\in\mathcal{V}_{r}$. Since $X^{*}\succ 0$ by constraint (21b), we have that $({e}_{i}-{e}_{j})^{\top}X^{*}({e}_{i}-{e}_{j})>0$ for all $i,j\in\mathcal{V}_{r}$. The latter provides $X^{*}_{ij}<\frac{1}{2}(X^{*}_{ii}+X^{*}_{jj})$. Upper bounding the diagonal entries with the bounds obtained earlier proves the upper bound in (25). The lower bound in (25) can be trivially obtained since $L(z^{*})$ is an M-matrix, and so its inverse has positive entries. ∎

The reduced Laplacian matrix $\tilde{L}_{e}$ of the existing network $\mathcal{G}_{e}$ is invertible as long as $\mathcal{G}_{e}$ is connected. If that is not the case, one could obtain bounds on $X^{*}_{ij}$’s by imposing a radial structure on the sought topology as discussed next.

IV-B Radial Topology Design

The setup considered here designs a network afresh, but under the requirement that it is radial. The analysis simplifies under the following mild assumption.

Assumption 2**.**

There exists a node in $\mathcal{V}$ that is incident to exactly one edge.

To derive bounds on the entries of $X^{*}$ minimizing (21) for the special case of $K=N$ (radial network), let us first construct the graph $\mathcal{G}_{r}=(\mathcal{V},\hat{\mathcal{E}}_{r})$, where $\hat{\mathcal{E}}_{r}$ consists of the edges in $\hat{\mathcal{E}}$, but with inverse weights $x_{ij}:=b_{ij}^{-1}$ for all $(i,j)\in\hat{\mathcal{E}}$. Based on $\mathcal{G}_{r}$, we define some additional properties that will be useful later. If one of the nodes in $\mathcal{V}$ satisfying Assumption 2 is selected as the reference node, then the weight (inverse susceptance) of its incident line is denoted by $x_{0}$. Let us also define the minimum of the inverse line susceptances

[TABLE]

Before solving (21), we find the maximum spanning tree of graph $\mathcal{G}_{r}$, and denote the sum of its edge weights by $f$. The maximum spanning tree can be found efficiently by finding the minimum spanning tree on $\mathcal{G}_{r}$ upon negating its edge weights [26]. Lastly, for each node $i\in\mathcal{V}_{r}$, we find its shortest path to the reference node [math] in $\mathcal{G}_{r}$. The sum of edge weights along this shortest path will be denoted by $h_{i}$.

Lemma 3**.**

Under Assumption 2, the entries of $X^{*}$, minimizing (21) for $K=N$, are bounded as

[TABLE]

Proof.

The bounds rely on a fundamental property of the inverse Laplacian matrix of a radial network: If $\tilde{L}$ is the reduced susceptance Laplacian of a radial network and $X:=\tilde{L}^{-1}$, then the entry $X_{ij}$ equals the sum of the inverse susceptances that are common to the paths of nodes $i$ and $j$ to the reference bus; see [27, Lemma 1]. Under Assumption 2, the common path of any pair of nodes $(i,j)\in\mathcal{V}_{r}$ must include at least the line incident to the reference bus, and hence, $X^{*}_{ij}\geq x_{0}$. By the definition of shortest path, the entry $X^{*}_{ii}$ is lower bounded by $h_{i}$ for all $i\in\mathcal{V}_{r}$, thus providing the lower bound in (26a).

Regarding the upper bound in (26a), recall that the $(i,i)$-th entry of $X^{*}$ equals the sum of weights on the path from $i$ to the reference node [math]. The sum of weights on the longest such path is still upper bounded by the sum weight $f$ of all edges in the maximum spanning tree. Because the entry $X^{*}_{ij}$ for $i\neq j$ must have at least one less edge than the longest path, the upper bound in (26b) holds as well. ∎

Note that the upper bounds can be tight in the setting where the maximum spanning and optimum trees are the same line graph.

IV-C Model Simplification and Bound Tightening

Exploiting the structure of $\hat{\mathcal{G}}=(\mathcal{V},\hat{\mathcal{E}})$ can provide additional information on the bounds of $X^{*}$ matrix entries to accelerate solving (21) for the special case of $K=N$ (radial grid). For example, if $\hat{\mathcal{G}}$ gets disconnected upon removing edge $\ell\in\hat{\mathcal{E}}$, this edge $\ell$ belongs to the sought tree topology and $z_{\ell}^{*}=1$ before solving (21). To identify such edges, we resort to the notion of a graph cutset $\mathcal{C}\subset\hat{\mathcal{E}}$, defined as a subset of edges that once removed, splits the graph $\hat{\mathcal{G}}$ into two or more connected components. The edges in this cutset will be termed as critical edges. Now, we present a simple algorithm to enumerate all the single critical edges $(|\mathcal{C}|=1)$ by initializing the graph $\hat{\mathcal{G}}$’s weights to unit values.

1. Step 1:

Solve the min-cut problem on $\hat{\mathcal{G}}$ with unit edge weights by using the standard max-flow min-cut algorithm  [28]. 2. Step 2:

Increase the weight for the edge labelled as critical to $1+\epsilon$ for some $\epsilon>0$. 3. Step 3:

If $|\mathcal{C}|>1$, quit; else, go to Step 1.

The second step ensures that every time we are identifying a new critical edge. Upon completing this process, the entries of $z^{*}$ corresponding to the critical edges can be safely set to $1$. This process not only reduces the binary search for $z^{*}$ in (21), but it further tightens the lower bounds on certain $X^{*}_{ij}$’s as discussed in Lemma 4; see also Fig. 2.

Lemma 4**.**

Suppose a critical edge $\ell=(i,j)\in\hat{\mathcal{E}}$ partitions the nodes of $\hat{\mathcal{G}}$ into two disjoint connected components $\mathcal{V}_{\ell}$ and its complement $\bar{\mathcal{V}}_{\ell}$. If $\mathcal{V}_{\ell}$ contains node $i$ as well as the reference bus, then (26b) can be tightened as

[TABLE]

Proof.

The edge $(i,j)\in\hat{\mathcal{E}}$ is the only edge that connects the nodes in $\bar{\mathcal{V}}_{\ell}$ to the rest of the network. Hence, the common path to the reference bus of any two nodes in $\bar{\mathcal{V}}_{\ell}$ must include the path of node $j$ to the reference. It follows that the shortest path weight $h_{j}$ is a valid lower bound for all nodes in $\bar{\mathcal{V}}_{\ell}$. ∎

Identifying cutsets of cardinality larger than $1$ offers additional information to tighten the bounds of entries of the $X$ matrix. If graph $\hat{\mathcal{G}}$ gets disconnected upon removing lines $\ell_{1},\ell_{2}\in\hat{\mathcal{E}}$, then at least one of these lines should be active. This logical conclusion translates to the constraint $z_{\ell_{1}}+z_{\ell_{2}}\geq 1$, which can be augmented to (21) to tighten the McCormick reformulation and possibly accelerate the MILP solver. Cutsets of larger cardinality, say $|\mathcal{C}|=k,k>1$, can be identified by iterating Steps 1 through 3 of the algorithm described earlier. In this case, we assign the weights of $k+\epsilon$ on the critical edges, where $\epsilon>0$.

V Numerical Tests

All tests were run on a 2.7 GHz Intel Core i5 laptop with 8GB RAM. The MILP formulations were solved using Gurobi v8.0.1 optimizer, written in Julia/JuMP  [25, 29].

The performance of the MILP in (21) was tested for augmenting an existing network as well as for designing a radial one afresh. For the augmentation setup, the IEEE $39$-bus system benchmark was used as the pre-existing connected network [30]. The set $\hat{\mathcal{E}}$ was selected by $10$ randomly picked additional lines. From these lines, we solved the restricted version of (21) for $K=\{2,3,4,5\}$. To satisfy Assumption 1, we assumed $M_{i}=10^{-4}$ on all buses that did not host generators, and $D_{i}=d=0.025$ for all $i\in\cal V$. To grade an arbitrarily constructed network, we compared its squared ${\cal H}_{2}$ norm to that of the optimal network of the same edge cardinality. The results are summarized in Table I. Additional lines are useful in minimizing disturbances, and so the budget constraint in (21) was always met with equality. For all cases, the augmentation design problem was solved in less than 5 seconds.

We next considered the radial topology design problem with $\hat{\mathcal{E}}$ composed of all edges in the IEEE $39$-bus network. The optimal cost obtained after solving the MILP in this case was $1.669$, and the time required to find the optimal tree was close to $2$ hours. Considering that the problem needs to be solved once off-line, this running time may not be of concern. Figure 4 shows the optimal radial topology that was identified.

Instead of using the bounds of Lemma 3 and the bound tightening procedure of Section IV-C, we attempted to solve the same MILP with the relatively looser bounds of $0\leq X^{*}_{ij}\leq 10$ for all $(i,j)\in\cal V$. In this case, the solver reached the optimality gap of 60% after running for $3$ hours. Clearly, having tighter bounds improves the computation time.

For a 1 per unit impulse input at bus $39$, Figure 3 compares the frequency response of the optimal tree (left panel) with an arbitrarily selected suboptimal tree (right panel). Notice that not only is the amplitude of oscillations lower in the left panel, but the generators also seem to drift together. These observations indicate that the optimal tree selected is much more effective at minimizing the effect of small disturbances.

VI Conclusions and Future Work

We have presented a general framework to study the effect of network topology on power grid dynamics. Using a system-theoretic approach, we have shown that the original topology design problem (MI-SDP) can be reformulated in tractable form as an MILP. To improve the computation time, graph-theoretic properties have been exploited to simplify the model and provide tighter bounds on the continuous optimization variables involved. Numerical tests on the IEEE 39-bus benchmark suggest that meshed networks exhibit improved coherence behavior. Current research efforts are focused on tackling the topology design task with non-uniform damping, exploring conditions for submodularity, and considering the design problem from an $\mathcal{H}_{\infty}$-norm perspective.