# Locating line and node disturbances in networks of diffusively coupled
dynamical agents

Robin Delabays1,2, Laurent Pagnier3, and Melvyn Tyloo4,5 1 Automatic Control
Laboratory, ETH Zurich, CH-8092 Zurich, Switzerland.
2 Center for Control, Dynamical Systems and Computation, UCSB, Santa Barbara,
CA 93106-5070 USA.
3 Department of Mathematics, University of Arizona, Tucson, AZ 85721, United
States.
4 Department of Quantum Matter Physics, University of Geneva, CH-1211 Geneva,
Switzerland.
5 School of Engineering, University of Applied Sciences of Western Switzerland
HES-SO, CH-1951 Sion, Switzerland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A wide variety of natural and human-made systems consist of a large set of
dynamical units coupled into a complex structure. Breakdown of such systems
can have a dramatic impact, as in the case of neurons in the brain or lines in
an electric grid, to name but a few. Preventing such catastrophic events
requires in particular to be able to detect and locate the source of
disturbances as fast as possible. We propose a simple method to identify and
locate disturbances in networks of coupled dynamical agents, relying only on
time series measurements and on the knowledge of the (possibly Kron-reduced)
network structure. The strength and the appeal of the present approach lies in
its simplicity paired with the ability to precisely locate disturbances and
even to differentiate between line and node disturbances. If we have access to
measurement at only a subset of nodes, our method is still able to identify
the location of the disturbance if the disturbed nodes are measured. If not,
we manage to identify the region of the network where the disturbance occurs.

Keywords: Line disturbance, node disturbance, detection and localization.

##  1 Introduction

Despite recent advances in network science [1, 2], the understanding of large,
networked dynamical system is still incomplete, even though such systems can
play a major role in our daily life. Large scale events therein should be
monitored or mitigated, if not completely avoided. Preventing catastrophic
events in a networked system requires, in particular, an accurate and reliable
assessment of its state in real time. Indeed, the early detection of a source
of disturbances and its precise location allow preventive measures to be
taken.

In networks, the most fundamental disturbances can affect either of their two
basic constituents, namely either nodes or lines. Nodal perturbations,
occurring mostly as variations of intrinsic velocities of the disturbed
agents, act as an additive perturbations, whereas disturbances on lines, which
are modifications of the coupling parameters, are represented as
multiplicative disturbances in the model, which are much harder to tackle
analytically. This partly explains the scarcity of the literature covering
line disturbances. Nevertheless, such perturbations are at least as important
as nodal perturbations and, in our opinion, should be investigated in details.

Systems that can be modeled as a set of diffusively coupled dynamical agents
where state estimation is needed for security, comfort, or diagnostic reasons
are extremely diverse. To name but a few, they range from neurons in the brain
to transportation networks and disease spreading models [3, 4, 5]. The case of
high voltage electrical grids is an important example of such systems and has
attracted a lot of work along the last years, especially in the context of the
ongoing energy transition in many parts of the world [6, 7]. We will
systematically refer to the power grid and often use its terminology in this
manuscript for illustrative purposes, even though our considerations extend to
a much larger class of systems.

In the context of power grids, recent applications of results from network
science and dynamical systems theory allowed the development of new techniques
aiming at a fast assessment of the system’s state. Such improvements are based
on the measurement of some quantities such as voltage amplitudes, phases, and
frequencies, which are nowadays widely accessible on a high time resolution
thank, for instance, to Phasor Measurement Units (PMUs). Part of the work on
state estimation, such as Refs. [8, 9, 10], is mostly focused on the
estimation of the current operating state of the whole system. This is
performed using various techniques, such as estimation of the eigenmodes of
the network through probing signals [8] or resonance methods [10], or based on
the measurement of the response of the system to ambient noise [9]. Closer to
our interest in this manuscript, Refs. [11, 12, 13, 14, 15] proposed some
approaches aiming at locating the source of a disturbance. For networks
composed of areas with weak inter-area connections and strong intra-area
connections, Ref. [11] uses the residues of an estimated transfer function in
order to locate a nodal disturbance. Other approaches rely on _Discrete
Wavelet Transform (DWT)_ [12, 13] or logistic regression [14] of PMU
measurements in order to identify the source of a disturbance. Most of these
methods are designed to locate a nodal disturbance only.

Despite their relevance in electrical networks [16, 17, 18], the case of line
disturbances led to fewer results compared to nodal disturbances [19, 20], and
less focus was put on assessing the impact of a line perturbation. Methods
based on the disturbance propagation [15] would be able to determine the area
where a disturbed line is located, but to the best of our knowledge, there are
no such methods specifically designed to locate line disturbances.

Recent works analyzed the case of line failures in power grids, e.g.,
following an adversary attack [21, 22, 23, 24, 25, 26, 27]. These approaches
assumes full knowledge of the system prior to the line failure, and analyzes
the mismatch between the expected and actual voltages and powers to recover
the set of failed lines, by optimizing over the potential fault locations. On
top of the line failures, the authors of [22, 23, 24, 25] assume that the
attacker is able to block data acquisition in the area of the fault, and
sometimes even sends some spurious data to the monitoring device. In the
approach of [26, 27], the authors assume only partial access to measurements,
and find the most likely location for the failure. However, in these
approaches, the authors assume that the system operator is aware of the
attack, and their analysis follows the line failure, which means that dramatic
events already happened.

In this manuscript, we propose to detect and locate disturbances whose impact
on the system is not (yet) threatening its operation, without any a priori
knowledge of the faults’ characteristics or approximate location. Our method
relies on the measurements of the state of the dynamical agents composing the
system and on the knowledge of the interactions between them when the system
is unperturbed. This means that, in theory, our method can be applied online,
without requiring to know if the system is perturbed or not. Conveniently, we
cover the case of disturbed nodes as well as disturbed lines in a unified
framework. Assuming that the amplitude and the rate of change of these
disturbances are not too large, we propose a way to locate the faulty elements
(nodes and/or lines) in the network, based on measurement of the nodes’
trajectories. In a similar spirit as Refs. [22, 23, 24, 25], but in a more
general framework, we locate the source of the disturbances by inspection of
the mismatch between expected and actual agents’ velocities. If we have access
to measurements at the faulty elements, we are able to precisely locate it,
even though such slow disturbance will spread throughout the whole network and
have a similar impact on each agent. In case of partial measurements, we are
nevertheless able to determine the area of the network where the faulty
element is located. Furthermore, our implementation is adapted to a continuous
monitoring of the system and can be used to preventively detect and locate
deficient elements before the occurrence of a more significant event, such as
a cascading failure, and that at a rather low cost as our method does not rely
on any optimization but only on a single matrix-vector product at each time
step.

The manuscript is organized as follows. We recall some preliminary tools and
give a description of the type of models considered in Section 2. In Sections
3 and 4, we detail the method to locate the faulty element for a single and
multiple disturbances, and numerical illustrations are given in Section 5.

##  2 Notations and framework

Throughout this manuscript, we denote by $\bm{e}_{i}$ the $i$th vector of the
canonical basis, with $1$ at index $i$ and zero everywhere else, and
$\bm{e}_{ij}=\bm{e}_{i}-\bm{e}_{j}$. We write ${\rm
diag}(\\{m_{i}\\})\in\mathbb{R}^{n\times n}$ for the diagonal matrix with
elements $m_{1},m_{2},...,m_{n}$ on its diagonal and $\bm{1}_{n}$ (resp.
$\bm{0}_{n}$) the vector of ones (resp. zeros) of length $n$, where the
subscript is omited when the dimension is clear from the context. Finally, we
denote by $\mathbb{I}$ the identity matrix.

###  2.1 Networked diffusive systems

We model the dynamics of a set of $n\in\mathbb{N}$ diffusively interacting
agents as

$\displaystyle
m_{i}\ddot{x}_{i}+d_{i}\dot{x}_{i}=\omega_{i}-\sum_{j}a_{ij}f(x_{i}-x_{j})\,,\qquad
i\in\\{1,...,n\\}\,,$ (1)

where the real numbers $m_{i}$, $d_{i}$, and $\omega_{i}$ are respectively the
inertia, the damping, and the natural velocity of agent $i$,
$a_{ij}\in\mathbb{R}$ is the $(i,j)$th element of the weighted adjacency
matrix of the interaction graph, and $f\colon\mathbb{R}\to\mathbb{R}$ is an
odd, differentiable coupling function, such that $f^{\prime}(0)>0$. We assume
the interaction graph to be undirected, i.e., $a_{ij}=a_{ji}$ (the directed
case is discussed in the Appendix). Note that we will be interested in fixed
points of Eq. (1), i.e., when the left hand side vanishes, meaning that all
our results apply to models with higher order time derivatives.

Assume Eq. (1) has a fixed point, denoted $\bm{x}^{*}\in\mathbb{R}^{n}$. When
the system is not too heterogeneous, i.e., the spread of $\\{\omega_{i}\\}$ is
small relatively to the magnitude of the couplings, the fixed point components
are small (in absolute value) and are well-approximated by the linearized
equation

$\displaystyle\bm{\omega}=\bm{L}\bm{x}\,,$ (2)

where $\bm{L}\in\mathbb{R}^{n\times n}$ is the Jacobian matrix defined as,

$\displaystyle L_{ij}=-\frac{\partial F_{i}(\bm{0})}{\partial x_{j}}\,,\qquad
i,j\in\\{1,...,n\\}\,,$ (3)

where $F_{i}(\bm{x})$ is given by the sum in the right-hand side of Eq. (1).
This Jacobian matrix satisfies the properties (i) $\sum_{j}L_{ij}=0$ and (ii)
$L_{ii}=-\sum_{j}L_{ij}$ $\forall i$, it is then a (weighted) Laplacian matrix
corresponding to the interaction graph of the system Eq. (1), with weights
related to the strength of the couplings between agents. This approximation is
common in the context of power grids, where it is called _DC power flow
equations_ [28, Sec. 9.7]. Obviously, for linear coupling functions, Eq. (2)
is not an approximation.

Remark. Throughout this manuscript, we will refer to dynamical agents as
_nodes_ and to the connections between them as _lines_. In the context of
graph theory, these two objects are usually called _vertices_ and _edges_
respectively.

###  2.2 Disturbances

We consider disturbances either at nodes, yielding the time varying natural
velocity

$\displaystyle\omega_{i}(t)=\omega_{i}^{0}+\xi_{\rm n}(t)\,,$ (4)

or on lines, described as varying coupling

$\displaystyle a_{ij}(t)=a_{ij}^{0}+\xi_{\rm l}(t)\,,$ (5)

