# The Largest Empty Sphere Problem in 3D Hollowed Point Clouds

Netzer Moriya

##  1 Abstract

We introduce a new approach for the adaptation of the Maximal Internal
Envelope method, extended to address the Largest Empty Sphere problem within
unstructured 3D point clouds. We explore the identification of the Largest
Empty Sphere by computing Convex Hull vertices and employing a Voidness Score
based on Minimal Distance Scoring for optimal segment selection. The
integration of Delaunay triangulation and Voronoi diagrams facilitates the
initial identification of potential Largest Empty Sphere candidates. Our
analysis reveals the method’s efficacy and efficiency, often locating the
Largest Empty Sphere in initial computational stages, suggesting a lower
complexity than initially projected.

##  2 Introduction

The general Largest Empty Sphere Problem (LESP), addresses the problem of
finding a hypersphere of largest radius in d-dimensional space whose interior
does not overlap with any given obstacles. In the case of a 3D point cloud,
the problem becomes finding the largest empty sphere within the point cloud.
This problem has broad applications across various disciplines, including 3D
facility location, molecular modeling, and robotics in three-dimensional
environments.

Historically, the LESP has been approached from various perspectives, each
contributing unique insights into the problem’s complexity. Aistleitner [1]
explored the concept in higher dimensions, examining the largest empty box and
its volume behavior in relation to dimensionality (see also [2]). In contrast,
Shi’s work [3] focused on a convex optimization approach for identifying the
largest empty sphere. Toussaint [4] proposed an algorithm for the largest
empty circle, integrating locational constraints, while Ryshkov [5] discussed
the empty sphere in the context of lattice coverings, adding a historical
perspective to the problem.

A common solution technique involves generating a mesh of evenly distributed
points throughout the space, then applying a nearest neighbors algorithm to
determine the furthest mesh point from any given point $P_{i}\in P$, thereby
identifying the center of the largest void. However, this approach can be
exteremly computationally demanding in higher dimensions. An alternative
strategy involves Delaunay triangulation, assessing the volume of each formed
tetrahedron. These tetrahedrons represent potential voids, with the largest
among them potentially constituting the largest empty sphere. A different
approach employs the Voronoi diagram on a sphere, combinatorially akin to a 3D
convex hull. In this model, each hull face correlates to a Voronoi vertex,
linked to an empty circle, leading to an efficient $O(n\log n)$ algorithm [6].

The LES problem in shell-like or hollowed structures presents a more intricate
challenge, demanding specialized algorithms attuned to their unique spatial
configurations [7]. Conventional methods such as Voronoi diagrams and Delaunay
triangulation by themselves might not be suffice in such cases. Instead,
stochastic search approaches [8] and evolutionary algorithms have proven
effective in handling multidimensional spaces, offering advantages in
execution time and memory usage [9]. Within the context of metallic hollow
spheres, extensive research has been conducted on uniaxial deformation
behavior under significant deformations, with computational models developed
for these scenarios [10]. Additionally, finite elastostatics within
pressurized hollow spheres for compressible materials have been thoroughly
investigated, yielding explicit closed-form solutions for deformation and
stress fields [11]. Further research has provided closed-form solutions for
hollowed sphere models with specific materials under isotropic loadings,
offering valuable reference solutions and critical evaluations of macroscopic
criteria [12].

This paper aims to extend these methodologies to identify the largest empty
sphere within shell-like or hollowed 3D point clouds. We introduce a semi-
heuristic approach, commencing with Delaunay-Voronoi based void detection,
followed by convex hull construction for guiding sphere selection. Our
methodology ensures no cloud members are enclosed within the identified
sphere, addressing the unique challenges posed by these structures.

##  3 Problem Statement

Given a set $P$ of $n$ points representing a hollowed structure in the
Euclidean space, we seek to identify the largest sphere within the hollow
region that does not enclose any point in $P$. Unlike traditional LESP
formulations, the hollowed structure presents a unique challenge as the
largest empty sphere is expected to be centered within the void of the
structure, rather than merely avoiding the outer points of $P$. Formally, the
problem seeks to maximize the radius $r$ of a sphere with center $c$, located
within the hollow region, ensuring that for all points $p\in P$, the distance
from $c$ to $p$ is at least $r$.

###  3.1 Mathematical Formulation

The problem is formalized in the Euclidean space $\mathbb{R}^{3}$. Considering
a finite set of points $P=\\{p_{1},p_{2},\ldots,p_{n}\\}$ that outlines a
hollowed structure, the aim is to determine the center $c$ and radius $r$ of a
sphere $S\in S^{LES}$, where $S^{LES}$ denotes the set of potential largest
empty spheres in the 3D point cloud, satisfying the following conditions:

  * •

The center $c$ of the sphere $S$ lies within the void of the structure formed
by $P$, adhering to the spatial constraints imposed by the hollowed
configuration.

  * •

For every $p\in P$, the condition $\|c-p\|\geq r$ holds, where $\|\cdot\|$
denotes the Euclidean norm. This ensures that $S$ is external to all points in
$P$.

  * •

The objective is to maximize $r$, subject to the aforementioned conditions,
representing the largest feasible sphere within the hollow.

Additionally, we introduce an auxiliary function
$\phi:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}$ defined as
$\phi(c,r)=\min_{p\in P}\|c-p\|-r$. The function $\phi$ quantifies the minimum
distance from the sphere boundary to any point in $P$, subtracting the radius
$r$. The LES problem then becomes finding
$(c^{*},r^{*})=\arg\max_{(c,r)}\phi(c,r)$, subject to $\phi(c,r)\geq 0$ and
other geometric constraints imposed by the hollow structure. To facilitate
this, the Minimal Distance Score (MDS) is employed to assign a ’Void Score’ to
various regions within the 3D point cloud. This score aids in pinpointing
areas that are in closer proximity to the internal cavity of the cloud.

##  4 Proposed Methodology

We propose a new approach to identify the largest empty sphere (LES) within a
3D point cloud, an extension of the Maximal Internal Envelope (MIE) method
outlined in [13] from 2D to 3D contexts. This adaptation addresses the unique
challenges posed by 3D structures, particularly when analyzing voids or
cavities.

Initially, the process involves computing the Convex Hull vertices
$P^{CH}_{i}$ to encapsulate the point cloud. The subsequent phase includes a
systematic selection of optimal line segments $L^{best}_{ij}$, which are
strategically chosen based on their likelihood of intersecting with internal
void regions by assigning a ’Voidness Score’ to each segment and choosing the
$L^{best}_{ij}$. The approach then integrates Delaunay triangulation with
Voronoi diagrams to identify initial points (referred to as DV points) for
constructing the sets of both MIE points and largest empty spheres.

This process is iteratively developed using a pseudo-recursive method,
enhancing higher-order MIE points detection and potential LES.

Each iteration refines the set of MIE points, facilitating the identification
of additional candidate points. Starting with a hypothetical sphere centered
at a DV point on segment $L^{best}_{ij}$, this sphere incrementally expands
until it first intersects with points in the cloud. The sphere’s center then
moves along $L^{best}_{ij}$, incrementally expanding in predefined step sizes
(controlled by a parameter $3\leq k\leq n$). Points first intersected by the
expanding sphere are aggregated into set $S^{LES}$ and recorded as candidate
LES spheres, alongside the sphere centers and corresponding radii. This
procedure is methodically repeated for all segments $L_{ij}$ associated with
$P^{CH}_{i}$, thereby ensuring thorough exploration and inclusion of
significant spheres within the void’s interior. As the simulation converges to
its final stage, the largest sphere, $\max_{\text{sphere}\in
S^{LES}}\text{radius}(\text{sphere})$, is then selected as the LES in the 3D
point cloud.

We applied this method to two distinct point set configurations: a single
sphere with random dispersion along its circumference, and a pair of congruent
spheres with a deliberately engineered internal void. In both cases, our
method successfully revealed the internal set of potential largest empty
spheres, identifying the sphere with the maximal radius as the LES.

###  4.1 A Short Description of the MIE method

Here we relate to the principle of the method, described in details in [13]
for the 2D case. In a following section, we extend the method to the 3D
scenario and describe the inclusion of the LES identification as mentioned
above.

The proposed methodology for analyzing voids in a three-dimensional point
cloud employs a two-stage process, integrating geometric and computational
techniques.

Stage 1: Voids Localization begins with convex hull analysis to identify
principal points of the cloud. Using Delaunay Triangulation, Voronoi Diagrams,
and the Minimal Distance Scoring (MDS) technique, segments likely to intersect
the void are determined. We then identify first-order Maximal Internal
Envelope (MIE) points based on proximity and angular relations to these
segments.

Stage 2: Iterative Refinement builds upon the MIE points, enhancing the
structural definition of the void iteratively. This computationally intensive
approach yields a detailed representation of the void’s geometry.

####  4.1.1 Key Steps and Equations

Minimal Distance Scoring (MDS) is central to our methodology. The MDS of a
segment $L_{ij}$ is given by:

$\text{MDS}(L_{ij})=\frac{1}{|P\setminus P_{CH}|}\sum_{p_{k}\in P\setminus
P_{CH}}d(L_{ij},p_{k})$ (1)

