Transient Stability of Droop-Controlled Inverter Networks with Operating Constraints

TL;DR

This paper analyzes the transient stability of droop-controlled inverter networks with operational constraints, providing criteria to ensure frequency synchronization and robustness against disturbances, using Lyapunov functions and cycle flow information.

Contribution

It introduces novel stability criteria incorporating cycle flows to reduce conservativeness and quantify robustness in inverter networks under constraints.

Findings

Less-conservative stability conditions achieved.

Cycle flow information improves stability analysis.

Quantified robustness to parameter disturbances.

Abstract

Due to the rise of distributed energy resources, the control of networks of grid-forming inverters is now a pressing issue for power system operation. Droop control is a popular control strategy in the literature for frequency control of these inverters. In this paper, we analyze transient stability in droop-controlled inverter networks that are subject to multiple operating constraints. Using a physically-meaningful Lyapunov-like function, we provide two sets of criteria (one mathematical and one computational) to certify that a post-fault trajectory achieves frequency synchronization while respecting operating constraints. We show how to obtain less-conservative transient stability conditions by incorporating information from loop flows, i.e., net flows of active power around cycles in the network. Finally, we use these conditions to quantify the scale of parameter disturbances to which the network is robust. We illustrate our results with numerical case studies of the IEEE 24-bus system.

Equations