where we take $\xi_{\rm n,l}(t)$ to be any function of time, whose rate of
change is typically smaller than any intrinsic time scale of the system Eq.
(1). Formally, we require that the disturbance’s rate of change is smaller
than the rate at which its effect is damped at the agents and smaller than the
speed at which it spreads throughout the network, i.e., for all $i,j$,

$\displaystyle\max_{t}|\dot{\xi}_{\rm n,l}(t)|$
$\displaystyle\ll\min\left\\{\frac{d_{i}}{m_{i}},\frac{\lambda_{j}}{\sqrt{m_{i}}},\frac{\lambda_{j}}{d_{i}}\right\\}\,,$
(6)

where $0=\lambda_{1}<\lambda_{2}<...<\lambda_{n}$ are the (real) eigenvalues
of the Jacobian matrix $\bm{L}$. This first requirement is reasonable in the
sense that a disturbance that varies too fast will not spread in the network,
and thus either its source is rather easy to locate or it has no influence on
the observed nodes.

In order to guarantee that the system remains in the vicinity of the initial
fixed point, we also require that the amplitude of $\xi_{\rm n,l}$ is not too
large,

$\displaystyle\max_{t}|\xi_{\rm
n}(t)|\ll\omega_{i}^{0}\,,\qquad\max_{t}|\xi_{\rm l}(t)|\ll a_{ij}^{0}\,.$ (7)

These last assumptions allow us to consider that the system remains in the
linear approximation regime. Again, this assumption is reasonable because a
disturbance with too large amplitude would compromise the normal operation of
the system, and its location then becomes a secondary problem. We will assume
Eqs. (6) and (7) to be satisfied throughout the manuscript.

###  2.3 The Sherman-Morrison-Woodbury formula

Our results rely heavily on the _Sherman-Morrison-Woodbury formula_ [29, Sec.
2.1.4], giving an explicit formulation for the inverse of the rank-$k$
perturbation of an invertible matrix $\bm{A}$,

$\displaystyle\left(\bm{A}+\bm{U}\bm{V}\right)^{-1}=\bm{A}^{-1}-\bm{A}^{-1}\bm{U}(\bm{I}+\bm{V}\bm{A}^{-1}\bm{U})^{-1}\bm{V}\bm{A}^{-1}\,,$
(8)

where $\bm{U}\in\mathbb{R}^{n\times k}$ and $\bm{V}\in\mathbb{R}^{k\times n}$
characterize the rank-$k$ perturbation. We emphasize that the Sherman-
Morrison-Woodbury formula applies to the pseudoinverse of Laplacian matrices
(even though such matrices are singular), provided that the columns (resp.
rows) of $\bm{U}$ (resp. $\bm{V}$) are orthogonal to the kernel of $\bm{A}$.

In the case where $k=1$, Eq. (8) reduces to

$\displaystyle\left(\bm{A}+\bm{u}\bm{v}^{\top}\right)^{-1}=\bm{A}^{-1}-\frac{\bm{A}^{-1}\bm{u}\bm{v}^{\top}\bm{A}^{-1}}{1+\bm{v}^{\top}\bm{A}^{-1}\bm{u}}\,,$
(9)

where $\bm{u}$ and $\bm{v}$ are vectors characterizing the rank-1
perturbation. Eq. (9) is usually referred to as the _Sherman-Morrison
formula_.

###  2.4 The Kron reduction

In the context of electrical networks, it is possible to rewrite the equation
governing the currents (the power flow equations [30]) with respect to a
restricted number of voltage variables through _Kron reduction_ [31]. This is
done by taking the Shur complement [29, Sec. 3.2.11] of the coupling matrix
with respect to a subset of nodes and adapting the velocities accordingly.
Indeed, this reduction is not restricted to power grids and can be applied to
the coupling matrix of any network. Partitioning the nodes in two sets
$I_{g}=\\{1,...,n_{g}\\}$ (the nodes that are not reduced) and
$I_{c}=\\{n_{g}+1,...,n_{g}+n_{c}\\}$ (the ones that are reduced), we write
the elements of Eq. (2) in block form (reordering indices if necessary)

$\displaystyle\bm{x}=\left(\begin{array}[]{c}\bm{x}^{g}\\\
\bm{x}^{c}\end{array}\right)\,,\quad\bm{\omega}=\left(\begin{array}[]{c}\bm{\omega}^{g}\\\
\bm{\omega}^{c}\end{array}\right)\,,\quad\bm{L}=\left(\begin{array}[]{cc}\bm{L}^{gg}&\bm{L}^{gc}\\\
\bm{L}^{cg}&\bm{L}^{cc}\end{array}\right)\,.$ (16)

The Kron-reduced coupling matrix is then the Schur complement of $\bm{L}^{cc}$
in $\bm{L}$,

$\displaystyle\bm{L}^{r}=\bm{L}^{gg}-\bm{L}^{gc}(\bm{L}^{cc})^{-1}\bm{L}^{cg}\,,$
(17)

and applying a similar reduction to the vector of natural velocities gives

$\displaystyle\bm{\omega}^{r}=\bm{\omega}^{g}-\bm{L}^{gc}(\bm{L}^{cc})^{-1}\bm{\omega}^{c}\,.$
(18)

This yields the Kron-reduced version of Eq. (2), which is restricted to the
non-reduced nodes

$\displaystyle\bm{\omega}^{r}=\bm{L}^{r}\bm{x}^{g}\,,$ (19)

which allows to solve the equations for the subset of nodes $I_{g}$.

Note that if the coupling matrix is Laplacian, then its Kron-reduced
counterpart is also Laplacian.

##  3 Locating a single disturbance

We assume the disturbance to vary sufficiently slowly, meaning that its effect
will propagate throughout the network. The detection of the disturbance can
then be performed at any point of the system (e.g., through Fourier transform
of the signal), but locating it requires a more elaborate procedure. We use
the same method to identify disturbances at nodes or on lines. In order to
perform this we need to know the Jacobian matrix $\bm{L}$ of the system,
possibly Kron-reduced if some nodes are not measured, i.e., we need to be
informed about the unperturbed state of the system. Under our assumption of
relatively small amplitude in ${\bm{\omega}}$, a reasonable approximation of a
stable fixed point of Eq. (1) is

$\displaystyle\bm{x}^{g}=(\bm{L}^{r})^{\dagger}\bm{\omega}_{r}\,,$ (20)

where $\dagger$ denotes the _Moore-Penrose pseudoinverse_ of the matrix [29,
Sec. 5.5.2]. Note that it then satisfies
$\bm{L}\bm{L}^{\dagger}=\mathbb{I}-n^{-1}\bm{1}\bm{1}^{\top}$ because the Kron
reduction $\bm{L}^{r}$ of $\bm{L}$ is a Laplacian matrix.

The effect of the disturbances we consider translates into time-varying
(reduced) Jacobian matrix $\tilde{\bm{L}}^{r}(t)$ and/or velocity vector
$\tilde{\bm{\omega}}^{r}(t)$. As we assumed the amplitude of the disturbance
to be relatively small, the evolution of the state of the system is reasonably
approximated as
$\bm{x}^{g}(t)\approx[\tilde{\bm{L}}^{r}(t)]^{\dagger}\tilde{\bm{\omega}}^{r}(t)$.
If we have access to the time series of $\bm{x}^{g}(t)$ and the knowledge of
$\bm{L}^{r}$, we can compute

$\displaystyle\bm{\psi}(t)=\bm{L}^{r}\bm{x}^{g}(t)=\bm{L}^{r}\left[\tilde{\bm{L}}^{r}(t)\right]^{\dagger}\tilde{\bm{\omega}}^{r}(t)\,,$
(21)

which we refer to as the _frequency mismatch_. We show now that the time
series of $\bm{\psi}(t)$ are able to locate the disturbed element with as much
accuracy as can be expected.

We distinguish three cases: a disturbance at a node (which can be reduced or
not); a line disturbance between two nodes that are not Kron-reduced; and a
line disturbance with at least one end-node that is Kron-reduced.

###  3.1 Nodal disturbance

This case is the easiest to treat. The Jacobian matrix is constant in time
$\tilde{\bm{L}}^{r}(t)=\bm{L}^{r}$ and only one component of the vector of
velocities is time-varying. There are now two possible cases. If the perturbed
node is not reduced, we can locate it exactly. Indeed, in this case, the
reduced velocity vector has the form

$\displaystyle\tilde{\bm{\omega}}^{r}(t)=\tilde{\bm{\omega}}^{g}(t)-\bm{L}^{gc}(\bm{L}^{cc})^{-1}\bm{\omega}^{c}=\bm{\omega}^{r}+\xi_{\rm
n}(t)\bm{e}_{i}\,,$ (22)

which, once plugged into Eq. (21), yields

$\displaystyle\bm{\psi}(t)=\bm{\omega}^{r}+\xi_{\rm
n}(t)\left(\bm{e}_{i}-n^{-1}\bm{1}\right)\,.$ (23)

Provided that $n$ is not too small, the amplitude of $\bm{\psi}(t)$ will be
significantly larger at the perturbed node, allowing then to identify it.

If the perturbed node is Kron-reduced, the reduced velocity vector is
expressed as

$\displaystyle\tilde{\bm{\omega}}^{r}(t)=\bm{\omega}^{g}-\bm{L}^{gc}(\bm{L}^{cc})^{-1}\tilde{\bm{\omega}}^{c}(t)\,.$
(24)

Again, plugging this into Eq. (21), we see that all non-reduced nodes that are
neighbors of the reduced component are significantly impacted by the
disturbance. We can then identify the area of the network in which the faulty
node is located, but we cannot identify the node exactly. This is not
surprizing as we do not have a good enough accuracy in the measurements.

###  3.2 Line disturbance between non-reduced end-nodes

This is the most interesting case because we are able to locate the faulty
line exactly. It also covers the case where no Kron reduction is applied. Here
the faulty line has no influence on the reduced velocity vector and is a
rank-1 perturbation of the reduced Jacobian matrix,

$\displaystyle\tilde{\bm{\omega}}^{r}(t)=\bm{\omega}^{r}\,,\qquad\tilde{\bm{L}}^{r}(t)=\bm{L}^{r}+\xi_{\rm
l}(t)\bm{e}_{ij}\bm{e}_{ij}^{\top}\,,$ (25)

where we adapted the size of $\bm{e}_{ij}$ and the indices $i$ and $j$ in
order to comply with the Kron reduction. Introducing Eq. (25) into the
computation of the frequency mismatch yields