Here, $d(L_{ij},p_{k})$ represents the distance from a line segment $L_{ij}$
to a point $p_{k}$ in the point set $P$ excluding the points in the Convex
Hull $P_{CH}$. This scoring mechanism is fundamental in identifying segments
within the point cloud that are likely candidates for forming the boundary of
the largest void.

In conjunction with MDS, we introduce an auxiliary function
$\phi:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}$, defined as
$\phi(c,r)=\min_{p\in P}\|c-p\|-r$. This function $\phi$ quantifies the
minimum distance from the boundary of a sphere with center $c$ and radius $r$
to any point in $P$, ensuring the sphere does not enclose any point in $P$.
Therefore, the LES problem is reformulated as an optimization problem where we
seek to maximize $r$ while ensuring $\phi(c,r)\geq 0$. The MDS aids in this
process by providing a ’Voidness Score’ to segments, helping in the selection
of the optimal center $c$ and radius $r$ for the largest empty sphere, while
adhering to the geometric constraints of the hollow structure.

Stages in Voids Identification include:

  1. 1.

Convex Hull Identification.

  2. 2.

Segment Analysis for intersecting the void using MDS.

  3. 3.

Delaunay Triangulation and Voronoi Diagrams to locate the optimal point on the
best segment.

Void Shape Construction uses an expanding polygon initiated from the best
segment identified by the MDS algorithm. This involves:

  1. 1.

Selecting vertices based on proximity to $L^{best}_{ij}$.

  2. 2.

Iterative expansion and adaptation of the polygon’s shape influenced by points
in $P$.

####  4.1.2 Complexity Analysis

The complexity analysis includes:

  * •

Convex Hull Calculation: $O(n\log n)$.

  * •

Pair Generation of Convex Hull Vertices: $O\left(\frac{n^{2}}{100}\right)$.

  * •

Overall Complexity: $O(n\log n)+4\times O\left(\frac{n^{2}}{100}\right)$.

This methodology effectively delineates the shape characteristics of voids in
point clouds, balancing detailed geometric analysis with computational
efficiency.

###  4.2 Identification of the Largest Empty Sphere in the Shell-Like 3D Point
Cloud

Given the above methodology for MIE points identification, the expanding
spheres intersecting with the MIE points, are aggregated into set $S^{LES}$
and recorded as candidate LES spheres, alongside the sphere centers and
corresponding radii. This procedure is methodically repeated for all segments
$L_{ij}$ associated with $P^{CH}_{i}$, thereby ensuring thorough exploration
and inclusion of significant spheres within the void’s interior. Further
orders of the calculations, based on higher order segments [13] are also
considered and added to the set $S^{LES}$. As the simulation converges to its
final stage, the largest sphere, $\max_{\text{sphere}\in
S^{LES}}\text{radius}(\text{sphere})$, is then selected as the LES in the 3D
point cloud.

The following schematic algorithm 1 detailing the computational steps for void
identification and the estimation of the LES in the cloud.

Table 1: Algorithm for the LES estimation Step | Action
---|---
1: | Identify Convex Hull $P^{CH}$ of $P$
2: | Select a starting point $P^{CH}_{i}$ from $P^{CH}$
3: | For each $P^{CH}_{i}$, create and analyze segments $L_{ij}$
4: | Calculate MDS($L_{ij}$) for each segment
5: | Select $L^{best}_{ij}$ with minimum MDS value
6: | Apply Delaunay Triangulation and Voronoi Diagrams to select DV point on $L^{best}_{ij}$
7: | Along each segment $L^{best}_{ij}$, define a dynamic Hypothetical Sphere that expands from the DV point.
8: | Identify Maximal Internal Envelope (MIE) points from $P$ that intersect with the expanding circle.
9: | Define the Hypothetical Sphere as a LES candidate
10: | Incorporate the Hypothetical Sphere into the subset $S^{LES}$.
11: | Update the set $S^{LES}$ and repeat the process for the new segments formed.
12: | Terminate the iterative process when the addition of new MIE points ceases.
13: | Select $\max_{\text{sphere}\in S^{LES}}\text{radius}(\text{sphere})$ as the LES in the cloud.

Figure 1: Schematic Algorithm for the LES estimation.

##  5 Results - Largest Empty Sphere Selection

In addressing the challenge of identifying the Largest Empty Sphere (LES)
within an unstructured 3D point cloud, we have adapted and applied the Maximal
Internal Envelope (MIE) method. This adaptation involved extending the
original MIE approach, as detailed in [13], to incorporate a third dimension.
The process entailed constructing a collection of potential candidates for the
largest empty sphere. Among these candidates, the sphere that most accurately
represents the LES was then determined and selected as the optimal solution.

###  5.1 Implementation and Analytical Insights

We applied the adapted MIE method to two distinct 3D point cloud
configurations: a shell-like sphere and a pair of offset congruent spheres.
The Convex Hull vertices ($P^{CH}_{i}$) for each configuration were computed,
revealing the geometric structure of the point clouds (Figure 2).

For each configuration, optimal line segments ($L^{best}_{ij}$) were
identified based on a ’Voidness Score’, derived from Minimal Distance Scoring
(MDS). Delaunay triangulation and Voronoi diagrams facilitated the
identification of initial points for LES estimation. The dynamic projection of
hypothetical spheres from these points and their interaction with the cloud
were analyzed to determine the potential LES candidates (Figures 3 and 4).

Figure 2: The 3D point distributions of a shell-like sphere (left) and the two
spheres construction with an internal cavity (right).

The process involved projecting a hypothetical, dynamic sphere from the
Delaunay-Voronoi (DV) point along the selected segment. Points intersecting
with the radius changing sphere, as shown in Figure 4, were included in the
subset $C^{MIE}\subseteq P$ with all potential largest empty spheres include
in $S^{LES}$. These points, determaning the radii of the potential Largest
Empty Spheres, were then subject to maximal radius analysis and upon
identifying the largest radius, the corresponding sphere was identified as the
best estimation for the LES.

Figure 3: Example of best segments selection based on MDS for the two clouds.
Every segment is connecting CH pairs of the respective 3D points distribution.

Figure 4: Example of MIE points identification based on best segment and VD
starting point on segment for the LES selection.

###  5.2 Maximal Radius Analysis and LES Selection

The collection of potential LES candidates ($S^{LES}$) was subjected to a
maximal radius analysis. The largest sphere from $S^{LES}$ was identified as
the LES (Figure 5).

Figure 5: The best sphere entitled LES for the shell-like sphere (left) and
the two-spheres configuration (right)

###  5.3 Computational Efficiency and Observations

Our analysis revealed the MIE method’s efficacy and efficiency in identifying
the LES within diverse 3D point cloud configurations. Interestingly, the LES
was often identified in the initial computation phase, indicating a lower
computational complexity than initially projected.

Key Insights:

  * •

Efficiency and Robustness: The method of Maximal Internal Envelope has been
proved both eficient and robust method to identify that Largest Empty Sphere
in cases where cavity is enclosed in a 3D point cloud.

  * •

Fast Convergence: In all our observations, over a vast number of simulations
obtained with numerious different 3D point cloud configurations, the LES was
identified during the first order of the calculationd, i.e., there was no need
for higher order iterations as the largest empty sphere was almost imediatly
identified during the CH stage.

  * •

Computational Complexity: The overall expected complexity analysis is
substantially lower than the $O(n\log n)+4\times
O\left(\frac{n^{2}}{100}\right)$ deduced for the MIE case. For the LES
problem, our assumption is that only the first term should be considered..

##  6 Conclusions

The Largest Empty Sphere Problem involving 3D hollowed structures, is an
important extension of the classical problem, representing key concepts in
spatial analysis and optimization of unstructured point clouds in 3D.

This study demonstrates the adaptation of the Maximal Internal Envelope (MIE)
methodology, explicitly tailored to tackle the challenge of identifying the
Largest Empty Sphere (LES) within unstructured 3D point clouds. The approach
described in this paper involves the computation of Convex Hull vertices,
followed by the application of a Voidness Score derived from the Minimal
Distance Scoring (MDS) algorithm, which is instrumental in the optimal
selection of segments. Furthermore, the incorporation of Delaunay
triangulation and Voronoi diagrams plays a pivotal role in the preliminary
identification of prospective LES candidates. Empirical analysis of this
method, proved efficient in identifying the LES within diverse 3D point cloud
configurations. Notably, the LES was often identified in the initial
computational phase, indicating a lower computational complexity than
previously estimated, prior to the initial simulation experiments. Future work
may focus on enhancing the algorithm’s efficiency and exploring its
applicability to larger, more intricate datasets in 3D facility location,
molecular modeling, and robotics in three-dimensional environments and related
fields.

## Declarations

All data-related information and coding scripts discussed in the results
section are available from the corresponding author upon request.