$$d_{\mathrm{cc}}(\alpha ,\beta )=\{\begin{array}{ll}|\alpha -\beta |,& \text{c.c. arc from α to β shorter than π}\\ -|\alpha -\beta |,& \text{otherwise}\end{array}$$

$$V_{<}^{-1}(c)=\{x{\in }^n\mid V(x)<c\}.$$

$$B_{i,e}=\{\begin{array}{ll}+1,& s(e)=i\\ -1,& t(e)=i\\ 0,& \text{otherwise}\end{array}$$

$$\Delta (\gamma )=\left\{\theta \in {\mathbb{T}}^n:|B^{\mathsf{T}}\theta |<\gamma \right\}$$

$${\dot{\theta }}_i={\omega }^{\ast }-\frac{p_{e,i}-p_i^{\ast }}{d_i},\forall i\in {\mathcal{V}}_I$$

$$p_{e,i}=p_i^{\ast }-d_i({\dot{\theta }}_i-{\omega }^{\ast }),\forall i\in {\mathcal{V}}_L$$

$$p_{ij}^{\mathrm{line}}={\tilde{a}}_{ij}+a_{ij}\sin ({\theta }_i-{\theta }_j-{\varphi }_{ij})$$

$$p_{e,i}=\sum _{j\in \mathcal{N}(i)}{\tilde{a}}_{ij}+a_{ij}\sin ({\theta }_i-{\theta }_j-{\varphi }_{ij}),\forall i\in \mathcal{V}$$

$$D\dot{\theta }=f(\theta )$$

$$f_i(\theta )=p_i-\sum _{j\in \mathcal{N}(i)}{a}_{ij}\sin ({\theta }_i-{\theta }_j-{\varphi }_{ij}),\forall i\in \mathcal{V}$$

$$\left|{\int }_0^{\infty }{p}_e(t)-p^{\ast }dt\right|\le \bar{s}$$

$$v(\theta )=\dot{\theta }-{\omega }^{\ast }\text{mathbbold}{1}_n=D^{-1}f(\theta )-{\omega }^{\ast }\text{mathbbold}{1}_n$$

$$V(\theta )=||v(\theta )|{|}_{\infty }$$

$${\gamma }_e^{\ast }=\frac{\pi }{2}-\max \left\{|{\varphi }_{ij}|,|{\varphi }_{ji}|\right\},\forall e=\{i,j\}\in \mathcal{E}$$

$$\underset{i\in I_{\max }}{\max }\{-\mathrm{sgn}(x_i)(Lx{)}_i\}\le 0$$

$$-\mathrm{sgn}(x_i)(Lx{)}_i$$

$$\le \sum _{j\ne i}{L}_{ij}(|x_i|-|x_j|)$$

$$J_{ij}(\theta )=\frac{\partial f_i(\theta )}{\partial {\theta }_j}=\{\begin{array}{ll}-\cos ({\theta }_i-{\theta }_j-{\varphi }_{ij}),& \{i,j\}\in \mathcal{E}\\ 0,& \text{else}\end{array}$$

$$V(t)=\underset{i\in I_{\max }}{\max }\left\{\mathrm{sgn}(v_i(t))v_i(t)\right\}$$

$${\text{D}}^{+}V(t)$$

$$=\underset{i\in I_{\max }}{\max }\left\{\mathrm{sgn}(v_i(t)){\left(D^{-1}J(\theta )v(t)\right)}_i\right\}$$

$$\begin{array}[]{rl}\text{minimize}:&V(\theta)\\ \text{variables}:&\theta\in\mathbb{T}^{n}\\ \text{subject to}:&\theta\in\partial S\\ &D^{-1}f(\theta)~{}\text{is pointed outward from}~{}S\end{array}$$

$${\delta }_0\le \frac{1}{2}{\bar{r}}_i{\left(\sum _{j\in \mathcal{N}(i)}{a}_{ij}\right)}^{-1}$$

$${\delta }_0\le \frac{{\lambda }_2(L)\cos ({\gamma }_{\max }){\bar{s}}_i}{d_i/d_{\min }}{\left(1+\frac{1}{2}\log \left(\frac{d_{\mathrm{sum}}}{d_{\min }}\right)\right)}^{-1}$$

$$|{\dot{p}}_{e,i}|=\left|\sum _{j\in \mathcal{N}(i)}{a}_{ij}\cos ({\theta }_i-{\theta }_j-{\varphi }_{ij})({\dot{\theta }}_i-{\dot{\theta }}_j)\right|\le 2{\delta }_0\sum _{j\in \mathcal{N}(i)}{a}_{ij}$$

$$\frac{d}{dt}v(\theta {)}^{\mathsf{T}}Dv(\theta )=2v(\theta {)}^{\mathsf{T}}J(\theta )\dot{\theta }=2v(\theta {)}^{\mathsf{T}}J(\theta )v(\theta )$$

$$\frac{d}{dt}v(\theta {)}^{\mathsf{T}}Dv(\theta )$$

$$\le -\cos ({\gamma }_{\max }){\lambda }_2(L)d_{\min }\left(v(\theta {)}^{\mathsf{T}}Dv(\theta )\right)$$

$$\min$$

$$f\in {\mathbb{R}}^n,y\in {\mathbb{R}}^m,{\eta }^{+}\in {\mathbb{R}}^m,{\eta }^{-}\in {\mathbb{R}}^m,$$

$$z^{+}\in \{0,1{\}}^m,z^{-}\in \{0,1{\}}^m$$

$$f=p-B^{+}{A}^{+}{\eta }^{+}-B^{-}{A}^{-}{\eta }^{-}$$

$${\eta }_e^{+}=\sin (y_e-{\varphi }_{s(e),t(e)}),\forall e\in \mathcal{E}$$

$${\eta }_e^{-}=-\sin (y_e+{\varphi }_{t(e),s(e)}),\forall e\in \mathcal{E}$$

Full text

I Introduction

Transient stability is a power systems problem of both practical importance and theoretical interest. The goal of transient stability analysis is to determine whether or not the system will return to a stable, frequency-synchronized operating point after a large disturbance. Transients are difficult to analyze: the governing differential equations are nonlinear, and linearization techniques are not useful for large-scale disturbances. Therefore, system operators typically rely on numerical simulation [13, Chapter 9.3] to study system behavior. Simulation is an effective tool for analyzing individual disturbance scenarios, but it has limitations. Simulating a comprehensive set of disturbances is computationally expensive, and it does not establish rigorous guarantees.

Direct methods of transient stability analysis address these limitations by establishing theoretical guarantees on transient behavior. Direct methods are not a substitute for simulation in real-world power system operation, since they rely on low-order, theoretically-tractable models. Instead, they provide significant theoretical insight into these simplified systems. Classical works on using Lyapunov-like methods to study transient stability include [2, 24, 8]. More recently, [25, 26] used set-theoretic control techniques to establish regions of attraction for the coupled swing equations. Direct methods are highly model-specific and provide conservative guarantees, so this topic is still the subject of active research.

Historically, the literature on direct methods has focused on networks of high-inertia synchronous generators. But the rise of distributed energy resources has sparked a growing interest in the stability of low-inertia inverter networks, particularly microgrids. Inertia is both a blessing and a curse from a control perspective—the same inertia that makes the system robust to disturbances also makes the system respond sluggishly to control inputs. A suitable fast-acting controller can make a low-inertia inverter network highly robust. Two broad classes of inverter controllers have emerged to exploit this low inertia: grid-following controllers, in which the inverter acts as a current source to track the local voltage signal; and grid-forming controllers, in which the inverter acts as a voltage source to stabilize voltage frequencies throughout the network. Both of these frameworks involve new models and require fresh approaches to direct transient stability analysis.

One of the most popular approaches to grid-forming control is proportional droop control, in which local voltage frequencies are modulated in proportion to the power drawn from neighboring buses. Recent work [23, 1, 29] has studied the dynamics of droop-controlled microgrids (DCMGs) via the inhomogeneous Kuramoto model. Under certain assumptions, equilibrium points of the Kuramoto model correspond to frequency-synchronized operating points of the DCMG, and regions of attraction around these equilibria provide a rigorous way to assess how robust DCMG operating points are to disturbances. Some progress has been made on estimating these regions of attraction [11], but these closed-form estimates tend to be very conservative and require stringent regularity assumptions on the topology or system parameters.

Another limitation of the literature is that few bounds on the transients are available. To a system operator, a guarantee of frequency synchronization alone is not satisfying, if the resulting transient will violate operating constraints (like constraints on line flows and nodal power injections). Recent work has begun to address transient stability in conjunction with other engineering constraints [20, 19]. If direct methods of transient stability are to provide more insight into the operation of DCMGs, then less-conservative regions of attraction, as well as bounds on quantities of engineering significance, are needed. This paper addresses these two needs.

Contributions

Our first contribution is to extend the transient stability problem. In addition to the classical notion of transient stability (asymptotic frequency synchronization), we impose five “desired properties” of transients, so as to enforce operating constraints on nodal frequencies, power flows across transmission lines, nodal power injections, nodal ramping rates, and reserves of stored energy.

Our second contribution is to provide two sufficient conditions for when a trajectory of a DCMG will exhibit transient stability and the five desired properties. Both certificates require only two pieces of information from from the initial condition (nodal frequencies and line angle differences) instead of the full (and harder to measure) vector of voltage angles. The first certificate can be viewed as a DCMG-specific form of Nagumo’s theorem, and it is intended as a theoretical basis for transient-stability-certifying algorithms. The second certificate, which consists of a tractable mixed-integer linear program (MILP), is built on top of the first certificate. These theoretical results use a physically-meaningful Lyapunov function called the “maximum frequency deviation,” which (to our knowledge) has not been used before to study power systems.

Our third contribution is to improve these two certificates using the winding partition of the $n$-torus. We introduced the winding partition in [17] to localize the multiple equilibrium points of network flows on the $n$-torus. This paper provides the first application of the winding partition to analyzing system dynamics (in contrast to its previous applications to statics problems). We show how to incorporate the “winding vector” of the initial condition (a quantity closely related to flows of active power around cycles in the network) into the two certificates, resulting in less-conservative conditions for transient stability and the other desired properties.

As a fourth contribution, we use our transient stability conditions to quantify the size of parameter disturbances with respect to which the DCMG is robust. We define a single number that quantifies the “size” of an arbitrary change in model parameters, and we compute a critical threshold such that post-fault transient stability is guaranteed in any disturbance “smaller” than the threshold. We examine particular disturbance modes, including changes in nominal power injections, voltage magnitudes, and branch admittance magnitudes. We illustrate all of these results numerically, using the IEEE 24-bus system as a case study.

Organization

The next four sections are organized around our four main contributions. After introducing our notation, model, and problem statement, Section II states the extended transient stability problem in Definition 1, formally introducing transient stability and the five desired properties of the transient. Section III presents our two certificates in Theorem 4 and Theorem 6, respectively. We review the winding partition in Section IV and improve both certificates by incorporating the winding vector in Theorem 9, and we show that these improved certificates are less conservative (Theorem 10). Finally, Section V uses the stability certificates to study how robust a DCMG is to changes in parameters, and it presents our numerical case studies.

II Preliminaries and Problem Statement

II-A Notation

The Circle and $n$-Torus

Let $\mathbb{S}$ be the circle, i.e., the set of phases or angles. For every pair of angles $\alpha,\beta\in\mathbb{S}$, we use $|\alpha-\beta|$ to denote the geodesic distance between them. The counterclockwise difference between two angles is the map $d_{\rm cc}:\mathbb{S}\times\mathbb{S}\to[-\pi,\pi)$, where

[TABLE]

In other words, we consider the clockwise and counterclockwise arcs from $\alpha$ to $\beta$. If the counterclockwise arc is shorter, then $d_{\rm cc}(\alpha,\beta)$ is the length of that arc. Otherwise, $d_{\rm cc}(\alpha,\beta)$ is the negated length of the clockwise arc. The $n$-torus, denoted $\mathbb{T}^{n}$, is the product of $n$ circles.

Sets

Given any set $S$ within $\mathbb{R}^{n}$ or $\mathbb{T}^{n}$, we refer to the interior of the set by $\mathrm{int}(S)$, the closure by $\mathrm{cl}(S)$, and the boundary by $\partial S=\mathrm{cl}(S)\setminus\mathrm{int}(S)$, with respect to the standard topologies on $\mathbb{R}^{n}$ and $\mathbb{T}^{n}$. Given a function $V:^{n}\to$ and a scalar $c\in$, we define a sublevel set

[TABLE]

Linear Algebra

The vector $\mathbbold{1}_{n}$ (resp. $\mathbbold{0}_{n}$) is a vector in n with all the entries equal to one (resp. zero). For every $v\in^{n}$, $\operatorname{diag}(v)\in^{n\times n}$ is a diagonal matrix with entries $\operatorname{diag}(v)_{ii}=v_{i}$ for every $i\in\{1,\ldots,n\}$. The $\infty$-norm of $v$ is $\|v\|_{\infty}=\max_{i}|v_{i}|$, and the $1$-norm of $v$ is $\|v\|_{1}=\sum_{i=1}^{n}|v_{i}|$. We define $v_{\rm sum}=\sum_{i=1}^{n}v_{i}$ and $v_{\rm min}=\min_{i}\{v_{i}\}$. For every $v,w\in^{n}$, we write $v\leq w$ (resp. $v<w$) if $v_{i}\leq w_{i}$ (resp. $v_{i}<w_{i}$), for every $i\in\{1,\ldots,n\}$. For a matrix $X\in^{n\times n}$, the Moore–Penrose pseudoinverse is denoted by $X^{\dagger}$.

Graph Theory

An undirected graph is a pair $G=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is a set of $n$ nodes, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of $m$ edges. The neighborhood of any node $i\in\mathcal{V}$ is denoted by $\mathcal{N}(i)$. While $G$ is undirected, we may enumerate and assign an arbitrary orientation to each edge $e\in\mathcal{E}$ by labeling one incident node as the “source” $s(e)$ and the other as the “sink” $t(e)$. The incidence matrix of the graph [6, §9.1] is the matrix $B\in\{-1,0,1\}^{n\times m}$ with entries

[TABLE]

Graphs and the $n$-Torus

We may assign a phase-valued state to every node in $G$, so that the full state of the graph is in $\mathbb{T}^{n}$. Given a state $\theta\in\mathbb{T}^{n}$, we use the abuse of notation $B^{\mathsf{T}}\theta$ to represent the vector in $\mathbb{R}^{m}$ of counterclockwise differences across each edge, i.e., $(B^{\mathsf{T}}\theta)_{e}=d_{\rm cc}(\theta_{i},\theta_{j})$, where $i$ is the source of $e$ and $j$ is the sink. Furthermore, given any $\gamma\in(0,\pi]^{m}$, we define the phase-cohesive set as the open set

[TABLE]

In contrast to the literature, where $\Delta(\gamma)$ takes a scalar-valued $\gamma$, the $\gamma$ we refer to in this paper is always a vector, allowing for inhomogeneous phase cohesion.

II-B Model

We consider a DCMG on an undirected topology $G=(\mathcal{V},\mathcal{E})$, with $n$ buses $\mathcal{V}=\{1,\dots,n\}$ and $m$ branches (or lines) $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. We assume that $G$ is connected, but otherwise we make no assumptions about its structure; both trees and cyclic graphs are acceptable.

Bus Model

Each bus has a complex voltage $E_{i}e^{j\theta_{i}}$, where $E_{i}>0$ is the voltage magnitude and $\theta_{i}\in\mathbb{S}$ is the phase. We assume that voltage controllers are operating at a much faster time scale than frequency controllers, so that $E_{i}$ is constant but $\theta_{i}$ is dynamic. We consider two types of buses: droop-controlled inverters and frequency-dependent loads. Buses in $\mathcal{V}_{I}\subset\mathcal{V}$ represent droop-controlled inverters, which produce a controllable voltage signal with constant magnitude $E_{i}$ and time-varying frequency $\dot{\theta}_{i}$. These inverters operate according to the frequency droop control law [7, 16]:

[TABLE]

Here $\dot{\theta}_{i}$ is the instantaneous AC frequency, $\omega^{*}$ is the nominal frequency (for example, 60 Hz), $p_{i}^{*}\geq 0$ is the nominal active power injection, and $d_{i}^{-1}>0$ is the droop coefficient. Buses in $\mathcal{V}_{L}=\mathcal{V}\setminus\mathcal{V}_{I}$ represent frequency-dependent loads [21, §9.1], where the instantaneous active power injection $p_{e,i}$ is

[TABLE]

Here $p_{i}^{*}\leq 0$ is the nominal active power load, and $d_{i}>0$. Note that (1) and (2) are algebraically equivalent.

Branch Model

For each branch $\{i,j\}\in\mathcal{E}$, we assume that the real power flow from node $i$ into the $\{i,j\}$ branch is

[TABLE]

where $\tilde{a}_{ij}\in$, $a_{ij}\geq 0$ and $\phi_{ij}\in(-\frac{\pi}{2},\frac{\pi}{2})$ are constants. These constants are not necessarily symmetric (i.e., $\phi_{ij}\neq\phi_{ji}$), so in general $p_{ij}\neq-p_{ji}$.

The AC steady-state active power flow across many types of branches can be written in the form (3). Transmission lines, for example, are typically represented by the nominal $\Pi$ model, which consists of a series admittance $Y_{ij}e^{j\varphi_{ij}}$ that is flanked by two shunt admittances $Y_{ii}e^{j\varphi_{ii}}$ and $Y_{jj}e^{j\varphi_{jj}}$ [21, §6.1]. Active power flow in the nominal $\Pi$ model is given by (3) with $\tilde{a}_{ij}=E_{i}^{2}(Y_{ii}\cos(\varphi_{ii})+Y_{ij}\cos(\varphi_{ij}))$, $a_{ij}=E_{i}E_{j}Y_{ij}$, and $\phi_{ij}=\varphi_{ij}+\frac{\pi}{2}$. In medium-length and short-length lines, the shunt admittance is typically purely capacitive or altogether negligible, either of which leads to the simplification $\tilde{a}_{ij}=E_{i}^{2}Y_{ij}\cos(\varphi_{ij})$. It is also typical that the series admittance is primarily inductive, so $\phi_{ij}\approx 0$. In the extreme case of lossless lines (with no shunt admittance), the active power flow reduces to the antisymmetric form $p_{ij}^{\rm line}=a_{ij}\sin(\theta_{i}-\theta_{j})$. For transformers, active power flow in the standard equivalent circuit model can also be written in the form (3). We omit the specifications of the parameters for brevity and refer the reader to [21, §6.2].

Dynamics

Due to conservation of energy, active power injections at each bus must balance against power outflows:

[TABLE]

Substituting this expression for $p_{e,i}$ into (1) and (2) leads to a differential equation in $\theta$ that captures the angle dynamics of the grid. We can write these dynamics compactly by defining a constant vector $p\in\mathbb{R}^{n}$ with entries $p_{i}=p_{i}^{*}+\omega^{*}d_{i}-\sum_{j\in\mathcal{N}(i)}\tilde{a}_{ij}$ for each $i\in\mathcal{V}$, as well as a matrix $D=\operatorname{diag}\{d_{i},\;i\in\mathcal{V}\}$. Then the system can be written as

[TABLE]

where $f:\mathbb{T}^{n}\to\mathbb{R}^{n}$ is a vector with entries

[TABLE]

Equation (4) is the model that we study in this paper. If the sine coefficients are homogeneous and the underlying graph is complete, this model is familiar in the physics community as the Kuramoto-Sakaguchi model [22], which has been used to study synchronization phenomena in coupled oscillator networks [28, 10, 5].

Limitations of the Model

Our model is based on several commonly-used simplifying assumptions that should be examined explicitly. Perhaps the most important simplification is that we neglect voltage dynamics and reactive power. This is particularly common in the controls community, and it is often justified by assuming fast-acting voltage controllers [26, 27]. If voltage control fails, possibly due to insufficient reactive power, then unmodeled dynamics of the $a_{ij}$ and $\tilde{a}_{ij}$ parameters may destabilize the system.

Another simplification is our use of steady-state AC models for branches in (3), which is very common in analysis of conventional power grids. These models assume sinusoidal nodal voltages at a constant frequency, an assumption that is technically contradicted by the dynamic frequencies in (4). But the purpose of this paper is to find sufficient conditions for “safe” transient stability, and a key aspect of safe power grid operation is a tight tolerance around the nominal frequency, typically under 1%. In other words, the trajectories that we are interested in certifying have only a small variance in frequency. Nonetheless, the effects of transmission line dynamics on inverter-based grids is a subject of recent interest, and we refer the reader to [14] for a rigorous study of this topic.

II-C Problem Statement

Under normal operation, nodal frequencies are synchronized at the nominal frequency $\omega^{*}$, and power injections $p_{e}$ are equal to the nominal power injections $p^{*}$. But contingencies, like failing transmission lines or a sudden change in power supply or demand, disrupt this equilibrium behavior. Droop control will stabilize the post-fault system about a new equilibrium, provided that this new equilibrium is sufficiently close to the pre-fault state. This local stability property is well-known and easily verified by eigenvalue analysis of (4).

Unfortunately, the dynamics of droop control after larger-scale disturbances are not as well understood, and local stability alone does not inspire confidence in a power system controller. Furthermore, the controller should ensure that the system’s critical engineering constraints are satisfied during the transient. In this paper, in addition to non-local transient stability, we consider five engineering constraints that are important in the context of inverter networks:

Definition 1** (Desired Properties).**

We define the following six properties that are desirable in a trajectory $\theta(t)$ of (4):

1. (P1)

Transient stability.* Nodal frequencies asymptotically synchronize, i.e., $\lim_{t\to\infty}\dot{\theta}(t)=\omega_{\textup{syn}}\mathbbold{1}_{n}$ for some synchronous frequency $\omega_{\textup{syn}}\in\mathbb{R}$.* 2. (P2)

Frequency constraint.* Nodal frequencies are bounded by $|\dot{\theta}(t)-\omega^{*}\mathbbold{1}_{n}|\leq\bar{\delta}$ for all $t\geq 0$, where $\bar{\delta}\geq\mathbbold{0}_{n}$ is a vector of frequency tolerances.* 3. (P3)

Angle difference constraint.* Voltage angle differences are bounded by $|B^{\mathsf{T}}\theta(t)|\leq\bar{\gamma}$ for all $t\geq 0$, where $\bar{\gamma}\in(0,\frac{\pi}{2}]^{m}$ is a vector of angle difference tolerances.* 4. (P4)

Power constraint.* Power injections are sufficiently close to the nominal injection, i.e., $|p_{e}(t)-p^{*}|\leq\bar{p}_{e}$ for all $t\geq 0$, where $\bar{p}_{e}\in\mathbb{R}_{\geq 0}^{n}$ is a vector of power tolerances.* 5. (P5)

Ramping constraint.* The rate of change in power injections is sufficiently small: $|\dot{p}_{e}(t)|\leq\bar{r}_{e}$ for all $t\geq 0$, where $\bar{r}_{e}\in\mathbb{R}_{\geq 0}^{n}$ is a vector of ramping tolerances.* 6. (P6)

Energy constraint.* The difference from nominal energy injection is bounded by*

[TABLE]

where $\bar{s}\in\mathbb{R}_{\geq 0}^{n}$ is a vector of nodal capacities to store or dump energy.

Each of these desired properties are necessary for safe operation of the power system. (P1) and (P2) are the standard objectives of primary and secondary frequency control, which keep nodal frequencies close to the rated frequency of grid components. (P3) protects transmission lines from overheating, since larger angle differences lead to larger current flow, and thus, more thermal dissipation. (P4) and (P6) ensure that the power and energy drawn from inverters are within a reasonable range. For example, an inverter powered by solar panels on a sunny afternoon is more flexible in its active power injection (via curtailment) than the same inverter on a cloudy morning. Finally, (P5) ensures that the rate at which power injections fluctuate is within the tolerance of the inverter.

The objective of this paper is to find computationally-tractable sufficient conditions on $\theta(0)$ for each of these six properties.

III Main Theoretical Results

We now proceed with our main results: two sets of sufficient conditions to certify that a trajectory satisfies transient stability and the five desired properties in Definition 1.

III-A Lyapunov Function

Our analysis is based on the frequency deviation vector, which measures the difference between instantaneous and nominal frequencies at each bus:

[TABLE]

From $v(\theta)$, we define our Lyapunov candidate function, the maximum frequency deviation

[TABLE]

If a trajectory $\theta(t)$ is clear from context, we will abuse notation and write $V(t)$ instead of $V(\theta(t))$. We will show that the min-max frequency deviation is non-increasing when voltage angle differences are sufficiently small; exactly how small depends on $\phi_{ij}$. For each branch, we define a critical arc length

[TABLE]

Due to the assumption that $\phi_{ij}\in(-\frac{\pi}{2},\frac{\pi}{2})$, the critical arc lengths satisfy the bound $\gamma_{e}^{*}\in(0,\frac{\pi}{2}]$, and the maximum value of $\frac{\pi}{2}$ is achieved by lossless transmission lines (for which $\phi_{ij}=0$). We collect the critical arc lengths into a vector $\gamma^{*}\in^{m}$.

As long as the angle difference across each branch is less than the critical arc length, the maximum frequency deviation is non-increasing:

Lemma 2** (Max Frequency Deviation is Non-Increasing).**

Let $\theta(t)$ be a trajectory of (4) such that $\theta(t)\in\Delta(\gamma^{*})$ on some interval $t\in[t_{0},t_{1}]$. Then $V(t_{1})\leq V(t_{0})$.

The proof of this property is based on the following lemma, which is (to our knowledge) novel:

Lemma 3** (Sign-Definiteness of Laplacian Matrices).**

Let $L\in\mathbb{R}^{n\times n}$ be a Laplacian matrix. For any $x\in\mathbb{R}^{n}$, let $I_{\rm max}=\{i:|x_{i}|=||x||_{\infty}\}$ be the set of nodes with maximal absolute value. Then

[TABLE]

Furthermore, if the digraph corresponding to $L$ is strongly connected, then equality holds if and only if $x\in\operatorname{span}(\mathbbold{1}_{n})$.

Proof of Lemma 3.

For every $i\in I_{\rm max}$, we compute

[TABLE]

The first line follows because $L$ has zero row sums and the second line because off-diagonal entries of $L$ are non-positive. But $i\in I_{\rm max}$ implies that $|x_{i}|-|x_{j}|\geq 0$, so we conclude that $-\operatorname{sgn}(x_{i})(Lx)_{i}\leq 0$. Equality clearly holds in the case where $x\in\operatorname{span}(\mathbbold{1}_{n})$. Now suppose that $\max_{i\in I_{\rm max}}\{-\operatorname{sgn}(x_{i})(Lx)_{i}\}=0$, which implies that $\sum_{j\neq i}L_{ij}(|x_{i}|-\operatorname{sgn}(x_{i})x_{j})=0$ for some particular $i\in I_{\rm max}$. But each summand is non-positive, so $|x_{i}|=\operatorname{sgn}(x_{i})x_{j}$ for each $j$ for which $L_{ij}\neq 0$; consequently, $x_{i}=x_{j}$ for each out-neighbor $j$ of $i$. It follows that $j\in I_{\rm max}$. Extending the same argument to $j$ and all of its neighbors, we see that if a directed path exists from $i$ to any node $k$, then $x_{k}=x_{i}$. But the graph is strongly connected, so we conclude that $x\in\operatorname{span}(\mathbbold{1}_{n})$. ∎

Proof of Lemma 2.

We first observe that $\ddot{\theta}=D^{-1}J(\theta)\dot{\theta}$, where $J(\theta)$ is the Jacobian matrix of $f(\theta)$. For $i\neq j$,

[TABLE]

If $\theta\in\Delta(\gamma^{*})$, then $\cos(\theta_{i}-\theta_{j}-\phi_{ij})>0$. Furthermore, evaluating the diagonal entries of $J(\theta)$ reveals that the matrix has zero row sums, so $-J(\theta)$ is the Laplacian matrix of a weighted, directed graph whose topology is identical to $G$ (treating each undirected edge in $G$ is a pair of directed edges). Note that this graph is strongly connected.

Let $I_{\rm max}=\{i:|v_{i}(t)|=||v(t)||_{\infty}\}$ be the set of buses with maximal frequency deviation, so we can write

[TABLE]

Using [6, Lemma 15.16(iii)] to compute the upper right Dini derivative of a pointwise-maximum function, we obtain

[TABLE]

where the last step follows because $\omega^{*}\mathbbold{1}_{n}\in\ker(J)$. But $D^{-1}J(\theta)$ is a Laplacian matrix corresponding to a strongly connected graph, so by Lemma 3, $D^{+}V(t)\leq 0$. Lemma 2 follows from this bound [6, Lemma 15.16(ii)]. ∎

In summary, the maximum frequency deviation is positive definite about the subspace of frequency-synchronized states, and it is non-increasing inside of $\Delta(\gamma^{*})$.

III-B Set-Theoretic Certificate

We now use the maximum frequency deviation to establish a set-theoretic transient stability and operating constraint certification. Our approach is to construct forward-invariant sets using $V$, based on the following optimization problem:

Problem 1** (Min-Max Frequency Deviation).**

Let $S\subseteq\Delta(\gamma^{*})$. We define $V^{*}(\partial S)$ to be the minimum value of the following:

[TABLE]

If the problem is infeasible, we define $V^{*}(\partial S)=+\infty$.

The min-max frequency deviation is the minimum value of $V(\theta)$ along the “outward boundary” of $S$, i.e., the portion of $\partial S$ where $\dot{\theta}$ is pointed away from the set.

Minimizing a Lyapunov function around a set boundary is a well-established technique for constructing forward-invariant sets—see, for example, Nagumo’s 1942 theorem [3, Theorem 4.7]. More recently, [26] applied this technique to a quadratic Lyapunov function for the coupled swing equations. In our case, the min-max frequency deviation is defined so that sets of the form $S\cap V^{-1}_{<}(V^{*}(\partial S))$ are forward invariant. This observation, together with the monotonicity of $V$, leads to the central theorem of the paper.

Theorem 4** (Set-Theoretic Certificate).**

Let $\theta(t)$ be a trajectory of (4). Let $\gamma_{0}=|B^{\mathsf{T}}\theta(0)|$ and $\delta_{0}=V(\theta(0))$ denote the initial angle differences and max frequency deviation. If there exist a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ and a set $\Delta(\gamma_{0})\subseteq S\subseteq\Delta(\gamma)$ such that $\delta_{0}<V^{*}(\partial S)$, then

1. (i)

$\theta(t)\in S\cap V_{<}^{-1}(\delta_{0})$* for all $t\geq 0$.* 2. (ii)

The transient stability property (P1) is satisfied.

Further conditions on $\gamma$ and $\delta_{0}$ lead to various desirable properties from Definition 1:

1. (iii)

The frequency constraint (P2) is satisfied for each bus $i$ if $\delta_{0}\leq\bar{\delta}_{i}$. 2. (iv)

The angle difference constraint (P3) is satisfied for each line $\{i,j\}$ if $\gamma_{ij}\leq\bar{\gamma}_{ij}$. 3. (v)

The power constraint (P4) is satisfied for each bus $i$ if $\delta_{0}\leq\bar{p}_{e,i}d_{i}^{-1}$. 4. (vi)

The ramping constraint (P5) is satisfied for each bus $i$ if

[TABLE]

Additionally, in the special case of lossless networks (where $a_{ij}=a_{ji}$ and $\phi_{ij}=\phi_{ji}=0$), the following is true:

1. (vii)

The energy constraint (P6) is satisfied for each bus $i$ if

[TABLE]

where $\lambda_{2}(L)$ is the smallest non-zero eigenvalue of the Laplacian matrix $L=B\left(\mathcal{\operatorname{diag}}\{a_{ij}\}_{\{i,j\}\in\mathcal{E}}\right)B^{\mathsf{T}}$.

Proof.

To prove statement (i), observe that any trajectory which escapes $S$ must cross through some point on $\partial S$ where $\dot{\theta}$ is pointed outward from $S$. By definition, $V^{*}(\partial S)\leq V(\theta)$ at such a point $\theta$. But Lemma 2 implies that $V(\theta)\leq V(0)$, which further implies that $\theta(0)\notin V_{<}^{-1}(\delta_{0})$ if the trajectory reaches this point. Forward invariance of $S\cap V_{<}^{-1}(\delta_{0})$ follows by contrapositive. Regarding statement (ii), recall from the proof of Lemma 2 that the frequency dynamics can be written $\ddot{\theta}=D^{-1}J(\theta)\dot{\theta}$, where $D^{-1}J(\theta)$ is the negated Laplacian matrix of a strongly connected digraph when $\theta\in\Delta(\gamma^{*})$. It follows from [6, Theorem 12.10] that $\dot{\theta}(t)$ converges to a consensus state.

Statements (iii) and (iv) follow trivially from statement (i). Statement (v) follows because droop control relates power injections to frequencies by $p_{e,i}=p_{i}^{*}-d_{i}(\dot{\theta}_{i}-\omega^{*})$ for each $i\in\mathcal{V}$. Therefore $|p_{e,i}(t)-p_{i}^{*}|\leq d_{i}V(t)\leq d_{i}\delta_{0}$ for all $t\geq 0$. To prove statement (vi), we observe for each bus $i$ that

[TABLE]

since $\cos(\theta_{i}-\theta_{j}-\phi_{ij})\in(0,1)$. To prove (vii), observe that

[TABLE]

where the last step follows because $\ker(J(\theta))=\operatorname{span}\{\mathbbold{1}_{n}\}$. Under the lossless assumption, $J(\theta)$ is negated symmetric Laplacian matrix, and the edge weights in the corresponding graph are $a_{ij}\cos(\theta_{i}-\theta_{j})$, which is lower-bounded by $a_{ij}\cos(\gamma_{\rm max})$. It follows from [6, Lemma 6.9(ii)] that $\lambda(-J(\theta))\geq\cos(\gamma_{\rm max})\lambda_{2}(L)$, so

[TABLE]

Therefore $v(\theta)^{\mathsf{T}}Dv(\theta)$ has an exponential upper bound, which decays in time at the rate $\cos(\gamma_{\rm max})\lambda_{2}(L)d_{\rm min}$. The integrand in (P6) can be upper-bounded using both $\delta_{0}$ and this exponential, yielding the condition in statement (vii). ∎

Theorem 4 simplifies transient analysis in two ways. First, the conditions depend on the quantities $\gamma_{0}$ and $\delta_{0}$, rather than the full initial state $\theta(0)$. A system operator can measure $\gamma_{0}$ through line flows and $\delta_{0}$ through nodal frequencies, rather than using state estimation to obtain $\theta(0)$. Second, the theorem recasts transient analysis as the search for a set $S\subseteq\Delta(\gamma^{*})$ with a sufficiently large min-max frequency deviation. The remainder of the section examines a computationally-efficient way to search for such a set.

III-C Groundwork for the MILP Certificate

It is impractical to repeatedly evaluate $V^{*}(\partial S)$ while searching for a set that satisfies Theorem 4. Fortunately, we can efficiently compute upper bounds on $V^{*}(\partial S)$ if we restrict our search to sets of the form $S=\Delta(\gamma)$. We obtain these upper bounds through a series of relaxations to Problem 1, and then we use these bounds to establish an easily-computable transient stability certificate. This subsection lays out the first of two relaxations that we make for this certificate.

When $S=\Delta(\gamma)$ for some $\gamma\in(0,\gamma^{*}]$, Problem 1 admits the following relaxation:

Problem 2** (Min-Max Frequency Deviation, Lower Bound).**

Let $\gamma\in(0,\gamma^{*}]$. We define $\widehat{V}(\gamma)$ to be the minimum value of the following:

[TABLE]

If the problem is infeasible, we define $\widehat{V}(\gamma)=+\infty$.

Recall that $s(e)$ and $t(e)$ represent the arbitrary “source” and “target” nodes of each $e\in\mathcal{E}$. To express constraint (6b) succinctly, we decompose the incidence matrix $B$ into two matrices $B^{+},B^{-}\in\{0,1\}^{n\times m}$, where $(B^{+})_{i,e}=1$ if and only if $s(e)=i$, and $(B^{-})_{i,e}=1$ if and only if $t(e)=i$, so that $B=B^{+}-B^{-}$. We also define two diagonal matrices $A^{+}=\operatorname{diag}\{a_{s(e),t(e)},\;\forall e\in\mathcal{E}\}$ and $A^{-}=\operatorname{diag}\{a_{t(e),s(e)},\;\forall e\in\mathcal{E}\}$. Constraints (6f) and (6g) are indicator constraints: if $z_{e}^{+}=1$, then the constraints $y_{e}=\gamma_{e}$ and $d_{s(e)}^{-1}f_{s(e)}-d_{t(e)}^{-1}f_{t(e)}\geq 0$ become “active,” but these constraints do not apply if $z_{e}^{+}=0$. Indicator constraints are easily encoded in the MILP framework [4], and many MILP solvers allow indicator constraints to be supplied explicitly.111CPLEX 12.9 supports explicit indicator constraints. Similarly, the Python interface to Gurobi 9 provides the method Model.addGenConstrIndicator(). Both links accessed 8/9/2020.

Problem 2 relaxes Problem 1 by optimizing the vector of counterclockwise differences $y=B^{\mathsf{T}}\theta$ directly, instead of optimizing $\theta\in\mathbb{T}^{n}$. Given this interpretation of $y$, constraints (6b)–(6d) ensure that $f=f(\theta)$ and that the cost function (6a) is equal to $V(\theta)$. Constraint (6e) guarantees that $\theta\in\Delta(\gamma)$, and (6f)–(6h) ensure that the underlying $\theta$ is on the “outward-pointing” boundary of $S$.

The most important property of Problem 2 is that it yields a lower bound to $V^{*}(\partial\Delta(\gamma))$:

Lemma 5** (Problem 2 is a Relaxation).**

Let $\gamma\in(0,\gamma^{*}]$, and let $S=\Delta(\gamma)$. The solutions to Problems 1 and 2 are related by $\widehat{V}(\gamma)\leq V^{*}(\partial S)$, where equality holds if the underlying graph $G$ is a tree.

The proof is contained in Appendix A. Due to this bound, we can replace the $\delta_{0}<V(\Delta(\gamma))$ condition in Theorem 4 with the stricter (but computable) condition $\delta_{0}<\widehat{V}(\gamma)$.

Theorem 6** (Computational Certificate).**

Consider a trajectory $\theta(t)$ of (4) on any connected graph $G$. Let $\gamma_{0}=|B^{\mathsf{T}}\theta(0)|$ and $\delta_{0}=V(\theta(0))$ denote the initial angle differences and initial max frequency deviation. If there exists a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ such that $\delta_{0}<\widehat{V}(\gamma)$, then statements (i)–(vii) from Theorem 4 hold, with respect to the set $S=\Delta(\gamma)$.

Proof.

Let $S=\Delta(\gamma)$. By Lemma 5, $\delta_{0}<\widehat{V}(\gamma)\leq V^{*}(\partial S)$, so $\gamma$ and $S$ satisfy the hypothesis of Theorem 4. ∎

Given a particular $\gamma\in[\gamma_{0},\gamma^{*}]$, Theorem 6 provides a certificate for transient stability and the other operating constraints in Definition 1, using only two properties of the initial condition: the initial angle differences $\gamma_{0}$, and the initial maximal frequency deviation $\delta_{0}$. Furthermore, Theorem 6 replaces Problem 1 with Problem 2, which is more readily solved by numerical methods.

But two issues still remain. The first problem with Theorem 6 is that it is conservative, since Problem 2 is a lower bound on Problem 1. This bound is only tight in acyclic networks, and the gap between these two problems tends to increase with the number of edges in the graph. In other words, denser graphs lead to more conservative certificates provided by Theorem 6 (compared to Theorem 4 applied to $S=\Delta(\gamma)$). Closing this gap with additional constraints requires some deeper analysis of the $n$-torus geometry, which we postpone to Section IV.

The second issue is that Problem 2 is difficult to solve. It is possible to tackle the problem using nonlinear programming techniques, but it is faster and safer to relax the problem to obtain a lower bound, which we examine next.

III-D MILP Certificate

Problem 2 is challenging to solve numerically, since it contains the nonlinear equality constraints (6c) and (6d). Furthermore, we must be careful about using nonlinear solvers to estimate $\widehat{V}$. If the solver obtains a sub-optimal solution, then this solution may exceed $V^{*}(\partial\Delta(\gamma))$, thereby invalidating Theorem 6. But any lower bound on $\widehat{V}(\gamma)$ can be used in place of $\widehat{V}(\gamma)$ in Theorem 6. We can get a lower bound—while simultaneously making the problem much easier to solve—by further relaxing Problem 2 into a MILP.

The relaxation is conceptually simple; all we need to do is replace (6c) and (6d) with linear and integer constraints. The general idea is to find a polytope or a union of polytopes that contain the sine curve. Then these bounds can be encoded within the MILP and substituted in for (6c) and (6d). The process of constructing these polytopes is a messy exercise in elementary geometry, so we do not go into details here. Instead, we present two examples in Figure 1, both of which relax the sine constraints with four linear bounds (although tighter and more complicated bounds are clearly possible.) In fact, one can achieve arbitrary precision in the relaxed sine constraints by using piecewise-linear bounds, at the expense of additional binary variables and slower computation.

Replacing (6c) and (6d) with the polytope relaxation turns Problem 2 into a MILP. This MILP can be solved with standard software like Gurobi, CPLEX, or MATLAB. The solution is a lower bound on $\widehat{V}(\gamma)$, which can safely be used in place of $\widehat{V}(\gamma)$ in Theorem 6. We also note that this MILP is computationally tractable, since the binary variables essentially split the problem into $2m$ linear programming sub-problems, each corresponding to one of the $2m$ faces of $\Delta(\gamma)$. Informally, Problem 2 is no more complex than a collection of $2m$ linear programs.

IV Improved Guarantees for Meshed Networks

In the previous section, we found transient guarantees based on two properties of the initial condition, $\gamma_{0}$ and $\delta_{0}$. But if $G$ is a cyclic topology (i.e., the graph is not a tree), we can make use of an additional property of the initial condition: its winding vector, $u_{0}$. Winding vectors have recently gained attention in the study of power transmission networks due to their relationship with loop flows [18, 9], i.e., net flows of power around cycles in the network. In [17], we partitioned the $n$-torus into equivalence classes of winding vectors, and we showed that a wide class of network systems with phase-valued states (including the Kuramoto model) have at most one equilibrium point within each equivalence class. Since winding vectors contain enough information to uniquely characterize equilibrium points of the system, it seems they may also provide information about the transient behavior.

In this section, we briefly review concepts related to the winding partition. We then apply these concepts to the transient stability problem at hand, using knowledge of the initial winding vector $u_{0}$ to obtain less-conservative certificates from Theorem 4 and completely close the gap between Problems 1 and 2. For simplicity, we assume throughout this section that $G$ contains at least one cycle. (Otherwise $G$ is a tree, in which case there is no gap between Problems 1 and 2 to begin with.)

IV-A Winding Partition of the $n$-Torus

Preliminaries

We start with some preliminaries on algebraic graph theory and its application to graph cycles. We refer the reader to [6, §9.3] for a more detailed discussion of these concepts. A simple cycle in $G$ is a sequence of consecutive nodes, where the first and last nodes are identical, but all other nodes are distinct. Given a simple cycle $\sigma=(i_{1},i_{2},\dots,i_{n_{\sigma}},i_{1})$, the cycle vector $v_{\sigma}\in^{m}$ is defined with respect to the incidence matrix $B$ by

[TABLE]

for each $e\in\mathcal{E}$. More formally, given an adjacent pair of nodes $i_{j},i_{j+1}$ in $\sigma$, we say that $\sigma$ traverses the edge $\{i_{j},i_{j+1}\}$ positively if $B_{i_{j},e}=+1$ and $B_{i_{j+1},e}=-1$; otherwise, it traverses the edge negatively. The set of cycle vectors for all simple cycles in $G$ span a vector space, called the cycle space of $G$. A set of simple cycles $\Sigma$ is called a cycle basis if the cycle vectors corresponding to elements of $\Sigma$ are a basis for the cycle space. A cycle basis $\Sigma=\{\sigma_{1},\sigma_{2},\dots,\sigma_{|\Sigma|}\}$ can be encoded in a cycle-edge incidence matrix $C_{\Sigma}\in\{-1,0,1\}^{|\Sigma|\times m}$, where

[TABLE]

Clearly $\operatorname{\mathrm{Img}}(C_{\Sigma}^{\mathsf{T}})$ is the cycle space. Furthermore, because $G$ is a connected graph, the dimension of the cycle space is $|\Sigma|=m-n+1$.

Winding Vectors and Winding Partition

We now review basic definitions regarding winding vectors and the winding partition. The partition divides the $n$-torus into equivalence classes induced by an underlying graph $G$. These equivalence classes are defined by how many times the phase differences across basis cycles of $G$ “wind” around the unit circle:

Definition 7** (Winding Numbers, Vectors, and Cells).**

Let $\theta\in\mathbb{T}^{n}$. Given any simple cycle $\sigma$ in $G$ with $n_{\sigma}$ nodes, the winding number of $\theta$ along $\sigma$ is

[TABLE]

where the nodes in $\sigma$ are indexed $\sigma=(1,\dots,n_{\sigma},1)$ and $\theta_{n_{\sigma}+1}=\theta_{1}$. Given a cycle basis $\Sigma$ of $G$, the winding vector of $\theta$ along $\Sigma$ is the vector

[TABLE]

For every winding vector $u\in w_{\Sigma}(\mathbb{T}^{n})$, the $u$-winding cell is the equivalence class

[TABLE]

We note that winding cells were introduced in [17], but winding vectors have been used in the study of power flows since [18]. The reader has likely encountered a similar concept of “winding number” from interpreting Nyquist plots.

Winding vectors are always integer-valued, a property which is analogous to Kirchoff’s voltage law (KVL). For real-valued nodal potentials, KVL guarantees that potential differences sum to zero around any cycle. Similarly, phase differences (in the sense of counter-clockwise arc length) sum to an integer multiple of $2\pi$ around any cycle. For example, suppose that $G$ is the triangle graph, consisting of a single cycle $\sigma=(1,2,3,1)$. Let $\theta\in\mathbb{T}^{3}$. If there is an arc of length $\pi$ that contains $\theta_{1},\theta_{2},\theta_{3}$, then $w_{\sigma}(\theta)=0$. Otherwise, $w_{\sigma}(\theta)=\pm 1$. Figure 2 (top) illustrates these three possible winding numbers, based on the configuration of phases around the cycle. Meanwhile, Figure 2 (bottom) illustrates $\Omega_{u}$ for each $u=-1,0,1$.

The cycle basis of a graph is often non-unique, and each valid cycle basis $\Sigma$ leads to a different definition of the winding vector $w_{\Sigma}(\theta)$. However, it turns out that the equivalence classes of $w_{\Sigma}$, namely the winding cells, do not depend on the particular choice of cycle basis. Given two cycle bases $\Sigma$ and $\Sigma^{\prime}$, every winding cell based on $w_{\Sigma}$ is identical to a winding cell based on $w_{\Sigma^{\prime}}$.

Physical Interpretation of Winding Vectors

Winding vectors are closely connected to net flows of active power around cycles in the network. These flows, called loop flows, are of considerable interest to the power systems community because they do not deliver useful power and can jeopardize system stability. For example, flows around the Lake Erie Loop were a major factor in the 2003 Northeast Blackout [15]. Rigorous connections between loop flows and winding vectors are discussed in detail in [18, 9, 17]. We here provide some basic intuition and show how to measure the initial winding vector using line flows, instead of the full state $\theta(0)$.

Consider the active power flows across each line given by (3). We assume that these flows are measurable, so that $p^{\rm line}$ is known, even if the state $\theta$ is not. If $\theta\in\Delta(\gamma^{*})$, then (7) can be written

[TABLE]

If we ignore shunt and series losses by setting $\tilde{a}_{i,i+1}=0$ and $\phi_{i,i+1}=0$, and we expand the arcsine function about the origin, we obtain

[TABLE]

The quantity $a_{i,i+1}^{-1}p_{i,i+1}^{\rm line}$ is a normalized line flow, scaled by the capacity of the line. Thus, up to second order, the winding number is a normalized loop flow (at least in the case of short, lossless transmission lines). While somewhat informal, this analysis suggests that winding vectors are a quantized measure of these loop flows. We can also use (10) to infer the winding vector from line flow data. Therefore, like the vector of angle differences $\gamma_{0}$, we can identify the initial winding vector $u_{0}$ using measurements of line flows instead of the full state $\theta(0)$.

IV-B Improved Certificates

We have previously seen (in Theorems 4 and 6) how to certify transient stability and other desirable properties from the initial condition, using the measurable quantities $\gamma_{0}$ and $\delta_{0}$ instead of the full state $\theta(0)$. As with $\gamma_{0}$ and $\delta_{0}$, the initial winding vector $u_{0}=w_{\Sigma}(\theta(0))$ provides additional information that we can exploit to make guarantees about the transient. In the remainder of this section, we will show how to use the winding vector to obtain better certificates out of Theorems 4 and 6, replacing the antecedents of these theorems with less-conservative conditions.

The new conditions are straightforward to state and prove. We modify Theorem 4 to search over sets of the form $\Delta(\gamma_{0})\cap\Omega_{u_{0}}\subseteq S\subseteq\Delta(\gamma)\cap\Omega_{u_{0}}$ instead of $\Delta(\gamma_{0})\subseteq S\subseteq\Delta(\gamma)$. Reducing the lower bound from $\Delta(\gamma_{0})$ to $\Delta(\gamma_{0})\cap\Omega_{u_{0}}$ directly incorporates $u_{0}$ and results in a larger search space for $S$, thereby expanding the set of cases that satisfy the antecedent of Theorem 4. Shrinking the upper bound from $\Delta(\gamma)$ to $\Delta(\gamma)\cap\Omega_{u_{0}}$ is not strictly necessary, but as we will see later on, we can always find an optimal $S$ within this smaller upper bound. Similarly, we will modify Theorem 6 by adding a constraint to Problem 2 that forces the optimum to reside within $\Omega_{u_{0}}$:

Problem 3** (Min-Max Frequency Deviation, Exact).**

Let $\gamma\in(0,\gamma^{*}]$ and $u\in\operatorname{\mathrm{Img}}(w_{\Sigma})$. We define $\widehat{V}(\gamma,u)$ as the minimum value of Problem 2, under the additional constraint $C_{\Sigma}y=2\pi u$. If the problem is infeasible, we define $\widehat{V}(\gamma,u)=+\infty$.

The additional $C_{\Sigma}y=2\pi u$ constraint confines the solution to $\Omega_{u}$, completely closing the gap between Problems 1 and 2:

Lemma 8** (Relations of Minima).**

Let $\gamma\in(0,\gamma^{*}]$, let $u\in\operatorname{\mathrm{Img}}(w_{\Sigma})$, and let $S=\Delta(\gamma)\cap\Omega_{u}$. The solutions to Problems 1, 2, and 3 are related by

[TABLE]

The proof is contained in Appendix A. We can now state the new transient stability certificates that account for the initial winding vector.

Theorem 9** (Certificates with Winding Vectors).**

Consider a trajectory $\theta(t)$ of (4), and let $\Sigma$ be a cycle basis of the underlying graph. Let $\gamma_{0}=|B^{\mathsf{T}}\theta(0)|$, $\delta_{0}=V(\theta(0))$, and $u_{0}=w_{\Sigma}(\theta(0))$ denote the initial angle differences, max frequency deviation, and winding vector. Consider the following two conditions:

1. (a)

There exist a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ and a set $\Delta(\gamma_{0})\cap\Omega_{u_{0}}\subseteq S\subseteq\Delta(\gamma)\cap\Omega_{u_{0}}$ such that $\delta_{0}<V^{*}(\partial S)$. 2. (b)

There exists a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ such that $\delta_{0}<\widehat{V}(\gamma,u_{0})$.

If either (a) or (b) are true, then statements (i)–(vii) from Theorem 4 hold, with respect to $S$ from condition (a) or $S=\Delta(\gamma)\cap\Omega_{u_{0}}$ from condition (b).

Proof.

The proof that statements (i)–(vii) follow from (a) is identical to the proof of Theorem 4, since the new upper and lower bounds on $S$ do not impact the argument for statement (i), and (because $S$ is still contained within $\Delta(\gamma)$) they have no bearing on statements (ii)–(vii). We can use this result from condition (a) to prove that (i)–(vii) follow from condition (b). Let $S=\Delta(\gamma)\cap\Omega_{u_{0}}$, and observe that $\delta_{0}<\widehat{V}(\gamma,u_{0})\leq V^{*}(\partial S)$ due to Lemma 8. Then $S$ satisfies (a), so all of the statements hold. ∎

It is straightforward to show that these new conditions which incorporate $u_{0}$ are valid certificates for transient stability, but this is not enough—if we are to go to the trouble of measuring $u_{0}$, we would like the assurance that this additional information actually leads to better transient stability certificates. With some simple but careful reasoning about the winding partition, we can see that (a) and (b) are less-conservative versions of the antecedents to Theorems 4 and 6, respectively:

Theorem 10** **(Theorem 9 is less conservative than Theorems

4 and 6).

Consider a trajectory $\theta(t)$ of (4), let $\Sigma$ be a cycle basis of the underlying graph, and let $\gamma_{0}=|B^{\mathsf{T}}\theta(0)|$, $\delta_{0}=V(\theta(0))$, and $u_{0}=w_{\Sigma}(\theta(0))$. The following are true:

1. (i)

If the hypothesis of Theorem 4 is satisfied, i.e., if there exist a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ and a set $\Delta(\gamma_{0})\subseteq S\subseteq\Delta(\gamma)$ such that $\delta_{0}<V^{*}(\partial S)$, then the set $S^{\prime}=S\cap\Omega_{u_{0}}$ satisfies condition (a) of Theorem 9. 2. (ii)

If the hypothesis of Theorem 6 is satisfied, i.e., if there exists a vector $\gamma\in[\gamma_{0},\gamma^{*}]$ such that $\delta_{0}<\widehat{V}(\gamma)$, then $\gamma$ satisfies condition (b) of Theorem 9.

Proof.

To prove (i), it is sufficient to show that $V^{*}(\partial S^{\prime})\geq V^{*}(\partial S)$, for which it is sufficient to show that $\partial(S\cap\Omega_{u_{0}})\subseteq\partial S$. Because the winding cells partition $\mathbb{T}^{n}$, each of the sets $\Delta(\gamma)\cap\Omega_{u}$ are disjoint. In fact, because $\gamma<\pi\mathbbold{1}_{m}$, the boundaries of these sets are non-overlapping. Since $S\subseteq\Delta(\gamma)$, we may conclude that $\partial S$ itself is partitioned into non-overlapping pieces $\partial(S\cap\Omega_{u_{0}})$; hence $\partial(S\cap\Omega_{u_{0}})\subseteq\partial S$. Similarly, for (ii) it is sufficient to show that $\widehat{V}(\gamma,u_{0})\geq\widehat{V}(\gamma)$, which we have from Lemma 8. ∎

The initial winding vector provides an additional bit of information about the initial state $\theta(0)$, and like the vector of initial angle differences $\gamma_{0}$, the initial winding vector $u_{0}$ can be inferred from measurements of active power flows. With knowledge of $u_{0}$, we can replace the set-theoretic certificate in Theorem 4 and the MILP certificate in Theorem 6 with the less-conservative conditions (a) and (b) of Theorem 9.

V Quantifying Robustness

An important task in power systems control is understanding how robust an operating point is to disturbances. It is straightforward to study the effects of a particular disturbance using simulation, but simulating a comprehensive set of contingencies (or combinations thereof) is time consuming. Our transient stability certificates can aid with this analysis by quantifying the scale of disturbances to which an operating point is robust.

In this section, we consider a DCMG that is operating at a synchronous state $\theta_{0}$. At time $t=0$, certain model parameters undergo an instantaneous perturbation—nominal injections change, for example, or a drop in nodal voltages or branch admittances occurs. The initial condition $\theta_{0}$ is no longer a synchronous state in the “post-fault” model. If the system is sufficiently resilient, then the post-fault transient will settle back down to a synchronous state, and none of the engineering constraints will be violated in the process—but this is not always the case. We will construct a sufficient condition for post-fault transient stability, based on the scale of the perturbations to model parameters.

Numerical Case Study

Throughout the section, we will illustrate our results using numerical examples from the IEEE-RTS 24-bus test case [12]. We parameterized (4) using branch and bus values from this test case, and we selected the initial voltage angles $\theta_{0}\in\mathbb{T}^{n}$, voltage magnitudes, and nominal power injections by solving for the optimal power flow in MATPOWER. For simplicity, we chose uniform droop coefficients of $10~{}\text{pu}\cdot\text{s}$ and a uniform nominal frequency of $60~{}s^{-1}$. We will refer to this model as the “pre-fault” model. Note that $\theta_{0}$ is a synchronous equilibrium of the pre-fault model (since it solves the active power flow equations with nominal injections), and the initial winding vector is $u_{0}=\mathbbold{0}_{11}$ (corresponding to the winding cell with minimal loop flows).

Code

The code that we used to generate numerical results in this section is publicly available at https://github.com/KevinDalySmith/DCMG-transient-stability. The optimization problems are implemented using the Python interface to Gurobi 9.0, so a local installation of Gurobi and an active license are needed to run it.

V-A Evaluating Post-Fault Transient Stability

We begin by examining how Theorems 6 and 9 apply to the problem of quantifying system robustness. Both of these theorems certify transient stability if the post-fault initial condition (i.e., the pre-fault synchronous state) is sufficiently close to a post-fault synchronous state, as measured by the initial max frequency deviation, $\delta_{0}=V(\theta_{0})$. Then transient stability is certified if there exists any $\gamma\in[\gamma_{0},\gamma^{*}]$ such that $\delta_{0}<\widehat{V}(\gamma)$ or $\delta_{0}<\widehat{V}(\gamma,u_{0})$, as in the following example.

Example 11** (Certificates in the 24-Bus System).**

To illustrate Theorems 6 and 9 in the 24-bus system, we randomly select a set of 100 test points $\Gamma\subset[\mathbbold{0}_{m},\;\gamma^{*}]$ and evaluate $\widehat{V}(\gamma)$ and $\widehat{V}(\gamma,u_{0})$ at each $\gamma\in\Gamma$, with $u_{0}=\mathbbold{0}_{11}$. Consider an arbitrary initial condition with a maximum angle difference $\bar{\gamma}=||B^{\mathsf{T}}\theta(0)||_{\infty}$ and initial frequency deviation $\delta_{0}=V(\theta(0))$. Theorem 6 certifies transient stability of the resulting trajectory if

[TABLE]

Under the additional assumption that $w_{\Sigma}(\theta(0))=\mathbbold{0}_{11}$, Theorem 9 certifies transient stability if

[TABLE]

Figure 3 plots both of these conditions. Each curve plots the right-hand side of the preceding inequalities as a function of $\bar{\gamma}$ (using polytopic relaxations to the sine constraints). For any initial condition corresponding to a point $(\bar{\gamma},\delta_{0})$ below the curve, transient stability is certified.

The curve maximizing $\widehat{V}(\gamma)$ is significantly lower than the curve maximizing $\widehat{V}(\gamma,\mathbbold{0}_{11})$, i.e., the transient stability condition from Theorem 6 is more conservative than that from Theorem 9 (as guaranteed by Theorem 10). This plot makes a strong case for incorporating information from the initial winding vector—without it, only very small angle disturbances are certified in the 24-bus system.

In order to certify transient stability after a fault, we must ensure that the initial post-fault frequency deviation is below the critical threshold. One approach is to follow the procedure of Example 11: generate a plot similar to Figure 3 using the post-fault parameters, and check whether or not $\theta_{0}$ corresponds to a point below the curve. But this approach is cumbersome when considering a large number of contingencies, and it offers little advantage over simulation.

A much more efficient approach, similar to that in [26], is to define the “size” of a general disturbance and establish a threshold below which transient stability is certified in all disturbances that are “smaller” than the threshold. A natural way to define the size of a disturbance is to quantify its effect on the frequency deviation vector from (5). Suppose that $v:\mathbb{T}^{n}\to\mathbb{R}^{n}$ is the frequency deviation vector field defined with the pre-fault model parameters, and similarly, let $\bar{v}$ be the vector field defined with the post-fault parameters. If we can bound the difference $\xi(\theta)=\bar{v}(\theta)-v(\theta)$, then we can bound the solutions to Problems 2 and 3 after the disturbance based on their solutions before the disturbance:

Theorem 12** (Robustness to Parameter Changes).**

Consider the model (4), and let $v:\mathbb{T}^{n}\to\mathbb{R}^{n}$ be the associated frequency deviation vector. Let $\theta_{0}\in\mathbb{T}^{n}$ be a state for which $v(\theta_{0})=0$, i.e., for which all nodal frequencies are identical to $\omega^{*}$. After some perturbation in model parameters, suppose that the new frequency deviation vector is given by $\bar{v}(\theta)=v(\theta)+\xi(\theta)$, and let $\gamma_{0}=|B^{\mathsf{T}}\theta_{0}|$ and $u_{0}=w_{\Sigma}(\theta_{0})$ in the post-fault model. If there exists $\gamma\in[\gamma_{0},\gamma^{*}]$ such that

[TABLE]

where either $S=\Delta(\gamma)$ or $S=\Delta(\gamma)\cap\Omega_{u_{0}}$, then the trajectory of the perturbed model starting from $\theta_{0}$ satisfies statements (i)–(vii) from Theorem 4 with respect to $S$.

Proof.

Let $\Theta\subseteq\partial S$ be the feasible set of Problem 1 evaluated on the perturbed model, i.e., the set of points $\theta\in\partial S$ such that $D^{-1}f(\theta)$ is pointed outward from $S$. Then the solution to Problem 1 (evaluated on the perturbed model) is

[TABLE]

Given the initial condition $\theta_{0}$ to the perturbed model, the initial frequency deviation is $\delta_{0}=||v(\theta_{0})+\xi(\theta_{0})||_{\infty}=||\xi(\theta_{0})||_{\infty}$, so applying (11) and the lower bound on $V^{*}(\partial S)$, we obtain

[TABLE]

Therefore $\gamma$ and $S$ satisfy the hypothesis of Theorem 4 in the perturbed model, and the theorem statements follow. ∎

Condition (11) bounds the scale of the perturbation $\xi(\theta)$. As we will see, in many cases, the left-hand size of the equation is straightforward to compute (or at least upper bound). The right-hand side of (11) is almost identical to either $\widehat{V}(\gamma)$ or $\widehat{V}(\gamma,u_{0})$ (depending on whether $S$ is intersected with the winding cell); the only difference is that the “$D^{-1}f(\theta)$ is pointed outward from $S$” constraint is removed. It is straightforward to obtain a lower bound on $\min_{\theta\in\partial S}||v(\theta)||_{\infty}$ with a minor relaxation to either $\widehat{V}(\gamma)$ or $\widehat{V}(\gamma,u)$: simply remove the $d^{-1}_{s(e)}f_{s(e)}-d_{t(e)}^{-1}f_{t(e)}$ constraints from (6f) and (6g). This lower bound can be used in place of the left-hand side of (11).

In the IEEE 24-bus test case, we use random sampling to identify a point $\gamma\in[\gamma_{0},\gamma^{*}]$ for which $\min_{\theta\in\partial S}||v(\theta)||_{\infty}\geq 0.0435$, with respect to the set $S=\Delta(\gamma)\cap\Omega_{u_{0}}$, $u_{0}=\mathbbold{0}_{11}$. The arc lengths in this particular $\gamma$ range from 18.5 to 22.1 degrees, with a median of 20.3 degrees. Therefore, Theorem 12 guarantees that the IEEE 24-bus steady-state is robust to any perturbations in parameters for which $||\xi(\theta_{0})||_{\infty}+\max_{\theta\in\partial S}||\xi(\theta)||_{\infty}<0.0435$.

In the remaining subsections, we will apply Theorem 12 to particular modes of disturbances: fluctuations in nominal power injections, changes in nodal voltages, and changes in branch admittances.

V-B Perturbations of Nominal Injections

Suppose that the perturbed model is identical to the original model, except the vector of nominal frequency deviations has been shifted to $p^{*}+\Delta p^{*}$. It is then clear from (5) that the vector of frequency deviations suffers the perturbation $\xi(\theta)=D^{-1}\Delta p^{*}$. This quantity is constant, and condition (11) reduces to the condition

[TABLE]

for some $S=\Delta(\gamma)\cap\Omega_{u_{0}}$, with $\gamma\in[\gamma_{0},\gamma^{*}]$. In the IEEE 24-bus test case, a sufficient condition is $||\Delta p^{*}||_{\infty}<0.217~{}\text{pu}$. The median bus in this test case has a nominal injection magnitude of $0.94~{}\text{pu}$, so our certificate guarantees that the system is robust to disturbances of $23\%$ in the nominal injection of this bus. Figure 4 plots the relative tolerance of all buses in the system.

V-C Perturbations of Voltage and Admittance Magnitudes

Next, we consider perturbations to nodal voltage magnitudes and branch admittance magnitudes. Both of these values are encoded in the $\tilde{a}_{ij}$ and $a_{ij}$ parameters, so these perturbations can be represented with perturbations $\Delta\tilde{a}_{ij}$ and $\Delta a_{ij}$. The entries of the corresponding perturbation vector are

[TABLE]

For simplicity, assume that $\Delta\tilde{a}_{ij}\leq 0$ and $\Delta a_{ij}\leq 0$ (i.e., there is a loss in voltage magnitudes or branch admittances). In order to compute $||\xi(\theta_{0})||_{\infty}$, we first compute $\eta_{ij}=\sin(\theta_{i}-\theta_{j}-\phi_{ij})$ using $\theta_{0}$, so that

[TABLE]

Similarly, defining $\bar{\eta}_{ij}=\max\{\sin(\gamma_{\{i,j\}}-\phi_{ij}),\;\sin(\gamma_{\{i,j\}}+\phi_{ij})\}$ as an upper bound on $|\sin(\theta_{i}-\theta_{j}-\phi_{ij})|$ for $\theta\in\Delta(\gamma)$, we can bound

[TABLE]

Then condition (11) is satisfied if the sum of these two quantities is less than $\max_{\theta\in\partial S}||v(\theta)||_{\infty}$. Note that this condition is a set of linear constraints on $\Delta\tilde{a}_{ij}$ and $\Delta a_{ij}$.

For a simple illustrative example, suppose that one particular bus $\ell\in\mathcal{V}$ suffers a loss in voltage magnitude, so that $E_{\ell}\to(1-\alpha)E_{\ell}$ for some $\alpha\in(0,1]$. Then $\Delta\tilde{a}_{\ell j}=-\alpha(2-\alpha)\tilde{a}_{\ell j}$, $\Delta a_{\ell j}=-\alpha a_{\ell j}$, and $\Delta a_{j\ell}=-\alpha a_{j\ell}$ for all $j\in\mathcal{N}(\ell)$, while the remaining perturbations are zero. In the IEEE 24-bus test case, we compute the largest value of $\alpha$ that satisfies the previous equation for each $\ell\in\mathcal{V}$, using 0.0435 for the right-hand side of the bound. These largest $\alpha$ are plotted in Figure 5. The median bus can tolerate a 1% loss of voltage magnitude before $\Delta a_{ij}$ and $\Delta\tilde{a}_{ij}$ violate the above bound. This is much more restrictive than the bound for nominal power injections, which is to be expected, given that we used a conservative upper bound on $\max_{\theta\in\partial S}|\xi_{i}(\theta)|$ instead of the exact value.

VI Conclusion

In this paper, we study transient stability in power networks consisting of droop-controlled inverters and frequency-dependent loads. We extend the notion of transient stability to include not only frequency synchronization but also operating constraints on nodal frequencies, angle differences, power injections, ramping, and storage reserves. To analyze the transients, we introduce a physically-meaningful Lyapunov-like function, and we re-cast the transient stability problem as an optimization problem that admits an efficient relaxation. We show that incorporating information from loop flows (in the form of the winding vector) can make these transient stability certificates less conservative. Finally, we show how these certificates can be used to quantify the size of parameter disturbances to which the system is robust.

The model we use in this paper is, of course, a highly simplified model for frequency dynamics. Nonetheless, we hope that this paper provides a step toward understanding the fundamental behavior of future low-inertia power grids. Extensions of this work may offer rigorous answers to open theoretical questions about these systems. At what scale do the ubiquitous grid-following inverters harm system stability, and how can power engineers use droop-controlled inverters to mitigate this effect? How do legacy high-inertia generators affect transient behavior in power networks that are dominated by inverters? Future research may enrich this work with models of different types of generators, to better-understand frequency dynamics as power grids transition to low-inertia power sources.

Appendix A Proof of Lemmas 5 and 8

Proof of Lemma 5.

Let $\theta$ be the minimizing argument of Problem 1, let $y_{e}=d_{\rm cc}(\theta_{i},\theta_{j})$ be the counterclockwise angle difference across each branch $e\in\mathcal{E}$, and define $f$, $\eta^{+}$, and $\eta^{-}$ according to (6b)–(6d). Because $\theta\in\partial\Delta(\gamma)$, (6e) holds because $\theta\in\mathrm{cl}(\Delta(\gamma))$, and there exists some edge $e\in\mathcal{E}$ such that $y_{e}=\pm\gamma_{e}$. In the case where $y_{e}=\gamma_{e}$, let $z_{e}^{+}=1$, so the two linear constraints in (6f) activate. The first constraint is just $y_{e}=\gamma_{e}$, which we have assumed is true. The second constraint, $d_{s(e)}^{-1}f_{s(e)}-d_{t(e)}^{-1}f_{t(e)}\geq 0$, holds because $D^{-1}f(\theta)=D^{-1}f$ points outward from $\Delta(\gamma)$. Thus (6f) is satisfied, and (6g) and (6h) are satisfied by setting all other entries of $z^{+}$ and $z^{-}$ equal to zero. In the case where $y_{e}=-\gamma_{e}$, we employ a complementary argument with $z^{-}_{e}=1$. In both cases, all constraints are satisfied, so $(f,y,\eta^{+},\eta^{-},z^{+},z^{-})$ is feasible. Finally, the cost function (6a) is equal to $V(\theta)$ (the cost function of Problem 1) when evaluated at $\theta$. Hence $\widehat{V}(\gamma)\leq V^{*}(\partial\Delta(\gamma))$.

To see that equality holds in the tree case, let us consider the argmin $(f,y,\eta,z^{+},z^{-})$ of Problem 2. Because $G$ is a tree and $|y|\leq\gamma$, there always exists $\theta\in\Delta(\gamma)$ for which $y_{e}=d_{\rm cc}(\theta_{s(e)},\theta_{t(e)})$ for all $e\in\mathcal{E}$. Then (6e)–(6h) ensure that $\theta\in\partial\Delta(\gamma)$, and (6a)–(6d) ensure that the cost function of Problem 2 is identical to $V(\theta)$. Hence $V^{*}(\partial\Delta(\gamma))\leq\widehat{V}(\gamma)$ if $G$ is acyclic. ∎

The proof of Lemma 8 proceeds similarly, but it uses some properties of the winding partition of the $n$-torus. Readers unfamiliar with the winding partition are encouraged to read IV-A and glance at [17]. The property that we need is the following lemma, which shows how the winding partition relates to the boundaries of phase-cohesive sets:

Lemma 13**.**

Let $\Omega_{u}$ be a winding cell, and let $\gamma\in(0,\pi\mathbbold{1}_{m})$. Consider the set $S=\Delta(\gamma)\cap\Omega_{u}$. Then $\partial S=\partial\Delta(\gamma)\cap\Omega_{u}$.

Proof.

We first argue that $\mathrm{cl}(\Delta(\gamma))\cap\partial\Omega_{u}=\emptyset$. If $\theta\in\mathrm{cl}(\Delta(\gamma))$, then $|\theta_{i}-\theta_{j}|<\pi$ for all $\{i,j\}\in\mathcal{E}$ (by the assumption that $\gamma<\pi\mathbbold{1}_{m}$). Then there is a neighborhood around $\theta$ within which the winding numbers around each cycle do not change, so $\theta$ belongs to the interior of a winding cell.

First, using elementary properties of topology, we have

[TABLE]

Because $\partial S\subseteq\mathrm{cl}(\Delta(\gamma))$, we have that $\partial S\cap\partial\Omega_{u}=\emptyset$. Therefore, from the elementary property $\partial S\subseteq\partial\Delta(\gamma)\cup\partial\Omega_{u}$, we obtain the result $\partial S\subseteq\Delta(\gamma)$. Furthermore, because $\partial S\subseteq\mathrm{cl}(\Omega_{u})$ but $\partial S\cap\partial\Omega_{u}=\emptyset$, we have that $\partial S\subseteq\mathrm{int}(\Omega_{u})\subseteq\Omega_{u}$. Hence $\partial S\subseteq\partial\Delta(\gamma)\cap\Omega_{u}$.

To show equality, we invoke the winding partition to write

[TABLE]

For $v\neq u$, the sets $\Omega_{u}$ and $\partial(\Delta(\gamma)\cap\Omega_{v})$ are disjoint, since we have shown that $\partial(\Delta(\gamma)\cap\Omega_{v})\subseteq\Omega_{v}$. Thus, intersecting both sides of the equation with $\Omega_{u}$, we obtain $\partial\Delta(\gamma)\cap\Omega_{u}\subseteq\partial(\Delta(\gamma)\cap\Omega_{u})=\partial S$. This completes the proof. ∎

We now prove Lemma 8.

Proof of Lemma 8.

The statement that $\widehat{V}(\gamma)\leq\widehat{V}(\gamma,u)$ is obvious, since Problem 2 is a relaxation of Problem 3.

To show that $\widehat{V}(\gamma,u)\leq V^{*}(\partial S)$, let $\theta$ be the minimizing argument of Problem 1. As in the proof of Lemma 5, select the values of the decision variables $(f,y,\eta,z^{+},z^{-})$ accordingly to satisfy (6b)–(6h), thereby ensuring that the cost function (6a) is equal to $V(\theta)$. Because $\theta\in\Omega_{u}$, [17, Theorem 3.5] guarantees the existence of $x\in\mathbbold{1}_{n}^{\mathsf{T}}$ such that $y=B^{\mathsf{T}}x+2\pi C_{\Sigma}^{\dagger}u$. Multiplying across by $C_{\Sigma}$, we obtain

[TABLE]

We have performed two simplifications in this equation. First, $\operatorname{\mathrm{Img}}(C_{\Sigma}^{\mathsf{T}})$ is the cycle space, which is identical to $\ker(B)$ [6, Theorem 9.5], so the $C_{\Sigma}B^{\mathsf{T}}x$ term vanishes. Second, the summation structure in (7) implies that $u\in\operatorname{\mathrm{Img}}(C_{\Sigma})$, so the orthogonal projection matrix $C_{\Sigma}C_{\Sigma}^{\dagger}$ has no effect on $u$. Thus $(f,y,\eta)$ satisfies all of the constraints of Problem 3, so $\widehat{V}(\gamma,u)\leq V^{*}(\partial S)$.

Next, we will show that $V^{*}(\partial S)\leq\widehat{V}(\gamma,u)$. Let $(f,y,\eta,z^{+},z^{-})$ be the minimizing argument of Problem 3. Consider the equation $y=B^{\mathsf{T}}x+2\pi C_{\Sigma}^{\dagger}u$. Note that $\operatorname{\mathrm{Img}}(B^{\mathsf{T}})=(\operatorname{\mathrm{Img}}(C_{\Sigma}^{\dagger}))^{\bot}$, so we can decompose $y=y_{1}+y_{2}$ with $y_{1}\in\operatorname{\mathrm{Img}}(B^{\mathsf{T}})$ and $y_{2}=\operatorname{\mathrm{Img}}(C_{\Sigma}^{\dagger})$, and the equation can be split into $y_{1}=B^{\mathsf{T}}x$ and $y_{2}=2\pi C_{\Sigma}^{\dagger}u$. The first equation has a unique solution $x\in\mathbbold{1}_{n}^{\bot}$, while the second equation is true because $C_{\Sigma}y=2\pi u$, so there is a unique point $x\in\mathbbold{1}_{n}^{\bot}$ that satisfies $y=B^{\mathsf{T}}x+2\pi C_{\Sigma}^{\dagger}u$. It follows from [17, Theorem 3.5] that there exists $\theta\in\Omega_{u}$ such that $y_{e}=d_{\rm cc}(\theta_{s(e)},\theta_{t(e)})$ for all $e\in\mathcal{E}$. From the remaining constraints in Problem 3, we can see that $\theta\in\partial\Delta(\gamma)\cap\Omega_{u}$, so it follows from Lemma 13 that $\theta\in\partial S$. Furthermore, the constraints imply that the velocity vector is pointed outward. Thus $\theta$ is within the feasible set of Problem 1, and the identical values of the cost functions imply that $V^{*}(\partial S)\leq\widehat{V}(\gamma,u)$. ∎