$\displaystyle\bm{\psi}(t)$ $\displaystyle=\bm{L}^{r}\left[\bm{L}^{r}+\xi_{\rm
l}(t)\bm{e}_{ij}\bm{e}_{ij}^{\top}\right]^{\dagger}\bm{\omega}^{r}$ (26)
$\displaystyle=\bm{L}^{r}\left[(\bm{L}^{r})^{\dagger}-\frac{\xi_{\rm
l}(t)(\bm{L}^{r})^{\dagger}\bm{e}_{ij}\bm{e}^{\top}_{ij}(\bm{L}^{r})^{\dagger}}{1+\xi_{\rm
l}(t)\bm{e}_{ij}^{\top}(\bm{L}^{r})^{\dagger}\bm{e}_{ij}}\right]\bm{\omega}^{r}$
$\displaystyle=\bm{\omega}^{r}-\alpha(t)\left[\bm{e}_{ij}^{\top}(\bm{L}^{r})^{\dagger}\bm{\omega}^{r}\right]\bm{e}_{ij}\,,$

where we used Sherman-Morrison Eq. (9) at the second line and defined

$\displaystyle\alpha(t)=\frac{\xi_{\rm l}(t)}{1+\xi_{\rm
l}(t)\bm{e}_{ij}^{\top}(\bm{L}^{r})^{\dagger}\bm{e}_{ij}}\,,$ (27)

where we used that $\bm{L}\bm{L}^{\dagger}\bm{e}_{ij}=\bm{e}_{ij}$ which
follows from the property of the pseudoinverse discussed below Eq. (20) . One
sees that the only time varying components of $\bm{\psi}(t)$ are at the two
end points of the perturbed line, which allows to locate it exactly.

This case is illustrated in Fig. 1(d), which shows the time series of
$\psi_{i}(t)$ for nodes $1$ to $9$ of the network shown in Fig. 1(a), when the
coupling of the green line is varying with time. One sees that the perturbed
line is unambiguously identified.

Figure 1: Example network with the five possible types and location of
disturbances in Kron-reduced network. We consider an oscillating disturbance,
whose frequency is much smaller than the time scales of the network. (a)
Network before Kron reduction. The three square nodes $\\{10,11,12\\}$ are not
measured and thus are reduced. The colored nodes are subject to disturbances
[panels (b) and (c)], as well as the three colored lines [panels (d), (e), and
(f)]. (b-f) Time series of $\bm{\psi}(t)$ as defined in Eq. (21). The
correspondence between indices and colors is given by
$(1,2,3,4,5,6,7,8,9)\leftrightarrow({\rm
blue,orange,green,red,purple,brown,pink,gray,yellow})$. (b) Disturbance at the
non-reduced node $6$. (c) Disturbance at the reduced node $11$. (d)
Disturbance on line $(7,8)$. None of the end-nodes are reduced. (e)
Disturbance on line $(2,10)$ where node $10$ is reduced. (f) Disturbance on
line $(11,12)$, where both nodes are reduced.

Remark. We notice that the larger the weight of the disturbed line (in the
Jacobian matrix ${\bm{L}}$), the more accurate the identification. Indeed, if
the value $L_{ij}$ is small, its variation will hardly be noticed on the
dynamics of the agents. This means that our method is the most reliable for
disturbances occurring on the most important lines of the system.

###  3.3 Line disturbance with at least one reduced end-node

In this situation, the exact location of the perturbed line cannot be
determined based on the measurements because at least one of the end points is
hidden to the observer. We detail the case where the two end points of the
faulty line are in the Kron-reduced set. The mixed case where only one of the
end points is reduced is very similar. All of our computations are inspired
from [19, Secs. V.1 and V.2] and rely on an appropriate use of the Sherman-
Morrison formula Eq. (9) and of the formulae for the Kron reduction Eqs. (17)
and (18).

Assuming that $i$ and $j$ are in the Kron-reduced set, one gets time-varying
reduced velocity vector and time-varying reduced Jacobian matrix

$\displaystyle\tilde{\bm{\omega}}^{r}(t)=\bm{\omega}^{r}-\beta(t)\left[\bm{e}_{ij}^{\top}(\bm{L}^{cc})^{-1}\bm{\omega}^{c}\right]\bm{v}\,,$
(28)
$\displaystyle\tilde{\bm{L}}^{r}(t)=\bm{L}^{r}-\beta(t)\bm{v}\bm{v}^{\top}\,,$
(29)

where we defined

$\displaystyle\beta(t)=\frac{\xi_{\rm l}(t)}{1+\xi_{\rm
l}(t)\bm{e}_{ij}^{\top}(\bm{L}^{cc})^{-1}\bm{e}_{ij}}\,,$ (30)

and

$\displaystyle\bm{v}=\bm{L}^{gc}(\bm{L}^{cc})^{-1}\bm{e}_{ij}\,,$ (31)

$\bm{e}_{ij}$ is adapted according to the Kron reduction, and we used that
$\bm{L}$ is symmetric, i.e., $\bm{L}^{gc}=(\bm{L}^{cg})^{\top}$. Now in the
same spirit as before, the frequency mismatch is

$\displaystyle\bm{\psi}(t)$
$\displaystyle=\bm{L}^{r}[\tilde{\bm{L}}^{r}(t)]^{\dagger}\bm{\omega}^{r}(t)$
(32)
$\displaystyle=\bm{L}^{r}\left[(\bm{L}^{r})^{\dagger}-\frac{\beta(t)(\bm{L}^{r})^{\dagger}\bm{v}\bm{v}^{\top}(\bm{L}^{r})^{\dagger}}{1+\beta(t)\bm{v}^{\top}(\bm{L}^{r})^{\dagger}\bm{v}}\right]$
$\displaystyle\phantom{=}\times\left[\bm{\omega}^{r}-\beta(t)\left(\bm{e}_{ij}^{\top}(\bm{L}^{cc})^{-1}\bm{\omega}^{c}\right)\bm{v}\right]$
$\displaystyle=\bm{\omega}^{r}+\gamma(t)\bm{v}\,,$

where again, we used the Sherman-Morrison formula, Eq. (9), and where we
gathered all time-varying quantities in $\gamma(t)$.

For the case where node $i$ is not reduced and node $j$ is, a similar
calculation gives

$\displaystyle\bm{\psi}(t)=\bm{\omega}^{r}+\tilde{\gamma}(t)\tilde{\bm{v}}\,,$
(33)

where

$\displaystyle\tilde{\bm{v}}=\bm{e}_{i}+\bm{L}^{gc}(\bm{L}^{cc})^{-1}\bm{e}_{j}\,,$
(34)

with $\bm{e}_{i}$ and $\bm{e}_{j}$ being adapted to the Kron-reduced system,
and where, again, we assembled all time-varying quantities in
$\tilde{\gamma}(t)$.

In these two cases, the effect of the perturbed line will be measured [through
$\bm{\psi}(t)$] on the non-reduced nodes connected to the reduced component
containing the perturbed line. This can be verified by inspection of the
submatrix $\bm{L}^{gc}$ and is seen in Figs. 1(e) and (f), showing the time
series of $\psi_{i}(t)$, when the disturbed line is at the blue and orange
lines respectively. One cannot identify precisely its location, but it is
possible to identify the reduced component to which it belongs, which is the
most we could expect from the available measurements. In Fig. 1(e), where the
perturbed line connects a reduced node to a non-reduced one, we observe that
the amplitude of the frequency mismatch $\bm{\psi}(t)$ is much larger at the
non-reduced end of the disturbed line than at any other non-reduced node. This
makes sense if line couplings are all similar, the $i$th component of
$\tilde{\bm{v}}$ [Eq. (34)] is likely to be significantly larger than its
other components, which will translate as a larger amplitude of variation in
$\psi_{i}(t)$. However, in full generality, we cannot guarantee that the
amplitude of $\psi_{i}(t)$ will be larger than all other components of the
frequency mismatch.

In summary, by simply computing the frequency mismatch, we are able to locate
node and line disturbances, at locations in the system that are not
necessarily measured.

##  4 Multiple disturbances

Interestingly, we observe that our approach allows to locate multiple
disturbances acting simultaneously, provided they have sufficiently different
characteristics so one is able to distinguish them. We treat here the case of
two simultaneous disturbances, the case with more disturbances extends
naturally from it. We show the cases where the disturbances occur in non-
reduced areas. The reduced case works similarly as in the single disturbance
case.

For two nodal disturbances $\xi_{1}(t)$ and $\xi_{2}(t)$, at non-reduced nodes
$i$ and $j$ respectively, the calculation of the frequency mismatch yields,

$\displaystyle\bm{\psi}(t)=\bm{\omega}^{r}+\xi_{1}(t)(\bm{e}_{i}-n^{-1}\bm{1})+\xi_{2}(t)(\bm{e}_{j}-n^{-1}\bm{1})\,.$
(35)

Provided the disturbances $\xi_{1}$ and $\xi_{2}$ have sufficiently different
characteristics, one can identify their location if $n$ is not too small. If
one of the nodes is reduced, only the area of the disturbance can be located.

Two simultaneous line disturbances can be modelled as the following rank-2
perturbation of the Jacobian matrix,

$\displaystyle\tilde{\bm{L}}(t)=\bm{L}+\left(\begin{array}[]{cc}\bm{e}_{ij}&\bm{e}_{k\ell}\end{array}\right)\cdot\left(\begin{array}[]{c}\xi_{1}(t)\bm{e}_{ij}^{\top}\\\
\xi_{2}(t)\bm{e}_{k\ell}^{\top}\end{array}\right)=\bm{L}+\bm{U}\cdot\bm{V}_{t}\,,$
(39)

where lines $(i,j)$ and $(k,\ell)$ are perturbed. As in the single disturbance
case, one can then calculate ${\bm{\psi}}(t)$, which yields, using the
Sherman-Morrison-Woodbury formula Eq. (8),

$\displaystyle\bm{\psi}(t)$
$\displaystyle=\bm{L}[\tilde{\bm{L}}(t)]^{\dagger}\bm{\omega}=\bm{L}\left[\bm{L}^{\dagger}-\bm{L}^{\dagger}\bm{U}\left(\mathbb{I}+\bm{V}_{t}\bm{L}^{\dagger}\bm{U}\right)^{-1}\bm{V}_{t}\bm{L}^{\dagger}\right]\bm{\omega}$
(40)
$\displaystyle=\bm{\omega}-\bm{U}\underbrace{\left(\mathbb{I}+\bm{V}_{t}\bm{L}^{\dagger}\bm{U}\right)^{-1}}_{=\bm{X}}\bm{V}_{t}\bm{\theta}^{*}\,.$

Recall that $\bm{\omega}$, the columns of $\bm{U}$, and the row of
$\bm{V}_{t}$ are orthogonal to the kernel of $\bm{L}$.