## References

  * [1] Christoph Aistleitner, Aicke Hinrichs, and Daniel Rudolf, “On the size of the largest empty box amidst a point set,” Discrete Applied Mathematics, vol. 230, pp. 146-150, 2017. DOI: https://doi.org/10.1016/j.dam.2017.06.008. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0166218X17302949
  * [2] Chazelle, B. Computing the Largest Empty Rectangle. in SIAM J. Comput. 1986
  * [3] Jianming Shi and Yamamoto Yoshitsugu, “A D.C. Approach to the Largest Empty Sphere Problem in Higher Dimension,” in Proceedings of the 1996 International Conference, 1996. [Online]. Available: https://api.semanticscholar.org/CorpusID:118666703
  * [4] Toussaint, G.T. Computing largest empty circles with location constraints. International Journal of Computer and Information Sciences 12, 347–358 (1983). https://doi.org/10.1007/BF01008046
  * [5] S. S. Ryshkov and R. M. Erdahl, “The Empty Sphere Part II,” Canadian Journal of Mathematics, vol. 40, no. 5, pp. 1058–1073, 1988. DOI: 10.4153/CJM-1988-043-5.
  * [6] Hyeon-Suk Na, Chung-Nim Lee, and Otfried Cheong, “Voronoi diagrams on the sphere,” Computational Geometry, vol. 23, no. 2, pp. 183-194, 2002. DOI: https://doi.org/10.1016/S0925-7721(02)00077-9. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0925772102000779
  * [7] Anirban Nag, Vikranth Reddy, Saurav Agarwal, and Sandipan Bandyopadhyay, “Identifying Singularity-Free Spheres in the Position Workspace of Semi-regular Stewart Platform Manipulators,” in Advances in Robot Kinematics 2016, Jadran Lenarčič and Jean-Pierre Merlet, Eds. Cham: Springer International Publishing, 2018, pp. 421–430. https://doi.org/10.1007/978-3-319-56802-7_44
  * [8] M. Hughes, M. Goerigk, and M. Wright, “A Largest Empty Hypersphere Metaheuristic for Robust Optimisation with Implementation Uncertainty,” Computers & Operations Research, vol. 103, pp. 64-80, 2019. DOI: [https://doi.org/10.1016/j.cor.2018.10.013](https://doi.org/10.1016/j.cor.2018.10.013). URL: [https://www.sciencedirect.com/science/article/pii/S0305054818302727](https://www.sciencedirect.com/science/article/pii/S0305054818302727).
  * [9] Jong-Seok Lee, Taeg Sang Cho, Jiye Lee, Myung-Kee Jang, Tae-Kwang Jang, Dongkyun Nam, and Cheol Hoon Park, “A stochastic search approach for the multidimensional largest empty sphere problem,” in Proceedings of the 2004 International Conference, 2004. [Online]. Available: https://api.semanticscholar.org/CorpusID:2790670
  * [10] Marco Speich, Wolfgang Rimkus, Markus Merkel, and Andreas Öchsner, “Large Deformation of Metallic Hollow Spheres,” Materials Science Forum, vol. 623, pp. 105 - 117, 2009. [Online]. Available: https://api.semanticscholar.org/CorpusID:136856535
  * [11] Rohan Abeyaratne and Cornelius O. Horgan, “The pressurized hollow sphere problem in finite elastostatics for a class of compressible materials,” International Journal of Solids and Structures, vol. 20, pp. 715-723, 1984. [Online]. Available: https://api.semanticscholar.org/CorpusID:121867436
  * [12] Ph. Thoré, Franck Pastor, J. Pastor, and Djimédo Kondo, “Closed-form solutions for the hollow sphere model with Coulomb and Drucker–Prager materials under isotropic loadings,” Comptes Rendus Mecanique, vol. 337, pp. 260-267, 2009. [Online]. Available: https://api.semanticscholar.org/CorpusID:123367929
  * [13] Netzer Moriya, “Void Shape Identification in a 2D Point Distribution,” 2023, eprint: arXiv:2312.17533 [cs.CG]. DOI: https://doi.org/10.48550/arXiv.2312.17533.


# Ambient Forcing: Sampling Local Perturbations in Constrained Phase Spaces

Anna Büttner<EMAIL_ADDRESS>Potsdam Institute for Climate Impact
Research  Jürgen Kurths  Potsdam Institute for Climate Impact Research  Frank
Hellmann<EMAIL_ADDRESS>Potsdam Institute for Climate Impact
Research

## Abstract

Ambient Forcing is a novel method to sample random states from manifolds of
differential-algebraic equations (DAE). These states can represent local
perturbations of nodes in power systems with loads, which introduces
constraints into the system. These states must be valid initial conditions to
the DAE, meaning that they fulfill the algebraic equations. Additionally,
these states should represent perturbations of individual variables in the
power grid, such as a perturbation of the voltage at a load. These initial
states enable the calculation of probabilistic stability measures of power
systems with loads, which was not yet possible, but is important as these
measures have become a crucial tool in studying power systems. To verify that
these perturbations are network local, i.e. that the initial perturbation only
targets a single node in the power grid, a new measure, the spreadability,
related to the closeness centrality [1], is presented. The spreadability is
evaluated for an ensemble of typical power grids. The ensemble depicts a set
of future power grids where consumers, as well as producers, are connected to
the grid via inverters. For this power grid ensemble, we additionally
calculate the basin stability [2] as well as the survivability [3], two
probabilistic measures which provide statements about asymptotic and transient
stability. We also revisit the topological classes, introduced in [4], that
have been shown to predict the basin stability of grids and explore if they
still hold for grids with constraints and voltage dynamics. We find that the
degree of the nodes is a better predictor than the topological classes for our
ensemble. Finally, ambient forcing is applied to calculate probabilistic
stability measures of the IEEE 96 test case [5].

##  1 Introduction

Currently, power grids all over the world are reconstructed to reduce CO2
emissions. Many sectors, like the mobility sector, are electrified which leads
to an increased electricity demand. There are more loads in the grid whose
energy demand needs to be fulfilled to assure a stable operation. Therefore,
the power grid, which is already a critical infrastructure, becomes even more
important and stability has to be assured at all times.

Realistic power grids have to face a vast amount of possible scenarios, where
the system is perturbed from their stable operation. As these perturbations
are not necessarily small an investigation beyond linear stability analysis is
necessary. Probabilistic methods have proven to be an important tool for
studying the stability of power systems. The Basin Stability method, in
particular, should be emphasized in this context. It has provided novel
insights into multi-stability and resilience for complex systems, but also for
power grids in particular. [2, 6]. In power grids it is especially interesting
to perturb single nodes in the network to find especially vulnerable nodes and
to deduce topologies or network motifs which endanger the stability of the
network [6, 4].

Typically, power system stability is divided into several categories, which
traditionally include frequency and voltage stability [7]. The stability of
the voltages in power grids is a problem that is likely to become even more
relevant in the future [8] especially due to the increased share of
renewables. As renewable energy sources are connected to the grid via
inverters, the IEEE has updated their definitions for power grid stability to
include converter-driven-stability[7], which is mostly concerning inverters in
the power grid. This paper will focus on investigating three stability
categories, frequency, voltage, and inverter-based stability using
probabilistic stability assessments.

Although it is often advocated to increase the share of grid forming
inverters, which can maintain grid stability even when no synchronous machines
are present, grid following inverters, which introduce constraints into the
system and can not sustain grid stability by themselves, are still the
majority of inverters today. Therefore, grid following inverters have to be
included in a stability assessment. In the case of power grids, the
mathematical problem then naturally appears in the form of differential-
algebraic equations (DAEs).

It should be noted that power grids are always subjected to constraints due to
Kirchhoff’s circuit laws, which ultimately result from fundamental physical
laws. The mesh rule is a special case of Faraday’s law without an exterior
magnetic field and the nodal rule is a result of charge conservation [9]. Thus
these laws must never be broken and the system reacts infinitely fast to
fulfill them. There is also another type of constraint that is less rigid but
also important; it arises in the network due to the separation of time scales.

In case of a perturbation in the power grid the demand can react in a variety
of ways, but typically does so much faster than the time scales which we look
at and in such a way to recover the energy drawn at the node. Hence we will
model the loads as constraints.

Probabilistic stability measures cannot simply be calculated for DAE systems
since their phase space is constrained and random valid initial states, i. e.
states that fulfill the constraints, cannot be sampled arbitrarily. The
process of calculating valid initial conditions for DAEs is called DAE
initialization and is a well researched topic in numerical mathematics and
considered a challenging numerical problem. Established DAE initialization
methods, like Brown’s method [10], can not be used to sample single node
perturbations, since only a part of the system is initialized, and it is not
possible to regulate what happens to the other variables during the
initialization process. Thus it can not be controlled if a perturbation is
only targeted to a single node.

In this paper, we will develop a novel method that allows us to sample valid,
locally perturbed initial conditions for power grids in DAE form. We will use
this method to calculate the Single Node Basin Stability and survivability of
an ensemble of power grids and of the IEEE 96 test case. Finally, we will
revisit network measures from previous studies and check if they are still a
good predictor for the probabilistic stability measures.

##  2 Probabilistic Stability Measures

###  2.1 Single Node Basin Stability

The probability that a power system remains in the operating point and does
not run to one of its other attractors depends on its stability against non-
small disturbances [2]. The analysis presented by Menck et al. [2] introduced
Basin Stability $\langle{\beta}\rangle$ as a non-local and nonlinear measure
to quantify the stability of a system. The basin of attraction $\mathcal{B}$
of a state $x^{*}$ is the set of all initial conditions leading back to
$x^{*}$. Analytical computation of basins is challenging and becomes
increasingly difficult for complex or higher dimensional systems. Therefore,
the authors [2] rely on a single property of the basin $\mathcal{B}$ \- its
volume. The volume of the basin $\mathcal{B}$ is easier to determine and
contains information about the probability of returning to $x^{*}$. Thus, it
can be interpreted as a measure of stability [2].

The Basin Stability $\langle{\beta}\rangle$ of a state $x^{*}$ in a dynamical
system is usually calculated by a numerical Monte Carlo method [6]. Disturbed
states $x_{0}$ are randomly drawn from the phase space $\mathcal{S}$. The
perturbed system in state $x_{0}$ is then integrated, and it is checked,
whether the system returns to $x^{*}$ or not. The Basin Stability
$\langle{\beta}\rangle$ is then defined as:

$\displaystyle\langle{\beta}\rangle=\frac{V_{x^{*}}}{V}$ (1)

where $V$ is the number of all randomly drawn states and $V_{x^{*}}$ is the
number of states that returned to $x^{*}$. In the case of power grids, Basin
Stability can be used to study the stability of the system after the
disturbance of a single node [6].

In Single Node Basin Stability, only the variables of a single node are
perturbed, while the rest of the nodes remain at the operating point.
$\langle{\beta_{a}}\rangle=[0,1]$ is then the probability that the system
returns to the operating point after node $a$ is perturbed.

###  2.2 Single Node Survivability

Supplementary to the basin stability, which examines the asymptotic behavior
of a system after a perturbation, the transient behavior can be analyzed and
whether a system will remain in a desirable regime, for some time $t$, after a
perturbation. Hellmann et. al. introduces a new stability-related measure, the
survivability $\sigma(t)$ of a dynamical system [3]. The survivability
$\sigma(t)$ is the fraction of trajectories that stay within a desirable
regime, after an initial perturbations, up to a given time $t$. Again Monte
Carlo simulations can be used to calculate $\sigma(t)$, similar to the
calculation of the basin stability. Perturbed states $x_{0}$ are randomly
drawn and the perturbed system in state $x_{0}$ is integrated, and it is
checked whether the system leaves the desirable region or not. The
survivability is then defined as:

$\sigma(t)=\frac{V_{d}(t)}{V}$ (2)

where $V$ again is the number of all randomly drawn states and $V_{d}(t)$ is
the number of states which stay inside of the desired region for the
simulation time $t$. The single-node survivability is studied by initially
perturbing a single node $a$, starting from a stable operating state.
$\langle{\sigma_{a}(t)}\rangle=[0,1]$ is then the probability that the system
remains in the desired region, for the time $t$, after node $a$ is perturbed.
The desirable frequency region is defined as $\forall
a:|\omega_{a}|<\omega_{d}$, the region where the frequency deviations at all
nodes stay below a threshold $\omega_{d}$. The threshold $\omega_{d}$ is
chosen identical to the maximal perturbation level of the frequency
$\omega_{d}=\Delta\omega$.

For the voltage condition, we have taken inspiration from the EN 50160 Report
[11] for power quality standard[12]. The EN 50160 Report states that the
supply voltage variations must remain inside of an interval of $\pm 10\%$ for
$95\%$ of a week. We have thus decided to set the bounds for the voltage
survivability in such a way that the whole time series is evaluated and it is
checked if the voltage of a node remains outside of the $\pm 10\%$ bounds for
more than 5 seconds. In this way, we give the system enough time to mitigate
the perturbation, which might be bigger than $\pm 10\%$.

###  2.3 Feasibility

Due to the introduction of constraints into the system, there is an additional
mechanism that could result in instabilities. In DAEs finite-time solutions
exist, even when no singularities are present. The system dynamics can run off
to a state where the constraints can no longer be fulfilled. After an initial
perturbation to the power grid, it might run into a regime, where the power
flow is infeasible, which then results in a system collapse.

In our Monte-Carlo experiments, we evaluate the constraint equations $g$ at
the last time step $t_{end}$ before a system collapse occurs. If the sum of
the absolute of all constraints $g_{i}$ (3) is larger than a threshold $\nu$,
we say that the constraints have been violated:

$\sum_{i}|g_{i}(z(t_{end}))|>\nu.$ (3)

After testing different values, $\nu=0.001$ was selected as the threshold, as
it captures the true incidences where $g$ has been violated, but no cases of
purely numerical errors in $g$ were counted. We introduce the feasibility $f$
as the probability that the power flow will stay feasible after a
perturbation:

$f=\frac{V_{f}}{V}$ (4)

where $V$ is again the number of all states and $V_{f}$ is the number of
states which stay feasible, according to equation (3).

##  3 Differential-Algebraic Equations

In this paper we will describe power grids in the form of differential
algebraic equations (DAEs). For the next sections we will use the the very
general description of explicit DAEs given in equation (5) and (6):

$\displaystyle\dot{x}$ $\displaystyle=f(t,x(t),y(t))$ (5) $\displaystyle 0$
$\displaystyle=g(t,x(t),y(t))$ (6)

where equation (5) is an ODE and equation (6) are the algebraic constraint
equations of the system. We will go into details about the modeling choices
for power grids in the sections 5.1 and 5.2.

In the equations (5) and (6) $x$ depicts the variables whose derivatives
$\dot{x}$ appear in the DAE and are called _differential_. The derivatives of
the variables $y$ do not appear in the system of equations, i. e. they are
referred to as _algebraic_. In the course of the following derivation, an
explicit separation between differential and algebraic variables is not
necessary, hence we abbreviate $(x,y)$ as $z$. Since all of our equations are
autonomous, we can write $g(t,x(t),y(t))$ as $g(z(t))$.

Solving DAE systems is generally much harder than solving pure ODEs because
not every point $z$ is a solution of $g$. The problem is hence separated into
two parts: first finding valid initial conditions $z_{0}$, i. e. conditions
that fulfill the constraint equations (6), and then calculating the
trajectories [13].

Explicit DAEs can be interpreted as ordinary differential equations, whose
solutions lie on a manifold $\mathcal{M}$, which is defined by the constraint
equations $g$. Additionally the tangent space at a point $z\in\mathcal{M}$,
where $dg(z)$ denotes the total differential is defined as [14]:

$\displaystyle\mathcal{M}$ $\displaystyle:=\\{z\in\mathbb{R}^{n}:g(z)=0\\}$
(7) $\displaystyle T_{z}\mathcal{M}$ $\displaystyle=\mathcal{N}(dg(z)).$ (8)

For the scope of this paper, we limit ourselves to smooth manifolds that are
embedded in a euclidean space, meaning $\mathcal{M}\subset\mathbb{R}^{n}$. The
Euclidean space which surrounds the manifold $\mathcal{M}$ is called the
ambient space $\mathcal{S}$ of $\mathcal{M}$. In the case of power grids, this
is not a restriction as they are naturally embedded in an ambient space. Each
direction in the ambient space depicts one of the variables $z$, which will
become helpful later on.

##  4 Probability Distributions on Manifolds

While the problem of random processes on manifolds is already a long-known
issue, for example, Brownian motions on manifolds [15], the authors are not
aware that a method to draw random states from a manifold, that are associated
with a certain direction of the ambient space, does exist. To calculate
probabilistic stability measures of single nodes, these states are needed.
Thus we introduce a novel method of sampling states from manifolds, ambient
forcing.

###  4.1 Ambient Forcing

Since the manifolds in our problem are generally unknown, initial conditions
can not be directly sampled from $\mathcal{M}$. But the algebraic equations
$g(z)$ which define $\mathcal{M}$ and the so-called operation point $z_{0}$ of
the power grid, which lies on $\mathcal{M}$ since it is a fixed point, are
known. Starting from the following theorem 4.1 we will develop our approach.

###### Theorem 4.1

If $\mathcal{M}$ is a smooth manifold, every vector $v\,\in\,T_{z}\mathcal{M}$
can be interpreted as the velocity of some smooth curve $\gamma$ in
$\mathcal{M}$ [16].

Therefore if the tangent vectors $v\,\in\,T_{z}\mathcal{M}$ are known, a
relation to the curves in $\mathcal{M}$ exists. Since all manifolds which we
are interested in are embedded and smooth, the derivative can be written by
using the Jacobi matrix [16]. The constraint $g(z(t))=0$ is a chained function
$g\circ z$ and the chain rule is used to determine the derivative of
$g(z(t))$:

$\displaystyle\frac{dg(z(t))}{dt}$ $\displaystyle=J({g\circ
z})(t)=Jg(z)\cdot\dot{z}$ (9)

To conserve the constraint equations $g(z(t))$, the derivative must vanish. If
$Jg(z(t))$ is known and it is possible to identify its null space
$\mathcal{N}(Jg)$, we find the following relation:

$\displaystyle\frac{dg(z(t))}{dt}$
$\displaystyle=Jg(z(t))\cdot\dot{z}\stackrel{{\scriptstyle!}}{{=}}\vec{0}$
(10) $\displaystyle\Rightarrow\dot{z}$
$\displaystyle\in\mathcal{N}(Jg(z(t))).$ (11)

Any trajectory that satisfies equation (11) will conserve the constraints
$g(z(t_{0}))=g(z(t^{\prime}))$ for all $t^{\prime}>t_{0}$ and will thus create
states which lie on the manifold if $g(z(t_{0}))=0$ is satisfied.

Starting from equation (11), we abbreviate the null space of the Jacobian
$\mathcal{N}(z)=\mathcal{N}(Jg(z))$ and write $P^{\mathcal{N}}(z)$ for the
matrix representation of the orthogonal projection onto $\mathcal{N}(z)$. Then
for an arbitrary dynamic $\dot{z}=h(z)$, the manifold preserving version can
simply be defined as:

$\displaystyle\dot{z}=P^{\mathcal{N}}(z)h(z)$ (12)

Any differential equation that satisfies $h(z)\in\mathcal{N}(Jg(z(t)))$ has
trajectories that stay on the manifold if the initial condition $z_{0}$ lays
on the manifold.

In Single Node Basin Stability, a random displacement $F_{rand}$ of the system
in the direction of the variables of a node is considered. The directions of
the ambient space $\mathcal{S}$ are associated with variables of the nodes.
This means that $F_{rand}$ can be drawn from the directions of $\mathcal{S}$
associated with the variables. These variables are now potentially subjected
to constraints. A constraint version of such a forcing can still be defined
using the definition of equation (12):

$\displaystyle\dot{z}=P^{\mathcal{N}}(z)F_{rand}$ (13)

where $F_{rand}$ is an appropriate random, but constant, vector in the ambient
space. The definition in equation (13) is not intrinsic in terms of the
constraint manifold $\mathcal{M}$, but depends on the ambient space
$\mathcal{S}$ itself.

$F_{rand}$ is a state change corresponding to a failure in the system. The
change of state $F_{rand}$ itself is localized at a single node. A
distribution of states is obtained, associated with failures localized at
single variables or nodes without violating the constraints. Both, the amount
of $F_{rand}$ and the integration time $t$ of equation (13) determine the
final strength of the random perturbation. A compromise has to be made between
sampling strong enough perturbations for a meaningful investigation and to
preserve the network locality, which is studied in detail in section 6. To
simplify the dependence of the perturbation strengths we have decided to fix
the interval for the random integration times to $[0,1]s$ and to tune the
strength of the perturbation by adapting the intervals for the components in
$F_{rand}$. If a power grid is not constraint in the sense of equation (6),
meaning that we are in the case of pure ODEs, $g(z)$ is non-existing and we
can simply use $z=F_{rand}$, since we do not have to perform a projection.

##  5 Power Grid Modeling

So far we have used the very abstract definition of power grids as DAEs, we
want to introduce more specific models which, were used, from the network
level 5.1 down to the nodal dynamics 5.2. The nodal dynamics define the
equations of the DAEs, while the network gives the connection between them and
the power flow.

###  5.1 Network Models

The underlying structure of the power grids will be represented in the form of
graphs $\mathcal{G}=\\{N,E\\}$, where $\mathcal{G}$ is a set of nodes $N$ and
edges $E$ [17]. The edges connect the nodes with each other. In power grids,
the nodes are the buses and the edges are the transmission lines connecting
them. The nodal admittance matrix $Y$ is used in power system analysis to
calculate the power flow within the network. It is defined as follows:

$\displaystyle Y_{kl}=\begin{cases}y_{k}+\sum_{l=1,k\neq l}^{N}y_{kl}&k=l\\\
-y_{kl}&k\neq l\end{cases}$ (14)

where the admittance between two nodes $k$, $l$ is described by
$y_{kl}=g_{kl}+\mathrm{i}b_{kl}$. With $g_{kl}$ being the conductance and
$b_{kl}$ the susceptance. The term $y_{k}$ stands for the admittance of linear
loads connected to $k$ and for the admittance-to-ground at $k$ [18]. The
current injected at the nodes is given by Ohms law: $I=YU$. The nodal
Kirchhoff’s law describes the current injections at the nodes:
$\sum_{j}^{m}i_{j}=0$ where $m$ is the total number of all in- and out-coming
currents at the node. The Kirchhoff mesh rule is given by:
$\sum_{j}^{m}u_{j}=0$ where $m$ is the number of all voltages around a closed
loop. We rewrite the voltages as complex phasors:
$u_{k}=v_{k}e^{\mathrm{k}\phi_{k}}$, where $v_{k}$ is the voltage magnitude
and $\phi_{k}$ is the voltage angle. Using Kirchhoff’s and Ohm’s laws, the
power flow of the network is given by equations (16) and (17) [19]:

$\displaystyle S_{k}$ $\displaystyle=u_{k}i_{k}^{*}=P_{k}+\mathrm{i}Q_{k}$
(15) $\displaystyle P_{k}$
$\displaystyle=\sum_{l=1}^{N}v_{k}v_{j}[g_{kl}\cos(\phi_{k}-\phi_{l})+b_{kl}\sin(\phi_{k}-\phi_{l})]$
(16) $\displaystyle Q_{k}$
$\displaystyle=\sum_{l=1}^{N}v_{k}v_{l}[g_{kl}\sin(\phi_{k}-\phi_{l})+b_{kl}\cos(\phi_{k}-\phi_{l})]$
(17)

where $S_{k}$ is the apparent power at node $k$ and $P_{k}$ and $Q_{k}$ are
the real and reactive power injected at $k$ respectively.

###  5.2 Models for nodal dynamics

####  5.2.1 A normal form of grid-forming components

In this paper, we will model our dynamical nodes as grid-forming inverters.
While synchronous machines and their modeling are a well-researched topic and
many textbooks, such as [19, 18], provide a good overview, such information is
still largely missing for grid-forming inverters. While many models to
describe grid-forming inverters, such as [20, 21], exist unfortunately, a
simple, unified model, which still reflects the important properties of grid-
forming inverters, does not yet exist.

Recently [22], R. Kogler et al. have presented a normal form (18) that can
capture the most basic nonlinearities of various power system dynamics but
especially for grid forming components. The normal form is given as:

$\displaystyle v_{k}^{2}$ $\displaystyle=u_{k}u_{k}^{*}$
$\displaystyle\frac{d\omega_{k}}{dt}$
$\displaystyle=(A_{x,k}+B_{x,k}\omega_{k}+C_{x,k}v_{k}^{2}+G_{x,k}P_{k}+H_{x,k}Q_{k}))\frac{1}{M_{x,k}}$
(18) $\displaystyle\frac{du_{k}}{dt}$
$\displaystyle=(A_{u,k}+B_{u,k}\omega_{k}+C_{u,k}v_{k}^{2}+G_{u,k}P_{k}+H_{u,k}Q_{k})\cdot
u_{k}$

where $\omega_{k}$ and $u_{k}$ are the frequency and voltage at node $k$
respectively, $P_{k}$ and $Q_{k}$ represent, as in equation (15), the active
and reactive power at $k$ and $v_{k}$ is the voltage magnitude. All the other
coefficients are modeling parameters.

In [22] the authors have demonstrated that a large number of established
models, such as in [20] and [21], can be mapped to the normal form.
Furthermore, the normal form was validated by numerical simulations as well as
by experimental data. Therefore it is the perfect candidate for the nodal
dynamics of the grid-forming inverters in this work. As the derivatives of the
voltage and frequency show up in the normal form model we will also refer to
them as _differential_ nodes, following the definition in section 3.

####  5.2.2 PQ Loads

The consumers in this paper will be modeled by the traditional PQ load model.
We could have chosen other models, such as the exponential recovery load,
which reacts to the applied voltage but is nevertheless a constraint. However,
the focus of this work is on the study of the effect of constraints on power
systems, so it is reasonable to choose the simplest static model. The PQ node
locally fixes the power output of node $k$:

$\displaystyle 0=(P_{k}+\mathrm{i}Q_{k})-u_{k}\cdot i^{*}_{k}$ (19)

The voltage $u_{k}$ and current $i_{k}$ at node $k$ are now subjected to the
constraint equation. As the derivatives of the voltage do not show up in the
PQ model, we will also refer to them as _algebraic_ nodes.