Parsing Eq. (40), we see that

$\displaystyle{\bm{V}_{t}}\bm{\theta}^{*}=\left(\begin{array}[]{c}\xi_{1}(t)(\theta_{i}^{*}-\theta_{j}^{*})\\\
\xi_{2}(t)(\theta_{k}^{*}-\theta_{\ell}^{*})\end{array}\right)\,.$ (43)

Furthermore, we see numerically that the matrix $\bm{X}$ is close to be
diagonal. Namely, the off-diagonal elements of $\bm{X}$ are smaller (in
absolute value) than its diagonal elements by at least one order of magnitude.
This yields

$\displaystyle\bm{\psi}(t)\approx\bm{\omega}-X_{11}\xi_{1}(t)(\theta_{i}^{*}-\theta_{j}^{*})\bm{e}_{ij}-X_{22}\xi_{2}(t)(\theta_{k}^{*}-\theta_{\ell}^{*})\bm{e}_{k\ell}\,,$
(44)

where $X_{pp}$ is the $p$th diagonal element of $\bm{X}$. The disturbances can
then be located unequivocally, provided they have sufficiently different
characteristics.

Figure 2: (a)-(c) 2nd, 10th and 20th eigenvectors of the European electric
transmission grid. (d) All of its eigenvalues, which vary over almost six
order of magnitudes. Purple and black dots display the locations of faulty
lines in first and second columns of Fig. 4 respectively.

Even though we have not been able to certify analytically that $\bm{X}$ is
almost diagonal, we explain it heuristically as follows. Analytically,
everything boils down to showing that $\bm{V}_{t}\bm{L}^{\dagger}\bm{U}$ is
(close to) diagonal. Let us forget the time dependence of $V$ for now
($\xi_{1}=\xi_{2}\equiv 1$) and compute the $(\alpha,\beta)$th term of this
matrix, where we denote by $(i_{\alpha},j_{\alpha})$ the disturbed lines,

$\displaystyle(\bm{V}\bm{L}^{\dagger}\bm{U})_{\alpha\beta}=\bm{e}_{i_{\alpha}j_{\alpha}}^{\top}\bm{L}^{\dagger}\bm{e}_{i_{\beta}j_{\beta}}=\sum_{\gamma\geq
2}\frac{(u_{\gamma,i_{\alpha}}-u_{\gamma,j_{\alpha}})(u_{\gamma,i_{\beta}}-u_{\gamma,j_{\beta}})}{\lambda_{\gamma}}\,,$
(45)

where $\lambda_{\gamma}$ is the $\gamma$th eigenvalue of $\bm{L}$ and
$u_{\gamma,k}$ denotes the $k$th component of the associated eigenvector. At
first sight, one could think that the main contribution to Eq. (45) comes from
eigenmodes with small eigenvalues, i.e., $\gamma$ close to $2$. However it
appears that, for the corresponding eigenvectors, the components of
neighboring nodes are rather similar (see Fig. 2, showing some eigenmodes on a
large-scale network modelling the European power transmission grid). It is
particularly true in large-scale networks where these low laying eigenmodes
spread over the whole network, at odds with high laying ones that are rather
localized on a few nodes. One could think then, that the contribution comes
from these high laying eigenmodes as the components of neighboring nodes might
differ significantly [32]. However, these contributions also remain small due
to the denominator in Eq. (45) where large values for $\lambda_{\gamma}$ bring
down the contribution. The diagonal terms in $\bm{X}$ however are not
vanishing due to the identity matrix appearing in Eq. (40).

We performed a numerical investigation of the matrix
$\bm{V}\bm{L}^{\dagger}\bm{U}$ for the European grid. On average over all
possible pairs of lines, the absolute value of its off-diagonal term is
$5.2\cdot 10^{-6}$ and it reaches $5.9\cdot 10^{-2}$ for the worst
configuration. This seems to validate that the assumption we made on $\bm{X}$
holds, at least for interaction graphs based on infrastructures.

Even though we have not been able to ground our observations in analytical
results, it is remarkable that such a naive approach is able to isolate
various disturbances in a complex system.

##  5 Implementation and numerical validation

In order to validate our approach, we simulated three instanciations of the
dynamics Eq. (1) with three different types of line disturbances, on three
different networks.

###  5.1 Algorithm

We propose here and summarize in Fig. 3 the course of an algorithm that would
implement our method.

First, with the Jacobian matrix $\bm{L}$ and the times series
$\\{\bm{x}(t)\\}_{t=1}^{T}$ as input, one can straightforwardly compute the
frequency mismatch $\\{\bm{\psi}(t)\\}_{t=1}^{T}$. Then, the algorithm needs
to identify groups of outlier nodes, whose frequency mismatch trajectory are
significantly different from the bulk of all measured nodes. On top of that,
these outlier nodes need to be grouped in sets ($I_{1},...,I_{p}$) of
qualitatively similar trajectories, e.g., blue and green (resp. orange and
red) curves are grouped together in Figs. 5 and 6. Defining what
”significantly different” means is not always trivial and might be problem
dependent, but it is necessary in order to identify disturbed nodes.
Furthermore, identifying ”significantly different” curves and classifying them
by similarity might by a challenging problem for data analysts, but these
issues are beyond the scope of our expertise in general and of this manuscript
in particular.

Once the nodes have been classified, each subset of nodes should correspond to
a disturbance in the system. In Fig. 3, we propose a way to determine the type
of disturbance associated to each subset $I_{j}$.

Figure 3: Block diagram of a suggested algorithm implementing our method. Note
that this manuscript covers the first block only, whose task is to ease the
work of downstream blocks. This does not mean that these subsequent blocks are
trivial.

###  5.2 PanTaGruEl

The PanTaGruEl model of the interconnected European grid [32, 33] is built on
publicly available data of geolocalization of power system elements, and
parameters (admittances, productions, loads,…) are reconstructed based on
standard assumptions. It consists of 3809 nodes and 7343 lines. For our
simulations, we consider an oscillating disturbance,

$\displaystyle\xi_{\rm l}(t)=\xi_{0}\sin(\omega_{\rm m}t)\,,$ (46)

satisfying the conditions detailed in Sec. 2.2. From a practical point of
view, this allows us to tune the time scale of the disturbance, which is the
oscillation frequency $\omega_{\rm m}$ in order to guarantee that it is
sufficiently slow with respect to all time scales of the network.

In the following, we compare the disturbance detection using our method
relying on $\bm{\psi}(t)$ or by simply measuring the agents’ frequencies
$\dot{\bm{x}}(t)$ and value $\bm{x}(t)$. The left panels of Fig. 4 shows the
maximal amplitude of $\psi_{i}(t)$ for the oscillating disturbance of five
different lines of the PanTaGruEl model, shown in Fig. 2. Lines are perturbed
one at a time, in different simulations. The right panels show the maximal
amplitude of the frequencies $\dot{x}_{i}(t)$ for the same time series. The
red circle indicate the indices of the two extremities of the disturbed line.
Our method identifies unequivocally the endpoints of the disturbed lines
whereas the time series of frequencies are not able to do so.

Figure 4: Comparison of the detection of perturbed lines based on
$\psi_{i}(t)$, $\dot{x}_{i}(t)$ or $x_{i}(t)$, in the PanTaGruEl, a model of
the European electric transmission grid. The two columns show two
configurations with different malfunctioning elements. In each panel, the end
nodes of the two disturbed lines are indicated by the red circles are their
locations are indicated in Fig. 2. Reasonable noise was added at each node in
order to challenge the robustness of our approach. The disturbed lines are
exactly identified as the ones whose end nodes have a significant amplitude
for $\psi_{i}$.

In the case of sparse (incomplete) measurements, we have tested, but not shown
in this works, our method for faults where both ends of the perturbed line
belong to the reduced (non-measured) nodes. In this case, our methods
correctly detect the nodes in the vicinity of the faulty elements. However, at
this task, our method is not significantly better than the detection
proceeding by inspection of the trajectories of $x_{i}(t)$ , $\dot{x}_{i}(t)$
.

###  5.3 Euroroad and US Airports

In order to illustrate the variety of couplings that our method covers, we
applied it to:

  * •

A set of $n=1039$ agents linearly coupled according to the largest component
of the European road network from [34, 35]. Two arbitrary lines are pertubed
by a noisy signal. Results are shown in Fig. 5;

  * •

A set of $n=1572$ agents interacting through a _higher-order Kuramoto
coupling_ [36] (couplings of order up to $q=3$ are used here), with
interaction structure given by the largest connected component of the network
of US airports from [37, 38]. One line is perturbed by a ramp signal and
another one is subject to a step signal. Results are shown in Fig. 6.

In each case, we added a reasonable amount of noise at each node, for sake of
generality.

Figure 5: Phases, velocities, and frequency mismatch for the European road
network [34, 35] with linear coupling and two noisy lines. The pairs of
trajectories blue-green and orange-red correspond to the end nodes of each
disturbed lines. The grey area covers the same quantities for all other nodes.
The disturbances cannot be seen in the phases and velocities trajectories,
whereas they are clearly identified and located with the frequency mismatch.
The panels on the right side are the snapshot at the end of the time series of
the left side.  Figure 6:  Phases, velocities, and frequency mismatch for the
US airports network [37, 38] with third order Kuramoto coupling [36]. One line
is disturbed by a step signal and the other by a ramp signal. The pairs of
trajectories blue-green and orange-red correspond to the end nodes of each
disturbed lines. The grey area covers the same quantities for all other nodes.
The disturbances cannot be seen in the phases and velocities trajectories,
whereas they are clearly identified and located with the frequency mismatch.
The panels on the right side are the snapshot at the end of the time series of
the left side.

In both Figs. 5 and 6, the left panels show the time series of $x_{i}(t)$,
$\dot{x}_{i}(t)$, and $\psi_{i}(t)$ from top to bottom. The right panels show
the snapshot of the same quantities at the end of the time series on the left,
at $t=3.0$.

We acknowledge that, even though it is straighforward to identify the end-
nodes of the disturbed lines by visual inspection of Figs. 5 and 6, it might
be an algorithmic challenge to automatize this process.

##  6 Conclusion

We proposed an elegant method to identify and locate disturbances in networked
dynamical systems. Our method relies on time series of the agents degree of
freedom and is able to differentiate between line disturbances and nodal
disturbances. In the case of partial measurements, we are able to locate
precisely the perturbation if the nodes where the it occurs are measured.
Otherwise, we can determine the area of the network where the faulty,
unmeasured nodes/lines are located, using the Kron-reduced network. The main
condition required for our method to work out well is that the disturbance
changes much slower than the intrinsic time scales of the system. However it
is not restricted to a particular shape for the perturbation or a specific
dynamical model, making it rather widely applicable. Interestingly, our method
is also able to locate multiple disturbance occurring at the same time, given
that they have different characteristics. This represents a remarkable feature
as in large scale networks made of thousands or millions of nodes, many
disturbances are likely to overlap in the time series.