###  5.3 Network Ensemble

To apply our method and explore the probabilistic stability of power grids, an
ensemble of 100 synthetic power grids is drawn by using a random growth
algorithm [23]. This algorithm generates synthetic networks which share their
topological properties with realistic power grids when a set of appropriate
parameters is chosen.

All of the networks have 100 nodes each, which are randomly turned into PQ
nodes (19) or into normal forms (18). Additionally, the nodes are randomly
turned into consumers or producers by drawing a random power output between
$[-1,1]p.u.$ for each node.

PQ nodes with a positive power output depict grid feeding inverters. When
their power output is negative, they represent traditional loads. Normal forms
with positive power output depict grid forming inverters. Those with a
negative power output depict aggregated loads, for example, a small village,
which produces solar power and is connected to the next higher grid layer via
a grid forming inverter but ultimately consumes more power than it can
produce.

In the investigation in section 6 random frequency deviation in $[-5,5]$ and
voltage perturbations in $[-2,2]$ are sampled using Ambient Forcing. Starting
from section 7.1 the frequencies of the nodes were randomly perturbed in
$[-3,2]$ and the voltages in $[-2,0]$, which depicts a probabilistic version
of a short circuit. We used $V=250$ initial states for each node in the
network.

##  6 Localized Perturbations

Ambient Forcing is intended to be network local, meaning that the initial
perturbation only targets a single node, but is not clear if it succeeds.
Unlike in ODEs, in DAE systems localized perturbations are not mathematically
well defined, since not all initial conditions are valid. We have to project
the localized perturbation vectors $F_{rand}$ from the ambient space
$\mathcal{S}$ onto the manifold.

In the case of power grids, this reflects a real-life behavior of the
constraints. When a constraint voltage $u_{a}$ is perturbed it might not be
possible for the node $a$ itself to fulfill the constraint as the current
$i_{a}$ can not be altered freely since it is additionally constraint as it
has to fulfill Kirchhoff’s. Therefore changing the voltage $u_{a}$ results in
a voltage shift at adjacent nodes, which results in an alerted current
$i_{a}$, to fulfill the desired power output.

In power grids, the perturbations will naturally spread to adjacent nodes when
constraints are present. We will characterize the network locality of the
random distribution of perturbations $\rho_{k}$ at node $k$ by introducing a
novel network measure that is related to the closeness centrality.

First the amplification $A(k)$ of an node $k$ is defined as the sum norm of
the mean displacement vector from the operation point $z_{0}$ $\Delta
z(k)=z_{0}-z_{rand}(k)$ after a perturbation on node $k$. The amplification
$A(k)$ is thus the mean total strength of a perturbation. $A(k)$ is used to
calculate the perturbation fractions $p_{a}(k)$ of all nodes $a$ in the
network:

$\displaystyle A(k)$ $\displaystyle=\sum_{l}|\Delta z_{l}(k)|$ (20)
$\displaystyle p_{l}(k)$ $\displaystyle=\frac{|\Delta z_{l}(k)|}{A(k)}$ (21)

The spreadability $s(k)$ is then defined as the sum over all fractions
$p_{l}(k)$ weighed by the shortest path length $d(k,l)$ between the perturbed
node $k$ and all the other nodes $l$. A slightly modified form of the
normalized closeness centrality $C_{c}(k)$ [1] is used to normalize $s(k)$.
The number of all nodes $n$ is used as the numerator in the modified version
$C^{\prime}_{c}(k)$ instead of the usual $n-1$:

$\displaystyle C^{\prime}_{c}(k)$ $\displaystyle=\frac{n}{\sum_{a}d(k,l)}$
(22) $\displaystyle s(k)$
$\displaystyle=C^{\prime}_{c}(k)\sum_{l}p_{l}(k)d(k,l).$ (23)

If only node $k$ is displaced the spreadability vanishes. If all nodes in the
network are equally perturbed then $s(k)=1$. The spreadability of each node in
the ensemble was calculated using equation (23). The data of the network
ensemble is categorized by the minimal constraint distance $\min(d_{c})$.
Constraint nodes themselves have a minimal constraint distance of 0. The data
is also sorted by the minimal distance to the _next constraint_
$\min(d_{cn})$, the type of the node itself is not considered. This means that
neighbors of constraint nodes will always have a distance $\min(d_{cn})$ of 1,
independently of their type.

The data is represented by a combination of a box and violin plot shown in the
figures 1(a) and 1(b). If the ambient forcing algorithm does not introduce any
further spreading of perturbations only constraints and nearest neighbors of
constraints will have a non-vanishing spreadability. This behavior can be
observed in figure 1(a), meaning that the ambient forcing algorithm can be
used for single node perturbations. Furthermore, in figure 1(b), it can be
seen that the spreadability decreases as the distance to the next constraint
in the power grid increases, independently of the node type.

(a)

(b)

Figure 1: Box plots of the spreadabilities using different sorting mechanisms
for the distance. Figure 1(a) shows that only constraint nodes or neighbors of
constraint nodes spread. The second plot 1(b) shows that the spreadability of
the nodes decreases with the distance to the next constraint in the system,
independently if the node itself is a constraint or not.

##  7 Results

###  7.1 Topological Properties

The algorithm, which was introduced in [4], classifies nodes according to
their topological properties and can identify vulnerable nodes in power grids.
In this algorithm nodes with degree $d=1$ are iteratively removed, introducing
a sub-graph at each level.

We used the topological classiﬁcation scheme to investigate whether it is
still possible to separate the nodes according to their stability properties
when constraints, as well as the normal form, are introduced into the network.
While we can report that the classes are still able to separate the nodes with
respect to their stability properties, we have found that the degree is a much
better predictor than the topological classes [4], for both $\sigma$ and
$\beta$ but especially for the survivability.

To illustrate the connection between the degree and the probabilistic
measures, we have plotted the mean $\beta$ and $\sigma$ with respect to the
degree in the figures 2(a) and 2(b). The solid lines represent the mean and
the shaded area depicts the $15.9\%$ and $84.1\%$ percentile, which depict one
standard deviation in the normal distribution.

The basin stability, in figure 2(a), shows a weak upwards trend regarding the
degree. The picture for the survivability on the other hand clearly exhibits
that for both the algebraic and differential nodes the mean survivability
increases with the degree. This is in stark contrast to the finding of the
original paper [4], where the survivability decreases with the degree. For the
algebraic nodes, we see that the size of the percentiles decreases with the
degree, while it increases with the degree for the differential nodes.

We have verified that the strength of the mean perturbations $A$ does not
decrease with the degree. The results of this investigation can be found in
the appendix Amplification.

(a) Basin Stability

(b) Survivability

Figure 2: Degree dependence of the basin stability 2(a) and survivability
2(b). Solid lines represent the mean and the shaded areas depict the $15.9\%$
and $84.1\%$ percentile. Both the basin stability and survivability increase
with the degree.

As the effect of the degree on the survivability is especially prominent, we
have investigated this behavior further by calculating the probability
densities of the survivability for all nodes in the ensemble. For the
algebraic nodes, whose results are given in Fig. 3(a), we find that nodes with
a degree of 8 or higher survive most often and their mode is $\sigma=1$. Nodes
with smaller degrees survive less frequently and the distributions start to
broaden.

The survivability of the differential nodes, shown in figure 3(b), exhibits
the most striking dependence on the degree. Increasing the degree in the nodes
results in a higher mode of the survivability. The distribution of $\sigma$
becomes broader and more long-tail with an increasing degree. This is in
contrast to the behavior of the algebraic nodes where the distributions become
sharper with an increasing degree.

Additionally, we want to draw attention to the fine structuring, which becomes
visible in the survivability of the differential nodes. Already starting at a
degree of $d=2$, a second peak with higher survivability emerges. This hints
that there might be a second effect that influences the survivability.

(a) Algebraic

(b) Differential

Figure 3: Density Plot of the survivability classified by the degree of the
nodes. The left figure shows all algebraic nodes, while the differential nodes
are shown on the right. It can be seen that the survivability increases with
the degree, especially for the differential nodes.

###  7.2 Types of Instabilities in the System

Different possible mechanisms could result in instabilities and it is
worthwhile investigating which of these mechanisms are most important. The
following classification for instabilities will be used. First, we count
events where the frequency of at least one node does not return to the
synchronous state, i. e. that the system has desynchronized. Additionally we
compute the states where the voltage of at least one node does not return to
the fixed point. We will also fine-grain these and measure the events where
_only_ the frequency is desynchronized and then measure the states where only
the voltage runs off to another fixed point. Finally, we will look at all the
nodes where infeasible power flows have occurred.

We investigate the percentage of nodes $p_{n}$ in the ensemble, where these
events have occurred at least once. The results are shown in the table 1. We
can see that the desynchronization only occurs at 19.1$\%$ of the nodes at
all, and never occur on their own, without a change in voltage. Two example
trajectories of these desynchronized states can be found in the appendix 10.
Voltage drops occur at 29.5$\%$ of the nodes in the network and also appear
without a desynchronization, at 11.4$\%$ percent of the nodes. Two examples of
voltage drops that occur without desynchronization can be found in the figures
9. Finally, infeasible power flows appear most commonly, at 37.1$\%$ of the
nodes.

Faults | Desynchronization | Voltage Drops | Only Voltage Drops | Infeasible Flows
---|---|---|---|---
$p_{n}$ | 19.1 $\%$ | 29.5 $\%$ | 11.4 $\%$ | 37.1 $\%$
Table 1: Percentage of nodes $p_{n}$ in the ensemble where faults have
appeared at least once.

As the infeasible power flows occur most commonly we have decided to compare
the feasibility $f$ to the basin stability and survivability. For this
purpose, we have calculated the probability densities of $\beta$, $\sigma$,
and $f$ for the entire ensemble and compared the distributions using the
Jensen-Shannon divergence[24], which is a symmetric and bounded version of the
Kullback-Leibler divergence $D_{KL}(P,Q)$ [25]. $D_{KL}(P,Q)$ is a statistical
distance measure between two probability distributions $P$ and $Q$, which are
defined in the same probability space $X$:

$\displaystyle D_{KL}(P,Q)=\sum_{x\in
X}P(x)\log\left(\frac{P(x)}{Q(x)}\right)$ (24)

and it is a measure for the information which is lost if we approximate $P(X)$
using $Q(X)$. The Jensen-Shannon-divergence for $P(X)$ and $Q(X)$ is then
defined as:

$\displaystyle JSD(P,Q)$
$\displaystyle=\frac{1}{2}D_{KL}(P,M)+\frac{1}{2}D_{KL}(Q,M)$ (25)
$\displaystyle M$ $\displaystyle=\frac{1}{2}(P(X)+Q(X)).$ (26)

The Jensen-Shannon-divergence is bounded on the interval $[0,1]$ and
$JSD(P,Q)$ is 0 if $P(X)$ and $Q(X)$ are equal and 1 if both probability
distributions are not correlated.

The Jensen-Shannon-divergence can be calculated for $\beta$, $\sigma$, and $f$
as they are all defined in the same probability space and the possible
outcomes $x\in X$ are defined on the interval of $[0,1]$. The Jensen-Shannon
divergence $JSD(\beta,f)$ for the algebraic nodes is 0.079 and for the
differential nodes, it is 0.013. $JSD(\sigma,f)$ is 0.15 for the algebraic
nodes and 0.248 for the differential nodes. This means that little information
is lost if we use $f$ to approximate $\beta$, while a lot of information is
lost if $f$ is used to approximate $\sigma$. From this, we can conclude that
the infeasible power flows are the driving mechanism that lead to basin-
unstable systems, while other effects dominate the survivability.

Finally, we can also infer that there are no differences between producers or
consumers, i.e. the sign of the power, concerning the basin stability or
survivability. This means that not only inverters are at risk of inducing
voltage drops or infeasible power flows but also loads. Therefore it is
necessary to pay close attention to all types of actors in the grid.

From this, we can conclude that although algebraic nodes themselves are more
stable and survive more frequently than differential nodes, which can be seen
in the figures 2(a) and 2(b), they are responsible for most of the
instabilities in the grid, the infeasible power flows. This is another
argument for pushing ahead with the expansion of grid-forming inverters in
order not to increase the number of constraints in the grid even further,
since loads, which introduce the same problems, will always be connected to
power grids.

A more detailed load model is needed to properly understand this phenomena but
a fast energy storage at the loads could be used to change the power flows in
time and mitigate the infeasabilities. In a future investigation _Monte Carlo
basin bifurcation analysis_ [26] could be used to investigate the different
asymptotic states in such systems and their basin of attraction.

##  8 IEEE 96 System

Finally, we will also apply our new method to a more realistic test case. The
IEEE 96 test system [5] was designed for bulk power system reliability
evaluation studies. This test system is not representative of any specific or
typical power system, as it should hold universal characteristics to be useful
as a reference for testing the impact of different evaluation techniques. The
system should thus be seen as a hybrid and atypical system. To update the
system and show the influence of inverters, we have decided to model the
generators as normal forms (18). The loads are still depicted by traditional
PQ nodes (19). The parameters used for the normal form can be found in the
appendix Phase Amplitude Oscillators Parameters.

###  8.1 EN 50160 Report Perturbations

The perturbation parameters were chosen following the EN 50160 Report [11] on
Power Quality Standards. In this standard, the grid frequency $\omega$ must
remain inside of the bounds $\Delta\omega=[-3.0,2.0]$ for 100$\%$ of the
observed time. The random frequency perturbations are thus sampled from this
interval and the thresholds for the survivability are chosen identical. We
have decided to first perturb the voltage of each node in the network between
$\pm 10\%$ from its operation point level. For the voltage survivability, it
is checked whether the voltage of a node remains outside of $\pm 10\%$ of the
operation point level for more than 5 seconds.

These results show that the minimal basin stability of a node is 0.983$\pm$
0.006 and the minimal survivability is 0.95 $\pm$ 0.01. Thus it was decided to
not include a more extensive evaluation of the data as there is little
additional information that can be gathered, since all nodes are almost
globally stable. Nevertheless, it is an interesting result as it demonstrates
that the IEEE test case with inverters fulfills the requirements of the EN
report.

###  8.2 Increased Perturbation Strength

To gain further information, we increase the strength of the frequency
perturbations to $\Delta\omega=[-5,5]$, and the strength of the voltage
perturbations to $\pm 25\%$ of the operation point level. The results of the
basin stability and survivability are shown in the figures 4(a)-7(b). The
histograms in the figures 4(a) and 4(b) show that all PQ Node are globally
stable meaning that they always return to the initial fixed point after a
perturbation and never leave the desired region. For the normal form, some
nodes also have basin stability and survivability of 1. But some of the
differential nodes also become less stable. Additionally, it was found that
for all nodes the survivability $\sigma$ is smaller or equal to the basin
stability $\beta$, which can be seen in Figure 6 in the appendix.

(a) Basin Stability

(b) Survivability

Figure 4: Stacked Histograms of the Basin Stability 4(a) and survivability
4(b), grouped by the function type. The basin stability of differential and
algebraic nodes are similarly distributed. The second figure shows a clear
distinction in the distribution of the survivability between differential and
algebraic nodes.

It is interesting to check which of the EN 50160 Report requirements, for the
voltage amplitude $v$ or the frequency $\omega$, were broken. It should be
noted that the different $\beta$ and $\sigma$ results form the same set of
initial perturbations, just the two different conditions were evaluated.

We found that the basin stability remains the same when we evaluate the
voltage or the frequency condition. This means that when the frequency does
not return to the initial fixed point, the voltage will also not and vice
versa. For the survivability, we found that only evaluating the frequency
condition will overestimate the survivability. This indicates that there are
many nodes whose voltages leave the desired region while the frequency does
not. This evaluation is shown in the figures 7(a) and 7(b) in the appendix.

###  8.3 Voltage Sweep

As we have not yet seen an effect of voltage perturbations, we increase the
maximal perturbations up to $\pm 500\%$ of the operation point level. Although
only the voltage was perturbed, we still evaluate the frequency bounds for the
survivability. The bounds are again chosen following the EN report, meaning
that the threshold is $[-3,2]$. The results of the voltage sweep are shown in
the figures 5(a) to 5(b).

Up to a perturbation level of $\pm 100\%$, we do not see an effect on basin
stability and survivability and all nodes are globally stable. Increasing the
perturbation further, we see that the basin stability and survivability
decrease. We discover that the dynamical nodes have a lower $\beta$ and
$\sigma$ than the algebraic nodes. Just as before we see that the
survivability $\sigma$ is smaller than the basin stability $\beta$.

(a) Basin Stability

(b) Survivability

Figure 5: Dependence of the basin stability and survivability on the strength
of the voltage perturbations. The basin stability first starts to react to the
perturbation with a level of $\pm 100\%$ from the operation point. The basin
stability is generally higher than the survivability. Differential nodes are
less stable and survive less frequently than algebraic nodes.

Although the voltage is crucial for the survivability of the IEEE96 test case,
it is not the realistic voltage perturbations which destabilize them, as the
voltage perturbations only start to play a role at perturbation strengths of
around $\pm 100\%$. This implies that the frequency perturbations, in the
regime of the realistic perturbations, pushes the voltages out of the desired
regime. This further illustrates that it is important to associate several
stability classes with each other.

##  9 Conclusion

In this paper, we have introduced a novel method for probability distributions
on manifolds that can generate random, initial conditions for probabilistic
stability methods like basin stability and survivability. By introducing the
new network measure, the spreadability, we have proven that the initial
conditions are network local in a DAE sense and that the states are thus
usable for single node basin stability and survivability.

We have utilized the method to calculate the basin stability and survivability
of an ensemble of random synthetic power grids, which includes grid following
as well as grid forming inverters. This was the first time basin stability and
survivability were calculated for systems with constraints. We have applied
the topological classification scheme, which was introduced by Nitzbon et. al.
[4], and investigate if it still holds. We can report that, while the classes
are still able to separate the nodes with respect to their stability
properties, the degree is a much better predictor for the basin stability and
especially the survivability. While the degree is a good predictor we want to
draw attention to the fine structure inside of the distributions of the
survivability which indicates that there might be another effect that affects
it.

Interestingly we were able to show the inverse effect of the degree on the
survivability than the original paper [4]. We found that the survivability of
the nodes increases with the degree, whereas before $\sigma$ decreased with a
higher degree.

Additionally, we do not see the frequent desynchronization events which are
normally observed when only the swing equation is used as the nodal model. An
important destabilizing mechanism is that the systems frequently run into
infeasible power flow solutions, meaning that the constraints can no longer be
fulfilled. This could potentially be avoided by adding fast energy storage to
the loads and grid following inverters to reach power balance. A more detailed
model of loads is needed for such an investigation.

We also found out that algebraic nodes are more stable and survive more often
after a perturbation than the differential nodes. But they still play a
crucial role in the destabilization of the system as they constraint the power
flow and thus create infeasible solutions. Generally, the basin stability is
higher than the survivability.