We believe that the idea raised in this manuscript will help the development
of more efficient tools to detect and locate disturbances more accurately.
Future work will aim at proposing a scheme for online detection and
localization of disturbances.

## Acknowledgments

RD and MT were supported by the Swiss National Science Foundation under grant
no. 200020_182050. RD was supported by ETH Zurich funding and by the Swiss
National Science Foundation under grant no. P400P2_194359.

##  Appendix A Directed networks

Although all mathematical objects involved in the derivations of Sec. 3 are
well-defined for directed graphs, the results are harder to desintricate. With
directed interactions, the nodal perturbation is similar to the undirected
case,

$\displaystyle\tilde{\bm{\omega}}(t)=\bm{\omega}+\xi_{\rm n}(t)\bm{e}_{i}\,,$
(47)

and the line disturbance has the form

$\displaystyle\tilde{\bm{L}}(t)=\bm{L}+\xi_{\rm l}(t)\bm{e}_{i}\bm{e}_{ij}\,.$
(48)

In any case, the computation $\bm{\psi}(t)$ involves a vector of the form

$\displaystyle\bm{w}=\bm{L}\bm{L}^{\dagger}\bm{e}_{i}\,,$ (49)

which is the othogonal projection of $\bm{e}_{i}$ on the image of $\bm{L}$
[29, Sec. 5.5.2]. One can verify that $\bm{w}$ highly depends on the position
of node $i$ in the network and the structure thereof. Desintricating the
effect of the network structure on the vector $\bm{w}$ is beyond the scope of
this manuscript and we defer it to future work.

## References

## References

  * [1] S. H. Strogatz.  Exploring complex networks.  Nature, 410(6825):268–276, 2001.
  * [2] F. A. Rodrigues, T. K. DM. Peron, P. Ji, and J. Kurths.  The Kuramoto model in complex networks.  Phys. Rep., 610:1–98, 2016.
  * [3] C. Castellano, S. Fortunato, and V. Loreto.  Statistical physics of social dynamics.  Rev. Mod. Phys., 81:591–646, May 2009.
  * [4] A.-L. Barabási.  Network science.  Cambridge University Press, Cambridge, England, 2016.
  * [5] T. Stankovski, T. Pereira, P V. E. McClintock, and A. Stefanovska.  Coupling functions: Universal insights into dynamical interaction mechanisms.  Rev. Mod. Phys., 89:045001, Nov 2017.
  * [6] C. D. Brummitt, P. D. H. Hines, I. Dobson, C. Moore, and R. M. D’Souza.  Transdisciplinary electric power grid science.  Proc. Nat. Acad. Sci. USA, 110(30):12159, 2013.
  * [7] J. Markard.  The next phase of the energy transition and its implications for research and policy.  Nature Energy, 3:628–633, 2018.
  * [8] L. Dosiek, N. Zhou, J. W. Pierre, Z. Huang, and D. J Trudnowski.  Mode shape estimation algorithms under ambient conditions: A comaprative review.  IEEE Trans. Power Syst., 28(2):779–787, 2013.
  * [9] A. Y. Lokhov, M. Vuffray, D. Shemetov, D. Deka, and M. Chertkov.  Online learning of power transmission dynamics.  In IEEE Power Systems Computation Conference (PSCC), pages 1–7. IEEE, 2018.
  * [10] S. Furutani, C. Takano, and M. Aida.  Network resonance method: Estimating network structure from the resonance of oscillation dynamics.  IEICE Trans. Commun., E102-B(4):799–809, 2019.
  * [11] T. R. Nudell and A. Chakrabortty.  A graph-theoretic algorithm for disturbance localization in large power grids using residue estimation.  In 2013 Amercian Control Conference (ACC), pages 3467–3472. IEEE, 2013.
  * [12] S. Upadhyaya and S. Mohanty.  Power quality disturbance localization using maximal overlap discrete wavelet transform.  In Annual IEEE India Conference (INDICON). IEEE, 2015.
  * [13] A. T. Mathew and M. N. Aravind.  PMU based disturbance analysis and fault localization of a large grid using wavelets and list processing.  In IEEE Region 10 Conference (TENCON), pages 879–883. IEEE, 2016\.
  * [14] H.-W. Lee, J. Zhang, and E. Modiano.  Data-driven localization and estimation of disturbance in the interconnected power system.  In IEEE International Conference on Communications, Control, and Computing Techologies for Smart Grids (SmartGridComm). IEEE, 2018.
  * [15] A. Semerow, S. Horn, B. Schwarz, and M. Luther.  Disturbance localization in power systems using wide area measurement systems.  In IEEE International Conference on Power System Technology (POWERCON). IEEE, 2016.
  * [16] S. K. Chakravarthy.  Nonlinear oscillations due to spurious energisation of transformers.  IEE Proceedings-Electric Power Applications, 145(6):585–592, 1998\.
  * [17] J. Horak.  A review of ferroresonance.  In Proceedings of the 57th Annual Conference for Protective Relay Engineers, pages 1–29, 2004.
  * [18] S. Hassan, M. Vaziri, and S. Vadhva.  Review of ferroresonance in power distribution grids.  In Proceedings of the International Conference on Information Reuse and Integration, pages 444–448. IEEE, 2011.
  * [19] T. Coletta and P. Jacquod.  Performance measures in electric power networks under line contingencies.  IEEE Trans. Control Netw. Syst., 7(1):221–231, 2020.
  * [20] R. Delabays, M. Tyloo, and P. Jacquod.  Rate of change of frequency under line contingencies in high voltage electric power networks with uncertainties.  Chaos, 29:103130, 2019.
  * [21] S. M. Brahma.  Fault Location in Power Distribution System With Penetration of Distributed Generation.  IEEE Trans. Power Del., 26(3):1545–1553, 2011.
  * [22] S. Soltan, A. Loh, and G. Zussman.  Analyzing and Quantifying the Effect of k-line Failures in Power Grids.  IEEE Trans. Control Netw. Syst., 5(3):1424–1433, 2017.
  * [23] S. Soltan, M. Yannakakis, and G. Zussman.  Power Grid State Estimation Following a Joint Cyber and Physical Attack.  IEEE Trans. Control Netw. Syst., 5(1):499–512, 2018.
  * [24] S. Soltan and G. Zussman.  EXPOSE the Line Failures Following a Cyber-Physical Attack on the Power Grid.  IEEE Trans. Control Netw. Syst., 6(1):451–461, 2019.
  * [25] S. Soltan, M. Yannakakis, and G. Zussman.  REACT to Cyber Attacks on Power Grids.  IEEE Trans. Netw. Sci. Eng., 6(3):459–473, 2019.
  * [26] M. Jamei, A. Scaglione, and S. Peisert.  Low-Resolution Fault Localization Using Phasor Measurement Units with Community Detection.  In 2018 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm), pages 1–6. IEEE, 2018.
  * [27] M. Jamei, R. Ramakrishna, T. Tesfay, R. Gentz, C. Roberts, A. Scaglione, and S. Peisert.  Phasor Measurement Units Optimal Placement and Performance Limits for Fault Localization.  IEEE J. Sel. Areas Commun., 38(1):180–192, 2020.
  * [28] J. J. Grainger and W. D. Stevenson.  Power System Analysis.  McGrow-Hill, 1st edition, 1994.
  * [29] G. H. Golub and C. F Van Loan.  Matrix Computations.  4th edition, 2013.
  * [30] J. Machowski, J. Bialek, and J. R. Bumby.  Power system dynamics: stability and control.  John Wiley & Sons, 2nd edition, 2008.
  * [31] F. Dörfler and F. Bullo.  Kron reduction of graphs with applications to electrical networks.  IEEE Trans. Circ. Syst., 60(1):150–163, 2013.
  * [32] L. Pagnier and P. Jacquod.  Inertia location and slow network modes determine disturbance propagation in large-scale power grids.  PloS ONE, 14(3), 2019.
  * [33] M. Tyloo, L. Pagnier, and P. Jacquod.  The key player problem in complex oscillator networks and electric power grids: Resistance centralities identify local vulnerabilities.  Science Advances, 5(11):eaaw8359, 2019.
  * [34] L. Šubelj and M. Bajec.  Robust network community detection using balanced propagation.  Eur. Phys. J. B, 81:353–362, 2011.
  * [35] T. P. Peixoto.  The netzschleuder network catalogue and repository.  https://networks.skewed.de/, 2020.  Accessed: 2020-12-16.
  * [36] P. S. Skardal, E. Ott, and J. G. Restrepo.  Cluster synchrony in systems of coupled phase oscillators with higher-order coupling.  Phys. Rev. E, 84:036208, 2011.
  * [37] R. Guimera, S. Mossa, A. Turtschi, and L. A. N. Amaral.  The worldwide air transportation network: Anomalous centrality, community structure, and cities’ global roles.  Proc. Nat. Acad. Sci. USA, 102(22):7794–7799, 2005.
  * [38] V. Colizza, R. Pastor-Satorras, and A. Vespignani.  Reaction–diffusion processes and metapopulation models in heterogeneous networks.  Nature Phys., 3(4):276–282, 2007.


# Modeling of porous lithium metal electrodes: turning the Li-dendrite problem
around

Giovanna Bucci
Palo Alto Research Center Inc.
Palo Alto, CA, 94304 USA
<EMAIL_ADDRESS>
Tushar Swamy
Our Next Energy Inc.
Novi, MI, 29050 USA
&W. Craig Carter
Massachusetts Institute of Technology
Cambridge, MA, 02139 USA
&Morad Behandish
Palo Alto Research Center Inc.
Palo Alto, CA, 94304 USA

###### Abstract

The properties of rechargeable lithium-ion batteries are determined by the
electrochemical and kinetic properties of their constituent materials as well
as by their underlying microstructure. Microstructural design can be leveraged
to achieve a leap in performance and durability. Here we investigate a porous
electrode structure, as a strategy to increase the surface area, and provide
structural stability for Li-metal anodes. The porous architecture consists of
a mixed electron/ion conductor that function as a scaffold for lithium metal
deposition. A new finite element model was developed to simulate the large
topological changes associated with Li plating/stripping. This model is used
to predict the current density distribution as a function of material and
structural properties. A dimensionless quantity that combines Li-ion
conductivity, surface impedance and average pore size is shown to be a good
indicator to predict the peak current density. Preventing current localization
at the separator reduces the risk of cell shorting. The analyses show that the
peak current scales as $(hG)^{1/2}$, where $h$ is the ratio between surface
and bulk conductivity and $G$ is the average pore size. Stability analyses
suggest that the growth is morphologically stable, and that confining Li-
plating into pores can enable high-energy density solid-state batteries. In
addition to optimizing porous electrodes design, this finite element method
can be extended to studying other Li-battery structures.

_K_ eywords solid state battery $\cdot$ lithium metal anode $\cdot$ porous
architecture $\cdot$ mesoscale modeling $\cdot$ mixed electron/ion conductor

##  1 Introduction

The need for fast charging, high energy density, long battery life is
currently driving the search for a robust solid-state electrolyte which can
enable Li metal electrodes. The attempts in using lithium as a negative
electrode in Li-ion batteries have revived in recent years, as the energy
density of graphite-based cells is reaching its limits, and the demand for
high-energy applications is growing. High-rate, high-capacity cycling of Li
metal anodes has been hindered by cell-shorting due to non-uniform Li
deposition (forming a dendritic structure), causing the solid electrolyte (SE)
layer to fracture [1, 2, 3, 4, 5, 6]. To address this problem, we investigate
how a porous electrode structure can self-regulate the local current density
and prevent cell shorting. The porous topology provides a scaffold for Li-
deposition and stripping, supporting both mechanical integrity and Li-
accessibility. The goal of this study is to understand what characteristics
(geometry and material properties) are required to achieve stable cycling at
current densities relevant to vehicle electrification.

To date, several strategies have been developed to mitigate the dendritic
growth of metallic Li. In particular, Li metal storage in three-dimensional
(3D) framework electrodes is considered an effective dendrite growth
mitigation strategy since these electrodes with large surface areas and high
porosity result in reduced local Li plating current densities and further can
confine metallic Li within their porous architectures resulting in increased
tolerance to volume expansion and associated particle fracture effects [1, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. Promising porous
electrode strategies include 3D CuZn current collectors [7], Cu current
collector and Li composite anodes [8], porous carbon electrodes embedded with
high surface area N-doped carbon crystals [9], and with lithiophilic ZnO
quantum dots [10]. Effective porous electrode strategies for all-solid-state
systems include 3D electron-ion conducting frameworks such as porous garnet
electrolytes with carbon nanotube coatings [11], and 3D printing enabled
porous garnet microstructures [11, 18]. Most of these strategies, however,
exhibit uneven Li plating-stripping because the ionic resistance of the
electrolyte filled pores leads to preferential plating-stripping closer to the
separator interface [13]. To address this issue, recent research has focused
on (i) microstructural engineering and (ii) lithiophillic surface treatment,
in order to facilitate more uniform plating-stripping across the entire
surface area of the porous electrode [14, 15, 16, 17, 19].

These examples of mesoscale architecture motivate the development of a
modeling framework to study the relationship between structure, performance
and durability in porous electrodes. Modeling of energy storage materials has,
so far, focused on composite electrodes [22, 23, 24, 25, 26, 27, 28, 29, 30].
Here we present a computational method developed to track material growth
driven by electro-deposition within a complex geometry. Employing a phase
field method to track the Li interface is a possible approach [31, 32, 33, 34,
35, 36, 37, 38, 39, 40]. However, we expect Li-metal regions to not simply
merge, but to come in contact, deform and exert stress on the electrode walls.
For this reason, we choose to implement a finite element method with a sharp
interface representation, and we anticipate adding stress and fracture in
future work. We also anticipate that as pores begin to fill, there will be
additional terms in the Li diffusion potential. Additionally, the volume
change as lithium comes out of solution and plates within the pores creates
difficulties in the phase field method.

##  2 Methods

Figure 1: a) Cartoon representation of the anode porous architecture
consisting of a mixed electron/ion conductor that function as a scaffold for
lithium metal deposition. At discharge, Li progressively fills the pores close
to the separator (blue region on the right), and over time the deposition
front moves towards the current collector. b) A representative section of the
lattice with a Voronoi structure is modeled with a finite element approach.
The first part of the algorithm calculates the deposition rate at the inner
boundary (marked in black). A fixed current is applied at the separator
interface (blue line) and fixed voltage is imposed at the interface with
current-collector (red line on the left). Zero-flux boundary condition is
imposed at the top and bottom edge of the domain.

A finite element model was implemented to simulate electro-deposition in a
porous electrode structure. In the simulation, the deposition front moves as
the cells is the microstructure are filled with lithium (Fig. 1). We assume
that the Li plates at the Li/Vapor interface and not the Li/electrolyte
interface. The cathode overpotential is the driving force for lithium
deposition and it is calculated by solving the electrochemical problem
described below. We consider the domain sketched on the left of Fig. 1 as
representing a porous solid mixed electron/ion conducting material. Steady-
state conduction of Li-ion across the electrode is modeled by the Laplace
equation as in Eq. 1a, with the potential of lithium $\phi(x,y)$ treated as
independent variable. A fixed current density is imposed at the separator
interface $\partial\mathcal{P}_{+}$ (blue right edge in the figure in Fig. 1).
In Eq. 1c, the constant $\overline{i}$ is the ratio between the current
density and the electrolyte bulk conductivity $\kappa$. The kinetics at the
inner Li/Vapor interface is described by a linear approximation of the Butler-
Volmer equation, with the convective coefficient $h$ representing the ratio
between surface and bulk conductivity. This is the Robin boundary condition in
Eq. 1b and applies to the charge-transfer interface. At the negative
electrode, the potential is held constant at $\phi_{\text{Li}}$ by the
Dirichlet boundary condition Eq. 1d. A zero-flux boundary condition is assumed
for the top and bottom edge of the domain.

$\displaystyle\Delta\phi=0\qquad\text{on }\mathcal{P}$ (1a)
$\displaystyle\nabla\phi\cdot\mathbf{n}=h(\phi-\phi_{Li})\qquad\text{on
}\partial\mathcal{P}_{\text{in}}$ (1b)
$\displaystyle\nabla\phi\cdot\mathbf{n}=\overline{i}\qquad\text{on
}\partial\mathcal{P}_{+}$ (1c) $\displaystyle\phi=\phi_{Li}\qquad\text{on
}\partial\mathcal{P}_{-}$ (1d)

A two-dimensional model representation is chosen because it provides
phenomenology that will translate to 3D. For a square lattice, in the limit of
an infinitely thin wall, the surface area $S$ scales linearly with the number
of cells in each dimension as $S=6n_{x}n_{y}n_{z}$. The current density
distribution changes with the distance from the separator (x axis in Fig. 1),
and it is on-average uniform in the other two directions.

For this study we consider microstructures derived from a centroidal Voronoi
tessellation, as illustrated in Fig. 1. The pore structure will generally
depend on the fabrication method. Here we choose an isotropic structure with
closed cells, to demonstrate the model capabilities and derive initial scaling
laws.

Electro-deposition is modeled as surface accretion, by assuming that material
growth is normal to the Li/Vapor interface. The numerical implementation of
the model was coded in the Wolfram Language and simulated using Mathematica
111 https://www.wolfram.com/mathematica. Notable functions used from Version
13.0.0 of Mathematica: NDSolve, BSplineFunction, BoundayMesh, ToElementMesh,
BoundaryMeshRegion. Some of the key features of the algorithm are listed
below.

###  2.1 Algorithm

  1. 1.

The first step involves solving Eq. 1a to compute the current density $i$
along the charge-transfer interface $\mathcal{P}_{\text{in}}$. The current
density determine the amount of material plated at the interface. The normal
velocity field is computed as $v_{in}=\Omega_{Li}(i/F)$, with $F$ being the
Faraday constant and $\Omega_{Li}$ the partial molar volume of Li. The left
image below shows electric potential isolines; the right image illustrate the
velocity field at the Li/Vapor interface.

  2. 2.

The product between the normal velocity and the time increment $\Delta t$
determines the new location of the boundary nodes (illustrated in pink on the
left image below). Some regions show high density of points, corresponding to
finite elements with high aspect ratio. Re-sampling of the boundary nodes with
uniform spacing is performed at each time step to prevent element inversion,
and maintain solution accuracy. The re-meshing algorithms is robust with
respect convex/concave geometries and non-uniform current distribution.

  3. 3.

Acute angles or large time steps may cause self-intersection of the boundary.
An algorithm was devised the check for self-intersection and remove the nodes
within the intersecting segment as shown in the two images below

##  3 Results

In Fig. 2 a, the sequence of four images illustrates how the deposition front
moves over time, with Li-metal progressively filling the cells of the porous
electrode structure. The cells closer to the separator (on the right) are
filled first. The current distribution depends on the bulk conductivity of the
foam $\kappa$ and on the impedance at the Li/vapor charge-transfer interface
$\rho$. For this example the following values were chosen $\kappa=0.1$ S/m and
$\rho=10\Omega$ cm2 [18]. The relative current distribution through the
electrode thickness is illustrated in Fig. 2 b, as a function of the depth of
discharge (DOD). (At the beginning of lithiation, the state of charge is 1 and
the DOD is 0.) Initially, the local deposition current density reaches the
maximum value of $\approx 7.2\%$ of the applied current, and it decays
exponentially away from the separator (blue points). This is due to the
increase in surface area provided by the porous electrode structure. Around
DOD=0.5, a steady state regime is reached with a linear distribution of the
current within the electrode thickness (yellow points in Fig. 2 b). The time
for establishing the linear regime depends on the conductivity and interfacial
impedance. For instance, Fig. 3 a shows that a 5-fold increase in interfacial
impedance (from 10 to $50\Omega$ cm2) accelerates the transition to the linear
regime. A moderately larger overpotential at the interface has the effect of
homogenizing the current. However, high interfacial impedance is detrimental
to the battery efficiency.

Figure 2: a) A representative 2D slice of a porous anode, shown in green, is
filled with lithium (gray region) during battery discharge. The porous
structure is made of a mixed electron/Li-ion conducting material with the
following properties bulk Li-ion conductivity $\kappa=1$ mS/m and surface
impedance $\rho=10\Omega$ cm2. The model is able to track the electro-
deposition front inside the pores and predict the evolution in current density
distribution over time. b) Distribution of the relative current density across
the electrode thickness. The color indicates the depth of discharge (DOD)
according to the legend. The maximum value is achieved at the separator
interface, where the horizontal axis indicates a relative electrode thickness
of 1. The current decays away from the separator. The trend changes from
exponential to linear around DOD=0.5.