Finally, we have applied the method to calculate the basin stability of the
IEEE 96 system. For smaller, more realistic perturbations we see that again
the algebraic nodes are more stable than the differential nodes. In fact that
all algebraic nodes have basin stability and survivability of 1. Thus we have
increased the strength of the voltage perturbation and found that voltage
perturbations only start to play a role if we perturb the system from its
fixed point $\pm 100\%$, meaning that realistic voltage perturbations do not
play a role in the destabilization of the system.

## Code availability

All code to reproduce the results and figures of this paper is available at
the GitHub repository https://github.com/PIK-
ICoNe/AmbientForcing_Paper_Companion or on zenodo
https://doi.org/10.5281/zenodo.5901260. The implementation of the Ambient
Forcing Method can be found in a separate GitHub repository
https://github.com/PIK-ICoNe/AmbientForcing.jl.

## Acknowledgments

A. Büttner acknowledges support by the German Academic Scholarship Foundation.
All authors gratefully acknowledge the European Regional Development Fund
(ERDF), the German Federal Ministry of Education and Research and the Land
Brandenburg for supporting this project by providing resources on the high-
performance computer system at the Potsdam Institute for Climate Impact
Research.

## References

  * [1] L.. Freeman  “Centrality in social networks conceptual clarification”  In _Social Networks_ 1, 1978, pp. 215–239
  * [2] P. Menck, J. Heitzig, N. Marwan and J. Kurths  “How basin stability complements the linear-stability paradigm”  In _Nature Phys_ 9, 2013
  * [3] F. Hellmann, P. Schultz, C. Grabow, J. Heitzig and J. Kurths  “Survivability of Deterministic Dynamical Systems”  In _Nature Scientific Reports_ , 2016  DOI: DOI: 10.1038/srep29654
  * [4] J. Nitzbon, P. Schultz, J. Heitzig, J. Kurths and F. Hellmann  “Deciphering the imprint of topology on nonlinear dynamical network stability”  In _New J. Phys. 19 033029_ , 2017
  * [5] Reliability Test System Task Force  “The IEEE Reliability Test System 1996”  In _IEEE Transactions on Power Systems_ 14, 1999
  * [6] P.. Menck, J. Heitzig and J.. J.  “How dead ends undermine power grid stability”  In _Nature Commun_ 5, 2014
  * [7] N. Hatziargyriou, J. Milanovic, C. Rahmann, V. Ajjarapu, C. Canizares, I. Erlich, D. Hill, I. Hiskens, I. Kamwa, B. Pal, P. Pourbeik, J. Sanchez-Gasca, A. Stankovic, T. Cutsem, V. Vittal and C. Vournas  “Definition and Classification of Power System Stability – Revisited and Extended”  In _IEEE Transactions on Power Systems_ , 2021  DOI: 10.1109/TPWRS.2020.3041774
  * [8] J. Hossain and H.R. Pota  “Robust Control for Grid Voltage Stability: High Penetration of Renewable Energy”  Springer, 2014
  * [9] R. Feynman  “Feynman Lectures on Physics Volume 2”  Addison-Wesley, 1964
  * [10] P.N. Brown  “Consistent Initial Condition Calculation for Differential-Algebraic Systems”  In _Society for Industrial and Applied Mathematics_ , 1998
  * [11] “EN 50160, Voltage characteristics of electricity supplied by public distribution systems”, 1999
  * [12] “Voltage characteristics of electricity supplied by public electricity networks; German version, DIN EN 50160:2020-11”, 2020  DOI: https://dx.doi.org/10.31030/3187943
  * [13] E. Hairer  “Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations”  Springer, 2004
  * [14] K. Königsberger  “Analysis 2”  Springer, 2004
  * [15] N.. Kampen  “Brownian motion on a manifold”  In _Journal of Statistical Physics_ , 1986
  * [16] John M. Lee  “Introduction to Smooth Manifolds”  Springer, 2012
  * [17] M. Newman  “Networks. An Introduction”  Oxford University Press, 2010
  * [18] P.. Sauer and M.. Pai  “Power System Dynamics and Stability”  Springer, 2007
  * [19] J. J., J.W. Bialek and J.R. Bumby  “Power System Dynamics - Stability and Control”  Wiley $\&$ Sons, 2008
  * [20] G. Seo, M. Colombino, I. Subotic, B. Johnson, D. Gross and F. Dörfler  “Dispatchable Virtual Oscillator Control for Decentralized Inverter-dominated Power Systems: Analysis and Experiments”  In _IEEE Applied Power Electronics Conference and Exposition_ , 2019, pp. 561–566
  * [21] J. Schiffer, R. Ortega, A. Astolfi, J. Raisch and T. Sezi  “Conditions for stability of droop-controlled inverter-based microgrids”  In _AUTOMATICA_ 50, 2014  DOI: 10.1016/j.automatica.2014.08.009
  * [22] R. Kogler, A. Plietzsch, P. Schultz and F. Hellmann  “A normal form for grid forming power grid components”, 2021  arXiv:2106.00644 [nlin.AO]
  * [23] P. Schultz, J. Heitzig and J. Kurths  “A random growth model for power grids and other spatially embedded infrastructure networks”  In _Eur. Phys. J. Spec. Top._ 223, 2014, pp. 2593–2610  DOI: 10.1140/epjst/e2014-02279-6
  * [24] B. Fuglede and F. Topsoe  “Jensen-Shannon divergence and Hilbert space embedding”  In _International Symposium on Information Theory, 2004. ISIT 2004. Proceedings._ , 2004  DOI: 10.1109/ISIT.2004.1365067
  * [25] S. Kullback and R.. Leibler  “On Information and Sufficiency”  In _The Annals of Mathematical Statistics_ , 1951  DOI: 10.1214/aoms/1177729694
  * [26] M. Gelbrecht, J. Kurths and F. Hellmann  “Monte Carlo basin bifurcation analysis”  In _New Journal of Physics_ 22 IOP Publishing, 2020  DOI: 10.1088/1367-2630/ab7a05
  * [27] A. Agresti and B. Caffo  “Simple and Effective Confidence Intervals for Proportions and Differences of Proportions Result from Adding Two Successes and Two Failures”  In _The American Statistician_ American Statistical Association, Taylor & Francis, Ltd., 2000  DOI: https://doi.org/10.2307/2685779

## Appendix

### Phase Amplitude Oscillators Parameters

| $A$ | $B$ | $C$ | $G$ | $H$
---|---|---|---|---|---
$u$ | 0.075 $\cdot Q$ | 1$i$ | -0.0625 | 0 | -0.0125
$\omega$ | 5 $\cdot P$ | -1 | 0 | -5 | 0
Table 2: Normal form coefficients in terms of the original parameters. $A$
depends on the active and reactive power sets points, $P$ and $Q$
respectively, at the node node.

### IEEE96

Figure 6: Comparison between the basin stability $\beta$ and the survivability
$\sigma$ of all nodes in the IEEE 96 test case. It can be seen that the
survivability is generally smaller or equal to the basin stability. The
confidence intervals were calculate using ”add two successes and failures”
[27].

(a) Basin Stability

(b) Survivability

Figure 7: Comparison between the basin stability and survivability of the IEEE
96 test case evaluated for the different conditions. Either only the voltage
condition, denoted by an $u$ index, is evaluated or only frequency condition,
denoted by the $\omega$ suffix. The right figure shows that the basin
stability does not change when we evaluate different conditions, meaning that
when the frequency does not return to the initial fixed point the voltage will
also not and vice versa. The left figure shows that for the survivability
evaluating only the frequency condition is not enough as there are many nodes
in which voltage leaves the desired region while the frequency does not. The
confidence intervals were calculate using ”add two successes and failures”
[27].

### Amplification

(a) Density Estimation

(b) Degree Dependence

Figure 8: Distribution of the amplification $A$. The mode of the distributions
lies at 1.1. The second smaller peak at $A=0$ occurs because of the slack
nodes whose voltage can not be perturbed. In the right figure, we can see that
the mean amplification slightly increases with the degree. For higher degrees,
we can see that the size of percentiles increases. This dependency on the
degree can not explain the strong increase of the survivability with the
degree. If this would have been the case we would need to see a decrease of
$A$ with the degree.

### Power Grids and Example Trajectories

In this section, we will show some exemplary trajectories which occur after
perturbations on single nodes, as well as the corresponding synthetic power
grids. The nodes in the power grid are classified using the algorithm from
[4]. Algebraic nodes are depicted as triangles while differential nodes are
shown as circles. The size of the nodes corresponds to the single node basin
stability.

(a) Node 94

(b) Node 13

(c) Power Grid

Figure 9: Example of perturbations on nodes that lead to voltage drops without
a frequency desynchronization. Node 77, which is also labeled, runs to the
same fixed point as node 94 after a perturbation.

(a) Perturbation on a Sparse Sprout

(b) Power Grid

Figure 10: Desynchronization and a short circuit after a perturbation on the
sparse sprout 62.

(a) Perturbation on an Inner Tree Node

(b) Power grid

Figure 11: Perturbation on an inner tree node that leads to a
desynchronization and voltage drops.