Ideally if the current is uniformly distributed throughout the microstructure,
the maximum value drops because the entire interface is equally active for Li-
deposition. Decreasing the value of Li-ion bulk conductivity $\kappa$ leads to
higher local currents and tend to concentrate Li-deposition to a few cells at
the time (Fig. 3 b).

The maximum relative current typically occurs at the beginning of lithiation
and it serves as a first metric to evaluate the porous electrode design. To
guide the design of resilient high-energy density systems, we need to identify
what combination of material and structural properties prevents current
localization at the separator interface, as this can have a negative impact on
durability. To answer this question, we condense the structural and material
properties into the dimensionless quantity $hG$, i.e., the product between the
average grain size $G$ and the convective term $h$ (having units of 1/length).

The plot in Fig. 3 c illustrates the scaling of the peak current density at
the beginning of discharge $i_{max}$ as a function of the dimensionless
parameter $hG$. For Voronoi and square lattices (black and gray markers,
respectively), we found that the maximum relative current scales approximately
as $(hG)^{1/2}$. Only a small difference due to the specific cell geometry was
observed. The two microstructures small differences in tortuosity and wall
thickness. The square lattice has tortuosity 1, because the walls are straight
and orthogonal to the separator.

The wall thickness does not affect the peak current at the separator, but it
modifies the current decay away from the separator, as shown in Fig. 3 d
(green solid and dashed line). Doubling the wall thickness favors a more
homogeneous current density. The curves in Fig. 3 d refer to the conditions at
the beginning of discharge. Over time, the electrode walls and previously-
deposited lithium are both a pathway for Li-ion conduction, so the wall
thickness has strongest impact at the beginning of lithiation. Changing the
values of pore size $G$ and of the convective coefficient $h$ in such a way
that their product remains constant only affects the trend of the current away
from the separator, but not its peak value. The orange and black lines in Fig.
3 d indicate that decreasing the pore size is generally a good strategy to
diffuse plating over a larger area.

Fig. 3 c suggests that reducing the average cell size and increasing Li-ion
conductivity both contribute to lowering $i_{max}$. This result is intuitive,
however modeling allows for quantitative metrics in design and material
selection. This analysis show that the requirement for pore-size distribution
depends on the conductivity of the solid-electrolyte material. For highly
conductive solid-electrolytes like LLZO ($\kappa\sim 1$ mS/cm), if the average
pore diameter is $10\mu$ m [41] and interfacial impedance in the order or
$10\Omega$ cm2, the local maximum current density is about 5% of the global
current density. This significant reduction in current explains why mesoscale
architectures can be effective in improving the cycle life of solid-state
batteries with Li-metal, if the solid-electrolyte material is highly
conductive. The results in Fig. 2 and Fig. 3 indicate that the local
electrochemical conditions change significantly with the SE properties and it
is difficult to engineer the microstructure without an accurate model.

Figure 3: a, b) The current density distribution is plotted, as a fraction of
the applied value in the galvanostatic test. The horizontal axis indicates the
relative electrode thickness, ranging from zero (current collector) to one
(where the separator is located). The effect of bulk and surface conductivity
of the porous microstructure is highlighted by the difference in the two plots
(the values of $\kappa$ and $\rho$ are listed on each figure). Low bulk
conductivity tends to concentrate the current close to the separator interface
with lithium progressively filling a few pores at the time (figure b).
Conversely figure a shows a linear distribution of the current, with a larger
surface area engaged in electro-deposition and a lower peak current. c)
Scaling of the peak current with the dimensionless quantity $hG$ spanning over
several order of magnitudes. The convective parameters $h$ is the ratio
between surface and bulk conductivity and it appears in the boundary condition
Eq. 1b. $G$ is the average pore size. In addition to the Voronoi structure
shown in Fig. 2 a square lattice morphology was analyzed, and the results are
marked by black circles and gray squares respectively. In both cases the peak
current shows a clear scaling with respect to the dimensionless parameter as
$hG^{1/2}$. d) The images shows the sensitivity of the current distribution
away from the separator with respect to the wall thickness and the specific
values of $h$ and $G$, assuming that the product $h\cdot G$ is constant. The
curves refer to DOD=0, and the baseline case is marked by a solid line.
Comparing solid the dashed lines of the same color shows that doubling the
wall thickness and decreasing the pore size tend to make the current more
uniform and less concentrated at the separator interface.

##  4 Discussion

A mixed-conducting foam is investigated as a method to provide structure and
resilience to high-energy density solid-state batteries with a Li-metal anode.
The porous electrode structure greatly expands the surface area of the charge-
transfer interface, reducing the local current density to a fraction of the
applied planar current. The anode microstructure provides a conductive
scaffold that maintains accessibility of active materials during charge and
discharge With closed cells, the deposition front is confined within the
cells, preventing dendrites from puncturing the separator and shorting the
cell. If the deposition front becomes unstable dendrites tend to grow inward
rather than in the direction of the separator. Nevertheless, it is currently
unknown whether tension will be imposed onto the porous structure as the cell
become full of lithium. This may result in a failure mode that we have not
treated here. If a filled pore causes rupture of the matrix, it would not have
an effect on capacity, but it would increase tortuosity and potentially create
a pathway for dendrite growth in the direction of the separator.

Surface irregularities in the boundary layer can be included to study the
electrochemical conditions for their amplification, and predict the velocity
of dendrite propagation. The adaptive re-meshing features make this model well
equipped to study whether the growing interface is stable or not. These
effects are difficult to embed in a commercial finite element software, as
they require changes in the topology frequently during the simulation.

A stability analyses of the Li/Vapor interface was carried out, with some of
the results shown in Fig. 4. In the case of the left and center image a
sinusoidal boundary is assumed at the beginning of the simulation. Driven by
the local electric field, lithium tends to fill the grooves because the
current density is higher in the regions of higher curvature. The amplitude of
the instability decreases over time, and the interface turns into a circle
(marked in yellow in Fig. 4). The circle is off-centered with respect to the
square lattice because the current is higher towards the separator (marked by
the orange boundary in Fig. 4). Finally, a rough interface was generated by
random displacement of the points along a square boundary, as shown on the
right image of Fig. 4. Also in this case, the instability decays over time and
the Li-/Vapor interface becomes smooth.

These analyses indicated that Li-dendrite growth can be suppressed by
microstructural design. However, considerations of mechanical robustness
should be taken into account. Failure mechanisms originated from surface flaws
and from the pressure exerted on the surrounding walls as the pores are filled
with lithium remain to be investigated. This requires combining deposition
front tracking (presented here) with fracture mechanics, while capturing the
constitutive behavior of Li-metal and of scaffold material to accurately
predict stress. Analyses of the mechanical response will allow to optimize the
electrode microstructure for both conductivity and resilience.

Figure 4: The images show three examples of rough Li/Vapor interface evolving
into a smooth boundary. On the left and center images, the initial boundary
was chosen to be a sinusoidal wave of different wave length. On the right
image, the instability is created by random displacement of the nodes at the
Li/Vapor boundary. In all the three examples, the analyses indicate that the
Li/Vapor interface is stable with respect to non-uniform plating and surface
roughness.

##  5 Conclusions

The model suggests that a porous electrode matrix may provide means to reduce
shorting by confining Li plating into pores. During discharging, Li plating
proceeds inwards to fill the pore. The simulation suggest that the growth is
morphologically stable–it is possible that deleterious effects of dendrite
growth might be avoided by microstructural design of the electrode matrix.

The capacity of the battery is determined by the availability of pore volume.
However, the kinetics of charging are affected by the tortuosity of the matrix
and the surface area to volume ratio of the pores. Simulation provides scaling
approaches to determine what features of a porous electrode are the most
important, and how to tune the design to promote uniform current distribution
throughout the electrode thickness.

The primary contribution of this paper is the demonstration of an finite
element method that can be used and extended to other Li-battery structures
[42, 43, 44].

## Author Contributions

GB conceived the presented idea and took the lead in writing the manuscript.
GB and WCC developed modeling framework. TS contributed ideas on fabrication
requirements and economic aspects of the battery market. MB advised on the
modeling approach. All authors discussed the results and contributed to the
final manuscript.

## References

  * [1] Jing Chen, Jiawei Wu, Xiaodong Wang, An’an Zhou, and Zhenglong Yang.  Research progress and application prospect of solid-state electrolytes in commercial lithium-ion power batteries.  Energy Storage Materials, 35:70–87, 2021.
  * [2] Daxian Cao, Xiao Sun, Qiang Li, Avi Natan, Pengyang Xiang, and Hongli Zhu.  Lithium dendrite in all-solid-state batteries: Growth mechanisms, suppression strategies, and characterizations.  Matter, 3(1):57–94, 2020.
  * [3] Xinyou Ke, Yan Wang, Liming Dai, and Chris Yuan.  Cell failures of all-solid-state lithium metal batteries with inorganic solid electrolytes: Lithium dendrites.  Energy Storage Materials, 33:309–328, 2020.
  * [4] Eric Kazyak, Regina Garcia-Mendez, William S. LePage, Asma Sharafi, Andrew L. Davis, Adrian J. Sanchez, Kuan-Hung Chen, Catherine Haslam, Jeff Sakamoto, and Neil P. Dasgupta.  Li penetration in ceramic solid electrolytes: Operando microscopy analysis of morphology, propagation, and reversibility.  Matter, 2(4):1025–1048, 2020.
  * [5] Lukas Porz, Tushar Swamy, Brian W. Sheldon, Daniel Rettenwander, Till Frömling, Henry L. Thaman, Stefan Berendts, Reinhard Uecker, W. Craig Carter, and Yet-Ming Chiang.  Mechanism of lithium metal penetration through inorganic solid electrolytes.  Advanced Energy Materials, 7(20):1701003, 2017.
  * [6] Tushar Swamy, Richard Park, Brian W. Sheldon, Daniel Rettenwander, Lukas Porz, Stefan Berendts, Reinhard Uecker, W. Craig Carter, and Yet-Ming Chiang.  Lithium metal penetration induced by electrodeposition through solid electrolytes: Example in single-crystal li6la3zrtao12garnet.  Journal of The Electrochemical Society, 165(16):A3648–A3655, 2018\.
  * [7] Duo Zhang, Alvin Dai, Min Wu, Kang Shen, Teng Xiao, Guangya Hou, Jun Lu, and Yiping Tang.  Lithiophilic 3d porous cuzn current collector for stable lithium metal batteries.  ACS Energy Letters, 5(1):180–186, 2020.
  * [8] Qi Li, Shoupu Zhu, and Yingying Lu.  3d porous cu current collector/li-metal composite anode for stable lithium-metal batteries.  Advanced Functional Materials, 27(18):1606422, 2017.
  * [9] Shengxuan Lin, Yifan Wang, Yuhang Chen, Zihe Cai, Jiajia Xiao, Tahir Muhmood, and Xiaobin Hu.  Dendrite-free and ultra-long-life lithium metal anode enabled via a three-dimensional ordered porous nanostructure.  ACS Applied Materials & Interfaces, 13(35):41744–41752, 2021.
  * [10] Chengbin Jin, Ouwei Sheng, Jianmin Luo, Huadong Yuan, Cong Fang, Wenkui Zhang, Hui Huang, Yongping Gan, Yang Xia, Chu Liang, Jun Zhang, and Xinyong Tao.  3d lithium metal embedded within lithiophilic porous matrix for stable lithium metal batteries.  Nano Energy, 37:177–186, 2017.
  * [11] Shaomao Xu, Dennis W. McOwen, Chengwei Wang, Lei Zhang, Wei Luo, Chaoji Chen, Yiju Li, Yunhui Gong, Jiaqi Dai, Yudi Kuang, Chunpeng Yang, Tanner R. Hamann, Eric D. Wachsman, and Liangbing Hu.  Three-dimensional, solid-state mixed electron–ion conductive framework for lithium metal anode.  Nano Letters, 18(6):3926–3933, 2018.
  * [12] Dennis W. McOwen, Shaomao Xu, Yunhui Gong, Yang Wen, Griffin L. Godbey, Jack E. Gritton, Tanner R. Hamann, Jiaqi Dai, Gregory T. Hitz, Liangbing Hu, and Eric D. Wachsman.  3d-printing electrolytes for solid-state batteries.  Advanced Materials, 30(18):1707132, 2018.
  * [13] Jonghyeok Yun, Beom-Kyeong Park, Eun-Seo Won, Seung Hyun Choi, Hyon Chol Kang, Jung Ho Kim, Min-Sik Park, and Jong-Won Lee.  Bottom-up lithium growth triggered by interfacial activity gradient on porous framework for lithium-metal anode.  ACS Energy Letters, 5(10):3108–3114, 2020.
  * [14] Hailong Qiu, Tianyu Tang, Muhammad Asif, Xiaoxiao Huang, and Yanglong Hou.  3d porous cu current collectors derived by hydrogen bubble dynamic template for enhanced li metal anode performance.  Advanced Functional Materials, 29(19):1808468, 2019.
  * [15] Yuhwan Hyeon, Jaewoo Lee, Hamzeh Qutaish, Sang A Han, Seung Hyun Choi, Sung Won Moon, Min-Sik Park, Dongmok Whang, and Jung Ho Kim.  Lithium metal storage in zeolitic imidazolate framework derived nanoarchitectures.  Energy Storage Materials, 33:95–107, 2020.
  * [16] Shu-Hua Wang, Ya-Xia Yin, Tong-Tong Zuo, Wei Dong, Jin-Yi Li, Ji-Lei Shi, Chang-Huan Zhang, Nian-Wu Li, Cong-Ju Li, and Yu-Guo Guo.  Stable li metal anodes via regulating lithium plating/stripping in vertically aligned microchannels.  Advanced Materials, 29(40):1703729, 2017.
  * [17] Chengzhi Wang, Aoxuan Wang, Lingxiao Ren, Xuze Guan, Donghong Wang, Anping Dong, Chanyuan Zhang, Guojie Li, and Jiayan Luo.  Controlling li ion flux through materials innovation for dendrite-free lithium metal anodes.  Advanced Functional Materials, 29(49):1905940, 2019.
  * [18] Eongyu Yi, Hao Shen, Stephen Heywood, Judith Alvarado, Dilworth Y. Parkinson, Guoying Chen, Stephen W. Sofie, and Marca M. Doeff.  All-solid-state batteries using rationally designed garnet electrolyte frameworks.  ACS Applied Energy Materials, 3(1):170–175, 2020.
  * [19] Yao Zhao, Shuaiguo Hao, Li Su, Zhipeng Ma, and Guangjie Shao.  Hierarchical cu fibers induced li uniform nucleation for dendrite-free lithium metal anode.  Chemical Engineering Journal, 392:123691, 2020.
  * [20] Wanfeng Yang, Jiawei Wang, Wensheng Ma, Chaoqun Dong, Guanhua Cheng, and Zhonghua Zhang.  Free-standing cuo nanoflake arrays coated cu foam for advanced lithium ion battery anodes.  Journal of Power Sources, 333:88–98, 2016.
  * [21] Gang Zhang, Bin Deng, Quan yao Liu, Hua Yang, Mao xiang Jing, Shahid Hussain, and Tahani Saad AlGarni.  Enhanced areal specific capacity and energy density of solid-state lithium battery by using porous aluminum foam.  Journal of Energy Storage, 33:102167, 2021.
  * [22] R Edwin Garcıa, Yet-Ming Chiang, W Craig Carter, Pimpa Limthongkul, and Catherine M Bishop.  Microstructural modeling and design of rechargeable lithium-ion batteries.  Journal of The Electrochemical Society, 152(1):A255, 2004.
  * [23] Giovanna Bucci, Tushar Swamy, Yet-Ming Chiang, and W Craig Carter.  Modeling of internal mechanical failure of all-solid-state batteries during electrochemical cycling, and implications for battery design.  Journal of Materials Chemistry A, 5(36):19422–19430, 2017.
  * [24] Davide Grazioli, Marco Magri, and Alberto Salvadori.  Computational modeling of li-ion batteries.  Computational Mechanics, 58(6):889–909, 2016.
  * [25] Giovanna Bucci, Tushar Swamy, Sean Bishop, Brian W Sheldon, Yet-Ming Chiang, and W Craig Carter.  The effect of stress on battery-electrode capacity.  Journal of The Electrochemical Society, 164(4):A645, 2017.
  * [26] Giovanna Bucci, Brandon Talamini, Ananya Renuka Balakrishna, Yet-Ming Chiang, and W Craig Carter.  Mechanical instability of electrode-electrolyte interfaces in solid-state batteries.  Physical Review Materials, 2(10):105407, 2018.
  * [27] Anja Bielefeld, Dominik A Weber, and Jürgen Janek.  Microstructural modeling of composite cathodes for all-solid-state batteries.  The journal of physical chemistry C, 123(3):1626–1634, 2018.
  * [28] Benedikt Prifling, Daniel Westhoff, Denny Schmidt, Henning Markoetter, Ingo Manke, Volker Knoblauch, and Volker Schmidt.  Parametric microstructure modeling of compressed cathode materials for li-ion batteries.  Computational Materials Science, 169:109083, 2019.
  * [29] Martin Finsterbusch, Timo Danner, Chih-Long Tsai, Sven Uhlenbruck, Arnulf Latz, and Olivier Guillon.  High capacity garnet-based all-solid-state lithium batteries: fabrication and 3d-microstructure resolved modeling.  ACS applied materials & interfaces, 10(26):22329–22339, 2018.
  * [30] David E Stephenson, Bryce C Walker, Cole B Skelton, Edward P Gorzkowski, David J Rowenhorst, and Dean R Wheeler.  Modeling 3d microstructure and ion transport in porous li-ion battery electrodes.  Journal of The Electrochemical Society, 158(7):A781, 2011.
  * [31] Adam A Wheeler, Bruce T Murray, and Robert J Schaefer.  Computation of dendrites using a phase field model.  Physica D: Nonlinear Phenomena, 66(1-2):243–262, 1993.
  * [32] Alain Karma and Wouter-Jan Rappel.  Phase-field model of dendritic sidebranching with thermal noise.  Physical review E, 60(4):3614, 1999.
  * [33] Toshio Suzuki, Machiko Ode, Seong Gyoon Kim, and Won Tae Kim.  Phase-field model of dendritic growth.  Journal of Crystal Growth, 237:125–131, 2002.
  * [34] Ingo Steinbach.  Phase-field models in materials science.  Modelling and simulation in materials science and engineering, 17(7):073001, 2009.
  * [35] Zijian Hong and Venkatasubramanian Viswanathan.  Phase-field simulations of lithium dendrite growth with open-source software.  ACS Energy Letters, 3(7):1737–1743, 2018.
  * [36] Jiawei Zhang, Yuanpeng Liu, Changguo Wang, and Huifeng Tan.  An electrochemical-mechanical phase field model for lithium dendrite.  Journal of The Electrochemical Society, 168(9):090522, 2021.
  * [37] Qiao Wang, Geng Zhang, Yajie Li, Zijian Hong, Da Wang, and Siqi Shi.  Application of phase-field method in rechargeable batteries.  npj Computational Materials, 6(1):1–8, 2020.
  * [38] Liting Gao and Zhansheng Guo.  Phase-field simulation of li dendrites with multiple parameters influence.  Computational Materials Science, 183:109919, 2020.
  * [39] Lei Chen, Hao Wei Zhang, Lin Yun Liang, Zhe Liu, Yue Qi, Peng Lu, James Chen, and Long-Qing Chen.  Modulation of dendritic patterns during electrodeposition: A nonlinear phase-field model.  Journal of Power Sources, 300:376–385, 2015.
  * [40] Chih-Hung Chen and Chun-Wei Pao.  Phase-field study of dendritic morphology in lithium metal batteries.  Journal of Power Sources, 484:229203, 2021.
  * [41] Hao Shen, Eongyu Yi, Stephen Heywood, Dilworth Y. Parkinson, Guoying Chen, Nobumichi Tamura, Stephen Sofie, Kai Chen, and Marca M. Doeff.  Scalable freeze-tape-casting fabrication and pore structure analysis of 3d llzo solid-state electrolytes.  ACS Applied Materials & Interfaces, 12(3):3494–3501, 2020.
  * [42] Benjamin Delattre, Ruhul Amin, Jonathan Sander, Joël De Coninck, Antoni P. Tomsia, and Yet-Ming Chiang.  Impact of pore tortuosity on electrode kinetics in lithium battery electrodes: Study in directionally freeze-cast LiNi0.8co0.15al0.05o2(NCA).  Journal of The Electrochemical Society, 165(2):A388–A395, 2018\.
  * [43] Ruhul Amin, Benjamin Delattre, Antoni P. Tomsia, and Yet-Ming Chiang.  Electrochemical characterization of high energy density graphite electrodes made by freeze-casting.  ACS Applied Energy Materials, 1(9):4976–4981, 2018.
  * [44] Linsen Li, Randall M. Erb, Jiajun Wang, Jun Wang, and Yet-Ming Chiang.  Fabrication of low-tortuosity ultrahigh-area-capacity battery electrodes through magnetic alignment of emulsion-based slurries.  Advanced Energy Materials, 9(2):1802472, 2019.