#  Stability and Uncertainty Propagation in Power Networks:
A Lyapunov-based Approach with applications to
Renewable Resources Allocation

Mohamad H. Kazma, Graduate Student Member, IEEE and Ahmad F.
Tah$\text{a}^{\diamond}$, Member, IEEE ⋄Corresponding author. This work is
supported by National Science Foundation under Grants 2152450 and 2151571. The
authors are with the Civil $\&$ Environmental Engineering and Electrical $\&$
Computer Engineering Departments, Vanderbilt University, 2201 West End Ave,
Nashville, Tennessee 37235. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The rapid increase in the integration of intermittent and stochastic renewable
energy resources (RER) introduces challenging issues related to power system
stability. Interestingly, identifying grid nodes that can best support
stochastic loads from RER, has gained recent interest. Methods based on
Lyapunov stability are commonly exploited to assess the stability of power
networks. These strategies approach quantifying system stability while
considering: (i) simplified reduced order power system models that do not
model power flow constraints, or (ii) data-driven methods that are prone to
measurement noise and hence can inaccurately depict stochastic loads as system
instability. In this paper, while considering a nonlinear differential
algebraic equation (NL-DAE) model, we introduce a new method for assessing the
impact of uncertain renewable power injections on the stability of power
system nodes/buses. The identification of stable nodes informs the
operator/utility on how renewables injections affect the stability of the
grid. The proposed method is based on optimizing metrics equivalent to the
Lyapunov spectrum of exponents; its underlying properties result in a
computationally efficient and scalable stable node identification algorithm
for renewable energy resources allocation. The proposed method is validated on
the IEEE 9-bus and 200-bus networks.

###### Index Terms:

Power system modeling, nonlinear differential algebraic models, Lyapunov
exponents, nonlinear stability

##  I Introduction and paper contributions

Power systems are becoming more reliant on renewable energy resources (RER)
for power supply, in a much needed effort to decarbonize our electricity
systems. From a stability perspective, the uncertainty and intermittent
generation of RER may impede power system operations. This is a result of the
range of intermittent perturbations acting on the synchronous states [1]. As
such, studying the transient stability and resilience of power systems with
increased amounts of integrated RER has become essential. A fundamental
approach for such stability analysis is typically provided by simulating the
impact of RER integration through load disturbances and inertia changes, and
then assessing the system settling time.

The requirement for performing transient stability analysis entails the
identification of potential disturbances that might lead to system
instabilities [2]. To assess system stability, Lyapunov stability methods are
widely adopted [3]. The methods are divided into two categories: (i) a direct
stability method that is based on quantifying a Lyapunov energy function [4]
and (ii) a Lyapunov exponents (LEs) method that is based on characterizing
infinitesimal separation rates of system trajectories [2]. The former method
is based on quantifying an energy function that is indicative of system
stability. The latter method is adapted from the field of chaos and was first
introduced by [5] in the context of power networks for studying transient
chaotic swings. The aforementioned direct stability method is well-developed
however it becomes complex to quantify when considering nonlinear dynamical
systems. On the other hand, quantifying LEs of a nonlinear system arises more
naturally [6].

For power system stability, the maximal Lyapunov exponent (MLE) is able to
depict general characteristics of transient stability [7]. A negative MLE sign
is indicative of a stability, whereas a positive MLE informs that the system
is unstable. Theoretically, LEs measure the average exponential rate of
convergence and divergence of a pair of nearby trajectories in a
multidimensional state-space. The exponents are key to the determination of
stability since the exponents are invariant to any initial condition
perturbation that are within the same stability region [8]. With that in mind,
the main objective of this paper is to provide a framework for quantifying
stability based on metrics related to the LEs of the system; the exponents are
computed after an uncertain RER power injection is applied to a bus. The
method quantifies power system stability and uncertainty propagation from
negative loads that model injections of uncertain renewables. As such, the
proposed method informs the operator about the impact of allocating renewable
injections, at network buses, on the overall power grid stability.

Literature Review. LE-based approach for power system stability and its
applications to quantifying network characteristics has gained recent
interest. The literature can be divided into (i) simply quantifying the
stability of a power network, and (ii) identifying network properties. From a
stability perspective, a power system’s stability is usually studied in terms
of voltages, angles, and frequencies [9]. A method based on LEs [10], is
utilized to study the swings that lead to instability by monitoring the rotor
angle of generator machines. In [11], the transient stability of power systems
post-fault is studied by analyzing the system’s LEs. Voltage stability is
considered from phasor measurement units (PMU) data in [12]. The study
considers a model-free approach that monitors voltage data and then estimates
the systems LEs for stability predictions. In [2], global system stability is
studied by computing LEs and their corresponding Lyapunov vectors. The system
machines are modeled according to the swing equations and then the proposed
approach was demonstrated on a $9$-bus network. A model-free estimation
framework is considered for assessing rotor angle stability in [13]. The
proposed method relies on least-squares state estimation for computing rotor
angle estimates from PMU measurements. Similarly, in [7] the rotor angle
stability from MLE estimation is considered while utilizing a recursive least
squares estimation framework that is less time-consuming.

The LE-based stability methods outlined above consider the power system
dynamics by (i) modeling simplified and reduced-order power machine models or
(ii) utilizing model-free approaches that rely on PMU measurement data.
Moreover, these stability methods consider only voltage stability and rely on
dynamic state estimation methods for estimating rotor angle stability. The
main drawback of such methods are as follows. Simplified machine models that
do not model the both differential equations and algebraic constraints are not
equipped to model RER loads along with their uncertainty. The reason is that
renewables injections do not explicitly appear in the differential equations
but appear in the algebraic power flow constraints. On the other hand, a
direct model-free approach is only appropriate for disturbance and noise free
systems [14]. The reason is that any disturbance or uncertainty would be
wrongly considered as chaotic behavior when computing LEs. As such, there is
no clear way to explicitly delineate how renewable injections would impact
grid dynamics and the different types of stability. Moreover, when relying on
estimation frameworks for rotor angle stability, the estimation error can also
be incorrectly assigned as chaos in the system. Accordingly, in this paper we
compute a power system’s LEs by considering a model-based approach that models
the complete nonlinear differential algebraic power system model and therefore
depicts both differential and algebraic states.

The second branch of research that goes beyond simply studying the stability
of power systems and has recently gained interest, is the application of LEs
methods that enable the identification of certain network properties. The
studies [15, 16] both utilize a LEs algorithm to identify coherent generator
groups that generate synchronous swings. The effect of chaos from wind power
is predicted from a model decomposition approach that is based on LEs
computations in [17]. An optimal PMU allocation approach for full network
observability and transient stability improved is proposed in [18]. The
approach relies on LEs to identify the critical buses that provide the most
useful information for each system node. Nonetheless, influential nodes in
power systems play a stable role in the structure and dynamics of the network
[19]. The allocation of RER is proposed in [20]. The proposed method offers a
heuristic placement algorithm that estimates LEs based on voltage measurements
and then ranks the candidate buses based on the values of the LEs computed
from voltage stability. However, the framework proposed in [20] offers a
candidate node rank that (i) is based on model-free LEs estimation, where
instability can result from measurement noise, (ii) computes LEs while
considering chaos resulting from voltage stability only, and (iii) quantifies
node stability by measuring the LEs at that node while not taking into account
uncertainty propagation onto the other nodes.

Paper Contributions. Motivated by the aforementioned limitations in the recent
literature, the main contributions are as follows

  * •

We introduce a framework that enables the node stability assessment of a NL-
DAE representation of a power network. In specific, we provide quantitative
stability measures for identifying stable and unstable nodes in a power
network from a Lyapunov stability perspective. The proposed stability measures
are based on the Lyapunov characteristic spectrum of exponents, which are
computed for a system under renewable disturbances. These stability metrics
allow for (i) computing the overall stability by considering frequency,
voltage and rotor angle stability and (ii) quantifying renewables uncertainty
propagation, that is applied to a single node, on the overall power network
stability.

  * •

Based on the introduced stability framework, we propose an algorithm for
identifying optimal stable nodes that can accommodate renewables injections.
This approach allows to inform the operator where to practically allocate
renewables while maintaining the overall stability of the power network that
is inundated with uncertain fluctuations from RER injections. The result is an
ordered set of nodes that reflect stable nodes that minimally impact overall
power system stability.

  * •

The validity of the approach is demonstrated on standard power systems. The
robustness and effectiveness of the proposed optimal allocation problem is
shown to be robust against uncertain RERs.

We note that in this paper we do not review methods for RER allocation or
sizing. Such methods can consider different operational objectives (efficiency
maximization, and power loss and cost minimization). Instead, we are concerned
with quantifying the effect of RER allocation on the overall stability of a
power network and thus provide a method to inform operators regarding the
practical allocation of RER from a stability perspective. Readers are referred
to [21] for a more comprehensive review on RER allocation and sizing.

Notation. Let $\mathbb{N}$, $\mathbb{R}$, $\mathbb{R}^{n}$, and
$\mathbb{R}^{p\times q}$ denote the set of natural numbers, real numbers,
real-valued row vectors with size of $n$, and $p$-by-$q$ real matrices
respectively. The symbol $\otimes$ denotes the Kronecker product. The
cardinality of a set $\mathcal{N}$ is denoted by $|\mathcal{N}|$. The
operators $\mathrm{det}(\bm{A})$ returns the determinant of matrix $\bm{A}$,
$\mathrm{trace}(\bm{A})$ returns the trace of matrix $\bm{A}$. The operator
$\\{\bm{x}_{i}\\}_{i=0}^{\mathrm{N}}\in\mathbb{R}^{\mathrm{N}n}$ constructs a
column vector that concatenates vectors $\bm{x}_{i}\in\mathbb{R}^{n}$ for all
$i\in\\{0,1,\ldots,\mathrm{N}\\}$.

Paper Organization. The paper is organized as follows: Section II provides the
power system model. Section III provides preliminaries on LEs computation.
Section IV introduces the proposed stable node identification framework and
the optimal RER allocation problem for NL-DAEs. The numerical results are
presented in Section V, and Section VI concludes this paper.

##  II Problem formulation

In this section, we present the underlying NL-DAE power system model and
represent the full model as a nonlinear ODE. The resulting nonlinear ODE
system, that includes the algebraic constraints, is discretized using the
trapezoidal implicit discretization method.

###  II-A NL-DAE power system state-space model

We consider a NL-DAE formulation of a power system. This system depicts the
standard two axes $4$-th order transient model of a synchronous generator [22,
Ch. 7]. The considered model excludes exciter dynamics and turbine governor,
meaning that each of the machines has four states and two control inputs. We
note that the controller (i.e., governor response) adjusts the active power
produced by a generator to stabilize the system after a disturbance from
renewables.

To this end, we borrow the notation from [23]. A power system
$(\mathcal{N},\mathcal{E})$ can be represented graphically, where
$\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of transmission
lines, $\mathcal{N}=\mathcal{G}\cup\mathcal{L}\cup\mathcal{R}$ is the set of
all buses in the network. The number of buses within the network is
$N:=|\mathcal{N}|$, while the number of generator, load and renewable buses
are $G:=|\mathcal{G}|$, $L:=|\mathcal{L}|$ and $R:=|\mathcal{R}|$. The
physics-based components of the electromechanical transients can be written in
a semi-implicit NDAE form as follows

$\displaystyle\textit{generator dynamics}:\;\;\dot{\bm{x}}_{d}$
$\displaystyle=\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u}),$ (1a)
$\displaystyle\textit{algebraic constraints}:\;\;\bm{0}$
$\displaystyle=\bm{g}(\bm{x}_{d},\bm{x}_{a}),$ (1b)

where the dynamic states of the synchronous machine can be defined as
$\bm{x}_{d}:=\bm{x}_{d}(t)\in\mathbb{R}^{4G}$, the algebraic states can be
defined as $\bm{x}_{a}:=\bm{x}_{a}(t)\in\mathbb{R}^{2G+2N}$ and the input of
the system can be defined as $\bm{u}:=\bm{u}(t)\in\mathbb{R}^{2G}$. Let
$n_{d}:=4G$ and $n_{a}:=2G+2N$ define the size of the differential and
algebraic states. Let $n_{u}:=2G$ define the size of the input state vector
$\bm{u}$. Nonlinear mapping function $\bm{f}(\cdot)$ depicts the synchronous
machine dynamics, where
$\bm{f}(\cdot):\mathbb{R}^{4G}\times\mathbb{R}^{2G}\times\mathbb{R}^{2G}\rightarrow\mathbb{R}^{4G}$.
Nonlinear mapping function $\bm{g}(\cdot)$ depicts the algebraic constraints
and the power balance equations, where
$\bm{g}(\cdot):\mathbb{R}^{4G}\times\mathbb{R}^{2G}\times\mathbb{R}^{2N}\rightarrow\mathbb{R}^{2G+2N}$.

The physics-based components of the electromechanical transients depicting the
$4$-th order transient model are summarized in Appendix A. Accordingly, the
differential state, algebraic state and input vectors are summarized as
follows

$\displaystyle\bm{x}_{d}$
$\displaystyle=\left\\{\\{\delta_{i}\\}_{i=0}^{G},\;\\{\omega_{i}\\}_{i=0}^{G},\;\\{E^{{}^{\prime}}_{i}\\}_{i=0}^{G},\;\\{{T}_{\mathrm{N}i}\\}_{i=0}^{G}\right\\}^{\top},\leavevmode\nobreak\
$ (2a) $\displaystyle\bm{x}_{a}$
$\displaystyle=\left\\{\\{{P}_{\mathrm{G}i}\\}_{i=0}^{G},\;\\{{Q}_{\mathrm{G}i}\\}_{i=0}^{G},\;\\{{v}_{i}\\}_{i=0}^{N},\;\\{{\theta}_{i}\\}_{i=0}^{N}\right\\}^{\top},\leavevmode\nobreak\
$ (2b) $\displaystyle\bm{u}$
$\displaystyle=\left[\\{{E}_{\mathrm{fd}i}\\}_{i=0}^{G},\;\\{{T}_{\mathrm{r}i}\\}_{i=0}^{G}\right]^{\top}.\leavevmode\nobreak\
$ (2c)

where the time varying components in (27) are: $\delta_{i}$ the rotor angle
$\mathrm{(rad)}$, $\omega_{i}$ generator rotor speed $\mathrm{(rad/sec)}$,
${E}^{{}^{\prime}}_{i}$ generator transient voltage $\mathrm{(pu)}$,
${T}_{\mathrm{N}i}$ generator mechanical torque $(\mathrm{pu})$. Generator
inputs are: $E_{\mathrm{fd}i}$ generator internal field voltage
$(\mathrm{pu})$, $T_{\mathrm{r}i}$ governor reference signal $(\mathrm{pu})$.
Time varying components in (29) are $\theta_{ij}=\theta_{i}-\theta_{j}$
$\mathrm{(pu)}$ is the bus angle and $v_{i}$ is the bus voltage
$\mathrm{(pu)}$. Real and reactive power for a generator are denoted as
$P_{\mathrm{G}i}$ and $Q_{\mathrm{G}i}$. It is noteworthy to mentioned that
the load injections that model RER are included in the power balance equations
(29) as real and reactive power $P_{\mathrm{R}i}$ and $Q_{\mathrm{R}i}$.

###  II-B NL-DAEs: existence and uniqueness of a solution

The uniqueness and existence of a solution for NL-DAEs, as compared to
nonlinear ODEs (NL-ODEs), for any initial condition $(\bm{x}_{0},\bm{u}_{0})$
is not always guaranteed [24]. The existence and uniqueness of solutions to
NL-DAEs holds true if and only if the system (1) is strangeness-free [25,
Hypothesis 4.2], whereby the strangeness index is equal to zero. For brevity,
refer to [25, 26] for the detailed hypothesis that defines the strangeness
index of NL-DAEs. We note here that this index is a generalization of the
differentiation index of DAEs, that is defined as follows.

######  Definition 1.

The differentiation index [27, 25, 28] more ref of nonlinear or linear DAE
systems is defined as the number of times the algebraic equations are
differentiated to obtain a set of ODEs.

A linear DAE can be written as $\bm{E}\dot{\bm{x}}=\bm{A}\bm{x}+\bm{B}\bm{u}$.
This representation can be formulated by linearizing the system (1) around an
operating point where the state vector $\bm{x}$ is defined as
$\bm{x}:=\left[\bm{x}_{d},\;\bm{x}_{a}\right]^{\top}\in\mathbb{R}^{n_{d}+n_{a}}$.
The singular mass matrix $\bm{E}\in\mathbb{R}^{n_{d}+n_{a}}$ has ones on its
diagonal for the differential equations and zeros for the algebraic equations.
The time-invariant state-space matrices are defined as
$\bm{A}:=\left[\bm{A}_{d},\;\bm{A}_{a}\right]^{\top}\in\mathbb{R}^{n_{d}+n_{a}}$
and $\bm{B}:=\left[\bm{B}_{d}\right]\in\mathbb{R}^{n_{u}}$. Here
$\bm{A}_{d}\in\mathbb{R}^{n_{d}}$ represents the linearized differential
equations state-space matrix and $\bm{A}_{a}\in\mathbb{R}^{n_{a}}$ the
linearized algebraic equations state-space matrix.

We note that a linear DAE systems is of index if and only if the system is
regular. Regularity is considered an essential property for DAEs. This
property ensures the existence of consistent unique solutions for every
initial condition $(\bm{x}_{0},\bm{u}_{0})$ [29].

######  Definition 2.

A linearized DAE around an initial state $\bm{x}_{0}$ is regular if and only
if there exists $s\in\mathbb{C}$ such that $\mathrm{det}(s\bm{E}-\bm{A})\neq
0.$ Regularity is therefore characterized by the matrix pair
$(\bm{E},\bm{A})$.

A DAE power network model is regular if and only if there exists a path
whereby each load follows this path to a generator bus [29]. This holds true
for the IEEE case networks considered in the case studies section in this
work. The regularity and differentiation index of linearized power system (1)
are assessed in [23]. It is shown that the power system has a differentiation
index of one and is regular; see [23] for additional information. Such
conditions guarantee that for each consistent initial condition, a unique
solution always exists. With the evidence that the linearized DAE system, we
can now infer that the NL-DAE is strangeness-free without the rigorous proof.
Hence, for certain rank conditions every regular and linear DAE with
sufficiently smooth $(\bm{E},\bm{A})$ matrix pair satisfies the hypothesis
regarding the vanishing strangeness index, i.e., the system is strangeness-
free; see [25].

With that in mind, in order to ensure that there exists a solution to the NL-
DAE power system (1), the following assumption holds true throughout this
paper.

######  Assumption 1.

The NL-DAE (1) is considered strangeness-free, regular and of differentiation
index one. As such, a unique solution exists for any initial conditions
$(\bm{x}_{0},\bm{u}_{0})$.

The aforementioned assumption implies that the partial derivatives of the
differential and algebraic state-space functions $\bm{f}(\cdot)$ and
$\bm{g}(\cdot)$ with respect to $\bm{x}_{d}$ and $\bm{x}_{a}$ are non-singular
[25]. We note here that this assumption is mild and holds true for the power
system cases considered in the numerical studies section. In practical terms,
the numerical solvability of the DAE system (1) directly implies the existence
and uniqueness of solution for any input and consistent initial conditions
[30]. The time-domain simulations presented in Section V provide explicit
proof and validate Assumption 1 for the NL-DAE (1).

###  II-C NL-DAE to NL-ODE system transformation

NL-DAEs are considered stiff dynamical systems with time constants that can
span several orders of magnitude [31]. The algebraic constraints are
considered to exhibit null time constants. Stability and control of such NL-
DAE systems is limited, while that of NL-ODEs is not. Typically, an ODE system
is considered when modeling power system dynamics by either (i) neglecting the
algebraic constraints or (ii) relying on a decoupled model [23]. This limits
the ability to assess the transient stability of power system, in particular,
power networks that have RER. With that in mind, we consider a structure
preserving NL-DAE to NL-ODE transformation of the power system model (1). This
transformation allows for the complete modeling of the differential and
algebraic equations (27)–(29) while being modeled as a set of ODEs.

Transformation relies on applying the Inverse function theorem (IFT) [32,
Theorem 3.3.1] to resolve the algebraic constraints into ODEs. This method
requires the differentiation of the algebraic constraints in (1b) with respect
to time variable $t$. With that in mind, the NL-DAE system (1) can be
rewritten as

$\displaystyle\dot{\bm{x}}_{d}$
$\displaystyle=\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u}),$ (3a)
$\displaystyle\dot{\bm{x}}_{a}$
$\displaystyle=\tilde{\bm{g}}(\bm{x}_{d},\bm{x}_{a},\bm{u})=-(\bm{G}_{\bm{x}_{a}})^{-1}\bm{G}_{\bm{x}_{d}}\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u}),\leavevmode\nobreak\
$ (3b)

where matrices defined as
$\bm{G}_{\bm{x}_{a}}:=\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{a}}$
$\in\mathbb{R}^{n_{a}\times n_{a}}$ and
$\bm{G}_{\bm{x}_{d}}:=\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{d}}$
$\in\mathbb{R}^{n_{a}\times n_{d}}$ are the Jacobian matrices of the algebraic
constraints with respect to $\bm{x}_{a}$ and $\bm{x}_{d}$, respectively.

The NL-ODE system (3) obtained using the IFT method depicts the full dynamic
and algebraic relationships that are inherent to the system while also being
modeled as a set of ODEs. We note that the resulting NL-ODE system
representation of the algebraic constraints (1b) depends on (1a); it therefore
depends on control input $\bm{u}$. Readers are referred to [33] for a more
though discussion regarding the validity and accuracy of the transformed NL-
DAE model.

###  II-D Discrete-time modeling of NDAE power systems

Nonlinear systems, in particular DAE power system models, exhibit stiff
dynamics. Networks systems with stiff dynamics can be characterized by time
constants on local nodes or subsystems that have a significant contrasting
magnitude. As such, the discretization method of choice relies on the stiff
dynamics and the desired accuracy of the discretization. In practice, in order
to obtain stable and computational efficient dynamics, power systems that
exhibit transient conditions are solved numerically using implicit
discretization methods [30]. Explicit methods such as the implicit Runge-Kutta
(IRK) method cannot deal properly with stiff dynamics [34].

Typically, in the context of power system discrete-time dynamic modeling, the
implicit discretization methods utilized are as follows: (i) backward Euler
(BE) method [34], (ii) backward differential formulas (BDF) known as Gear’s
method [35], and (iii) trapezoidal implicit (TI) method [36, 34]. In [33], BDF
and TI methods are investigated for the dynamical system in (1). The results
show an accurate depiction of the transient dynamical states. With that in
mind, for the purpose of characterizing the LEs of the NL-DAE system (3), we
refer to the use of the TI discretization method. The main advantage of
implicit TI method is that it can handle a large class of stiff nonlinear
dynamical systems with large discretization time steps. To such end, the
nonlinear power system dynamics (3) can be rewritten in discrete-time as

$\displaystyle\bm{x}_{d,k}$
$\displaystyle=\bm{x}_{d,k-1}+\tilde{h}\left(\bm{f}(\bm{z}_{k})+\bm{f}(\bm{z}_{k-1})\right),$
(4a) $\displaystyle\bm{x}_{a,k}$
$\displaystyle=\bm{x}_{a,k-1}+\tilde{h}\left(\tilde{\bm{g}}(\bm{z}_{k})+\tilde{\bm{g}}(\bm{z}_{k-1})\right),$
(4b)

where vector
$\bm{z}_{k}:=[\bm{x}_{d,k},\;\bm{x}_{a,k},\;\bm{u}_{k}]^{\top}\in\mathbb{R}^{n_{d}+n_{a}+n_{u}}$
and
$\bm{x}_{k}:=[\bm{x}_{d,k},\;\bm{x}_{a,k}]^{\top}\in\mathbb{R}^{n_{d}+n_{a}}$
for time step $k$. The discretization time step size $\tilde{h}$ is defined as
$\tilde{h}:=0.5h$, where $h$ is the simulation time step size.

The Newton-Raphson (NR) algorithm is implemented to solve a system of implicit
discrete-time NL-DAEs; see [37]. For time-domain simulations, the algorithm
computes the unknown future states $\bm{x}_{k}$ using an iterative framework
until solution convergence, and then it advances to the next time step. Refer
to Appendix B for a detailed description of the NR method.

Moving forward we shall imply the dependence on control input, that is,
$\bm{f}(\bm{z}_{k}):=\bm{f}(\bm{x}_{k})$ and
$\tilde{\bm{g}}(\bm{z}_{k}):=\tilde{\bm{g}}(\bm{x}_{k})$. The discrete-time
nonlinear power system (4) can be written succinctly as

$\bm{x}_{k}=\bm{x}_{k-1}+\left[\tilde{\bm{f}}(\bm{x}_{k})+\tilde{\bm{f}}(\bm{x}_{k-1})\right],$
(5)

where nonlinear function
$\tilde{\bm{f}}(\bm{x}_{k}):=\left[\bm{f}(\bm{x}_{k}),\;\tilde{\bm{g}}(\bm{x}_{k})\right]^{\top}\in\mathbb{R}^{n}$,
such that $n:=n_{d}+n_{a}$ represents the differential and algebraic states.

##  III A Lyapunov-based stability approach

This section provides a brief background on the stability of dynamical systems
from a chaos and ergodic theory perspective; it describes the rationale and
practical implementation of LEs for stability of nonlinear systems.

###  III-A Lyapunov stability: maximal Lyapunov exponent

A well-known method for assessing the stability of nonlinear system trajectory
stems for Lyapunov’s indirect method on stability [38]. This method, common to
the fields of chaos and ergodicity, involves computing the Lyapunov spectrum
of exponents. These exponents represent the exponential convergence or
divergence of nearby perturbed trajectories of an attractor in the state-space
[13]. LEs denoted by $\bm{\lambda}$ are an indication of the asymptotic
behavior dynamical systems [39]. As such, this method presents a tool to
assess the systems sensitive dependence of initial conditions, i.e., analyze
the system’s response to transients induced due to perturbed initial
conditions. To calculate the LEs of nonlinear power system models, the use of
the variational equation of the original dynamical system becomes necessary
[40].

With that in mind, consider two nearby trajectories $\bm{x}_{k}$ and
$\bm{x}_{k}+\bm{\delta}{\bm{x}_{k}}$ resulting from initial states
$\bm{x}_{0}$ and $\bm{x}_{0}+\bm{\delta}{\bm{x}_{0}}$. Note that
$\bm{\delta}{\bm{x}_{0}}\in\mathbb{R}^{n}$ is an infinitesimal perturbation
$\varepsilon>0$ to initial conditions $\bm{x}_{0}$ and its exponential decay
or growth for $k\in\\{0\;,1\;,\;\ldots\;,\mathrm{N}-1\\}$ is denoted as
$\bm{\delta}\bm{x}_{k}\in\mathbb{R}^{n}$. The initial conditions are the
operating conditions obtained be solving the power flow equations.
Consequently, the uncertainty from renewables are modeled by perturbing the
initial loads in power flow equations. Note that $\mathrm{N}$ is the discrete-
time interval from the simulation, that is, $\mathrm{N}:=t/h$. To that end,
the variational system of the discrete-time nonlinear system (5) can be
written as

$\leavevmode\nobreak\
\bm{\delta}\bm{x}_{k}=\bm{\Phi}_{0}^{k}(\bm{x}_{0})\bm{\delta}\bm{x}_{0},$ (6)

where
$\bm{\Phi}_{0}^{k}(\bm{x}_{0}):=\left(\bm{I}_{n}+\tfrac{\partial\tilde{\bm{f}}(\bm{x}_{k},\bm{x}_{k-1})}{\partial\bm{x}_{k}}\right)\tfrac{\partial\bm{x}_{k}}{\partial\bm{x}_{0}}\in\mathbb{R}^{n\times
n}$ defines the variational mapping function, such that
$\bm{\Phi}_{0}^{0}(\bm{x}_{0})=\bm{I}_{n}$ and matrix
$\bm{I}_{n}\in\mathbb{R}^{n\times n}$ is an identity matrix. Readers are
referred to [41, 42] for the derivation of the variational system (6). The
variational system (6) depicts how small system disturbances evolve along the
system trajectory. That is, (6) allows the quantification the propagation of
uncertain load disturbances within the power system. For ease of notation,
moving forward we remove the dependency of
$\bm{\Phi}_{0}^{k}(\bm{x}_{0})=\bm{\Phi}_{0}^{k}$ on $\bm{x}_{0}$.

######  Remark 1.

Notice that $\bm{\Phi}_{0}^{k}$ represents the derivative of (5) with respect
to $\bm{x}_{0}$ for $k\in\\{0\;,1\;,\;\ldots\;,\mathrm{N}-1\\}$. This being
said, the transition matrix $\bm{\Phi}_{0}^{k}$ requires the knowledge of
$\bm{x}_{k}$ for all $k$. As such, we can apply the chain rule to evaluate
$\bm{\Phi}_{0}^{k}$ for any time-index $k$ as

$\bm{\Phi}_{0}^{k}=\bm{\Phi}_{k-1}^{k}\bm{\Phi}_{k-2}^{k-1}\;\ldots\;\bm{\Phi}_{0}^{1}\bm{\Phi}_{0}^{0}=\prod^{i=k}_{1}\bm{\Phi}^{i}_{i-1}.$
(7)

The subsequent assumption is introduced to ensure that the system’s response
to perturbations or uncertainties remains bounded along the system trajectory.
This is equivalent to considering a bound on the Jacobian of the nonlinear
dynamical system, that is,
$\left\lVert\bm{J}(\tilde{\bm{f}}(\cdot))\right\rVert=\left\lVert\frac{\partial\tilde{\bm{f}}(\cdot)}{\partial\bm{x}_{k}}\right\rVert<\infty$;
see [43].

######  Assumption 2.

Consider the variational system (6), the state-transition matrix
$\bm{\Phi}_{0}^{k}$ is bounded, i.e.,
$\sup\left\\{\;\|\bm{\Phi}_{0}^{k}\|:k\in\mathbb{N}\;\right\\}<\infty.$

It follows that for every $\bm{\delta}\bm{x}_{0}\in\mathbb{R}^{n}$ there
exists a unique solution for
$\bm{\delta}\bm{x}_{k}\;\forall\;k\in\\{0\;,1\;,\;\ldots\;,\mathrm{N}-1\\}$.
Assumption 2 is a direct implication of the regularity condition of the
nonlinear system (1) mentioned in Assumption 1. To that end, the maximal LE in
finite-time stating from initial state $\bm{x}_{0}\in\mathbb{R}^{n}$ can be
defined as

$\leavevmode\nobreak\
{\lambda}_{\mathrm{MLE}}:=\sup\left\\{\;\lim_{k\rightarrow\mathrm{N-1}}\frac{1}{k}\log\left(\dfrac{\left\lVert\bm{\delta}\bm{x}_{k}\right\rVert}{\left\lVert\bm{\delta}\bm{x}_{0}\right\rVert}\right)\;\right\\},$
(8)

where ${\lambda}_{\mathrm{MLE}}$ is the largest eigenvalue of time-evolution
of variational state-vector $\bm{\delta}\bm{x}_{k}$ with respect to the
initial perturbation $\bm{\delta}\bm{x}_{0}$, i.e., the spectral norm of
Cauchy-Green tensor defined as (9) under the action of the logarithm function;
see [44, 45].

$\bm{\Xi}:={\bm{\Phi}_{0}^{k}}^{\top}{\bm{\Phi}_{0}^{k}};$ (9)

where matrix $\bm{\Xi}\in\mathbb{R}^{n\times n}$ represents the deformation of
the initial perturbation vector $\bm{\delta}\bm{x}_{0}$ along the trajectory
of the system. The sign of the MLE ${\lambda}_{\mathrm{MLE}}$ indicates the
exponential divergence or convergence about a state orbit; see Figure 1. The
exponent is related to the average expansion and contraction along the
directions of the state-space, such that, for any initial perturbation
$\bm{\delta}\bm{x}_{0}$ the stability of the system can be defined according
to the following

  * •

Unstable if ${\lambda}_{\mathrm{MLE}}>0$, i.e., system is chaotic and
diverges.

  * •

Stable if ${\lambda}_{\mathrm{MLE}}\leq 0$, i.e., system trajectory is
attracted to a stable or periodic orbit.

The maximal LE throughout the literature is typically applied to study a
single aspect of stability, that being frequency, voltage, or rotor angle
stability. In this paper, we present a method for stability assessment that
considers all the above aspects simultaneously. As such, the next sections
introduce LEs-based stability through computing the full spectrum of exponents
resulting from different types of stability.

$({a})$$\bm{x}_{0}$$\bm{x}_{0}+\bm{\delta}\bm{x}_{0}$$\bm{x}_{t}$$\bm{x}_{t}+\bm{\delta}\bm{x}_{t}$$\bm{\lambda_{MLE}}\leq
0$ (a)

$({b})$$\bm{x}_{0}$$\bm{x}_{0}+\bm{\delta}\bm{x}_{0}$$\bm{x}_{t}$$\bm{x}_{t}+\bm{\delta}\bm{x}_{t}$$\bm{\lambda_{MLE}}>0$
(b)

Figure 1: Trajectory of nearby orbits starting from initial state $\bm{x}_{0}$
and perturbed initial state $\bm{x}_{0}+\bm{\delta}\bm{x}_{0}$: $(a)$
converging and $(b)$ diverging trajectories.

###  III-B Computation of Lyapunov spectrum of exponents

The MLE is sufficient to determine the asymptotic stability of a dynamical
system. In addition to asymptotic stability assessment, computing the spectrum
of LEs provides a characterization of the behavior of system along all
direction of the phase-space. That being said, the spectrum provides an
identification of the direction or nodes of a system that are sensitive to a
perturbation in initial conditions, i.e., it allows for the understanding of
uncertainty propagation along the dynamic trajectory of a system. The spectrum
of LEs can be computed as follows

$\displaystyle\bm{\Lambda}$
$\displaystyle:=\mathrm{spec}\left\\{\lim_{k\rightarrow\mathrm{N-1}}\frac{1}{k}\log\left(\left\lVert{\bm{\Phi}_{0}^{k}}\right\rVert\right)\right\\},\leavevmode\nobreak\
$ (10a)
$\displaystyle=\mathrm{spec}\left\\{\lim_{k\rightarrow\mathrm{N-1}}\frac{1}{2k}\log\left(\underbrace{{\bm{\Phi}_{0}^{k}}^{\top}{\bm{\Phi}_{0}^{k}}}_{\bm{\Xi}}\right)\right\\},$
(10b)

where $\bm{\Lambda}\in\mathbb{R}^{n\times n}$ is a diagonal matrix
representing the spectrum of LEs. Note that for a finite dimensional vector
space the spectrum coincides with the eigenvalues of the matrix. This
indicates that the spectrum of LEs are exactly the eigenvalues corresponding
to the covariant Lyapunov vectors. For any regular dynamical system, these
vectors represent the local decomposition of the phase space [46]. The local
covariant vectors along with the exponents inform us on the local stability
behavior along the phase trajectory of the system [47, 2].

######  Remark 2.

The assumption pertaining to the regularity of the dynamical system (1) is
mild. The existence of the full spectrum of LEs is well-established as a
result of the multiplicative ergodic theorem (Oseledets Theorem) proven in the
1960s; see [48].

The existence of the full spectrum of LEs guarantees regularity based on
Oseledets’ Theorem. The reason is that computing LEs requires the computation
of well-defined Jacobian matrices of the nonlinear system. These Jacobian
matrices correspond to the variational mapping function in (6), whose
existence is ensured by the properties guaranteed by Oseledets’ Theorem.
Therefore, regularity, defined by the existence and smoothness of such
Jacobian matrices, is implied. It follows that for general discrete-time
systems regularity is provided as a result of the aforementioned theorem [49,
50, 51, 52]. This result eliminates the technical challenges of computing the
LEs while requiring to verification of the regularity assumptions.

Computation of the spectrum of LEs is well-established [44, 53, 54, 55, 56].
However, generally the fundamental mapping function $\bm{\Phi}_{0}^{k+1}$ for
the variational system (6) is ill-conditioned. Therefore, the error in
computing the LEs increases and thus the solution tends to converge in the
direction of the MLE [57, 50, 58]. To avoid such problem when practically
computing the spectrum of exponents, one can rely on re-orthogonalization of
local directions of the mapping function $\bm{\Phi}_{0}^{k+1}$ along the
system trajectory. This results in a transformation of $\bm{\Phi}_{0}^{k+1}$
to an upper rectangular matrix. The existence of such triangular matrix for a
regular dynamical system is a direct result of Perron’s lemma [48, Lemma
1.3.3].

There are two orthogonal factorization classes for estimating the LEs from the
triangular matrix form (i) continuous QR- and (ii) discrete QR-method. We
compute the spectrum of LEs by utilizing the discrete QR-method; see [54]. The
matrix $\bm{\Phi}^{k}_{0}$ and its triangular factor $\bm{R}^{k}_{0}$ are
directly evaluated by a re-orthogonalization integration via discrete-QR
orthonormalization. For $k\in\\{0,\;1,\ldots,\;\mathrm{N}-1\\}$ the discrete-
QR factorization of transition matrix $\bm{\Phi}^{k}_{0}$ (7) can be written
as

$\bm{\Phi}^{k}_{0}=\bm{Q}_{k}\bm{R}^{k}_{0},\quad\bm{\Phi}^{0}_{0}=\bm{Q}_{0}\bm{R}^{0}_{0}=\bm{I}_{n},$
(11)

where $\bm{Q}_{k}\in\mathbb{R}^{n\times n}$ is an orthogonal matrix and
$\bm{R}^{k}_{0}\in\mathbb{R}^{n\times n}$ is a positive upper triangular
matrix. Due to the uniqueness of the solution of the QR-factorization, as a
result of regularity, the triangular matrix $\bm{R}^{k}_{0}$ can be written as

$\bm{R}_{0}^{k}=\bm{R}_{k-1}^{k}\bm{R}_{k-2}^{k-1}\ldots\bm{R}_{0}^{1}\bm{R}_{0}^{0}=\prod^{i=k}_{1}\bm{R}^{i}_{i-1}.$
(12)

With that in mind, to attenuate the ill-conditioned matrix computations from
(10), the spectrum of LEs can be expressed as

$\displaystyle\bm{\Lambda}$
$\displaystyle:=\mathrm{spec}\left\\{\lim_{k\rightarrow\mathrm{N-1}}\frac{1}{k}\log\left\lvert\prod^{i=k}_{1}\bm{R}^{i}_{i-1}\right\rvert\right\\},\leavevmode\nobreak\
$ (13a)
$\displaystyle=\mathrm{spec}\left\\{\lim_{k\rightarrow\mathrm{N-1}}\frac{1}{k}\sum^{i=k}_{1}\log\left\lvert\bm{R}^{i}_{i-1}\right\rvert\right\\}.$
(13b)

######  Remark 3.

Notice that the diagonal elements of $\bm{R}_{0}^{k}$ are only required for
the computation of the spectrum of LEs. It follows that to compute LEs, we
evaluate $\left\lvert\bm{R}_{0}^{k}\right\rvert$ instead of
$\left\lVert\bm{R}_{0}^{k}\right\rVert$; see [55].

Computing the LEs enables the formulation of a scalable method that assesses
stability of NL-DAE system against negative load injections renewables and to
potentially inform the power system operators of what renewables injections
could cause in terms of stability.

##  IV Stable node identification and optimal renewables allocation

The previous section introduced Lyapunov’s method for the stability analysis
of dynamical systems. In this section, we present a framework based on the
aforementioned stability quantification method for identifying stable nodes in
NL-DAE power systems. The proposed framework enables the allocation of
uncertain and intermittent loads from RER while considering the impact on
power system overall stability.

###  IV-A Stable node identification

The high penetration of renewables, along with their dynamic behavior and
lower system inertia, impose a challenge on the transient stability of the
system [9]. Typically, methods that assess stability based on the LEs of a
power system consider model-free computations. This implies that the LEs are
computed based on the voltage stability of the buses. However, although some
methods mentioned in this paper consider rotor angle stability, they rely on
state estimation. This indeed results in computational errors when estimating
the LEs that measure the infinitesimal trajectory of state perturbations. In
this paper, having provided a model-based approach for computing the LEs of a
NL-DAE power system representation, we present a framework to evaluate node
stability by considering a more comprehensive approach.

With that in mind, we assess the stability of the nodes, i.e., buses, of a
power system based on the several stability criteria: (i) voltage stability,
(ii) rotor angle stability, and (iii) frequency stability. The stable node
identification framework is based on computing the LEs while considering
system dynamic response after a RER load disturbance. To study the stability
of a power network after a load disturbance resulting from such intermittent
and uncertain RER injections, we introduce a quantitative measure equivalent
to computing the spectrum of LEs of the power system. The following
proposition establishes this expression.

######  Proposition 1.

The parametrized tensor matrix (9) representing the state-deformation along
the trajectory of the nonlinear discrete-time power system (5) can be
expressed as follows

$\tilde{\bm{\Xi}}(\bm{\gamma}):=\sum_{j=1}^{n}\gamma_{j}\left(\sum_{i=0}^{k}\left(\bm{\varphi}^{i}_{0}\right)^{\top}\bm{\varphi}^{i}_{0}\right)\;\in\;\mathbb{R}^{n\times
n},$ (14)

where $\bm{\varphi}^{i}_{0}$ represent the column vectors of matrix
$\bm{\Phi}_{0}^{k}$. The parametrization $\gamma_{j}$ determines the states
that are required for the stability assessment. That is, if $\gamma_{j}=1$
then the state is considered for LEs computations, else $\gamma_{j}=0$. The
parametrized vector $\bm{\gamma}\in\mathbb{R}^{n}$ that represents the
selected states is denoted as $\left\\{{\gamma}_{j}\right\\}_{{j}={1}}^{n}$.

###### Proof.

From (9) being multiplied with parametrization vector $\bm{\gamma}$ it follows
that

$\displaystyle\tilde{\bm{\Xi}}(\bm{\gamma})$
$\displaystyle={\bm{\Phi}^{k}_{0}}^{\top}\begin{bmatrix}\bm{I}\otimes\bm{\gamma}\end{bmatrix}^{\top}\begin{bmatrix}\bm{I}\otimes\bm{\gamma}\end{bmatrix}\bm{\Phi}^{k}_{0},\leavevmode\nobreak\
$ (15a)
$\displaystyle=\sum_{i=0}^{k}\left({\bm{\varphi}^{i}_{0}}\right)^{\top}\bm{\gamma}^{2}\bm{\varphi}^{i}_{0}=\sum_{i=0}^{k}\sum_{j=1}^{n}\gamma_{j}\left(\bm{\varphi}^{i}_{0}\right)^{\top}\bm{\varphi}^{i}_{0},\leavevmode\nobreak\
$ (15b)

where (15a) is a result of applying the dot product and the Kronecker product,
and (15b) holds true since $\bm{\gamma}^{2}=\bm{\gamma}$. The proof is
complete since (15b) is equivalent to (14). ∎

We note that throughout this section we refer to the deformation matrix (9)
instead of (13) for developing the overall stability quantification measures;
although numerically in the case studies we implement the algorithms by
computing the spectrum of LEs based on the QR-factorization mentioned above.
The reason is to simplify the proofs related to how we develop the stability
measures and their relation to LEs. Such that under the QR factorized
matrices, an unnecessary layer of complexity will be added to the proofs.

The above parametrized deformation matrix (14) represents the deformation
along the trajectory for the selected states. Such a parametrization allows
for choosing, at each bus, the states that contribute to the stability
computation. If a bus is a generator bus then we can choose from the
differential states (2a) and algebraic state vector (2b). If it is a load bus
we can choose from the algebraic state vector (2b). This allows for a given
bus to choose a subset of states for the stability quantification, however
when selecting more than one state, we obtain several LEs for each of the
buses. To present a singular LE for a bus, we first introduce some properties
for LEs that are necessary for the unified quantification.

######  Property 1 ([48]).

LEs $\bm{\lambda}$ exhibit the following properties

$\displaystyle(\mathrm{P}1.1)$
$\displaystyle\quad\bm{\lambda}(\beta\bm{A})=\bm{\lambda}(\bm{A})\;\forall\;\beta\in\mathbb{R}\backslash\\{0\\},$
$\displaystyle(\mathrm{P}1.2)$
$\displaystyle\quad\bm{\lambda}(\bm{A}_{1}+\bm{A}_{2})\leq\max\\{\bm{\lambda}(\bm{A}_{1}),\bm{\lambda}(\bm{A}_{2})\\},$

where $\bm{A}_{1}$ and $\bm{A}_{2}$ are of dimension $\mathbb{R}^{n\times n}$;
see [48, Theorem 2.1.2].

Based on the above properties, we now present the Lyapunov-based bus stability
identification in power networks. The following proposition establishes such
stability quantification method.

######  Proposition 2.

The stability of a node in a power network characterized by its LE can be
expressed as (16). Let $i\;\in\;\mathcal{N}$ denote the bus index, such that
$i_{\mathcal{G}}\;\in\;\mathcal{G}$ and
$i_{\mathcal{L}}\;\in\;\mathcal{L}\cap\mathcal{R}$ represent the generator and
load/renewable buses.

$\leavevmode\nobreak\
\lambda_{i}:=\begin{cases}\lambda_{i_{\mathcal{G}}}=\underset{k\rightarrow\mathrm{N-1}}{\lim}\left[\frac{1}{2k}\log\left(\tilde{\bm{\Xi}}(\bm{\gamma})\right)\right],&\gamma_{j}\in\left\\{\delta_{i},\omega_{i},v_{i}\right\\}_{i\;\in\;\mathcal{G}},\\\
\lambda_{i_{\mathcal{L}}}=\underset{k\rightarrow\mathrm{N-1}}{\lim}\left[\frac{1}{2k}\log\left(\tilde{\bm{\Xi}}(\bm{\gamma})\right)\right],&\gamma_{j}\in\left\\{v_{i}\right\\}_{i\;\in\;\mathcal{L}},\end{cases}$
(16)

where $j\in\mathbb{R}^{n}$ represents the index for the differential and
algebraic states defined by state vector $\bm{x}_{k}$.

###### Proof.

For the proof we consider a generator bus $i\;\in\;\mathcal{G}$. Let
$\tilde{\bm{\Xi}}(\bm{\gamma})$ define the deformation matrix of the power
grid. From Proposition 1, we can write the following

$\tilde{\bm{\Xi}}(\bm{\gamma})=\sum_{i=0}^{k}\sum_{j=1}^{n}\gamma_{j}\left(\bm{\varphi}^{i}_{0}\right)^{\top}\bm{\varphi}^{i}_{0},\quad\gamma_{j}\in\left\\{\delta_{i},\omega_{i},v_{i}\right\\}_{i\;\in\;\mathcal{G}},$
(17)

where applying the $\mathrm{log}$ function to the above matrix we obtain the
following

$\lambda_{i_{\mathcal{G}}}=\frac{1}{2k}\mathrm{log}\left(\sum_{i=0}^{k}\sum_{j=1}^{n}\gamma_{j}\left(\bm{\varphi}^{i}_{0}\right)^{\top}\bm{\varphi}^{i}_{0}\right).$
(18)

Taking the limit of the $\mathrm{log}$ and applying properties
$(\mathrm{P}1.1)$ and $(\mathrm{P}1.2)$ we obtain the Lyapunov exponents of
the power network buses for $i\;\in\;\mathcal{G}$. This holds true for
load/renewables buses $i\;\in\;\mathcal{L}\cap\mathcal{R}$. ∎

The above proposition establishes the stability of a node in a power grid; it
is based on the type of bus and thus the states for the stability assessment
represented by parametrization $\gamma_{j}$. For a generator bus
$i\in\mathcal{G}$, $\gamma_{j}$ represents rotor angle, frequency and voltage
stability, whereas for a load/renewables bus $i\in\mathcal{L}$, $\gamma_{j}$
represents voltage stability.

######  Remark 4.

We note here that for generator buses, i.e., $i\in\mathcal{G}$, the state
variables
$\left\\{\delta_{i},\;\omega_{i},\;v_{i}\right\\}_{i\;\in\;\mathcal{G}}$
contribute to the computation of the nodes’ LEs. These exponent are bounded by
the maximum exponents computed from each of the variables, a result of
Property $(\mathrm{P}1.2)$. That being said, the following upper bound holds
true

$\leavevmode\nobreak\
\lambda_{i_{\mathcal{G}}}\leq\max\\{\lambda_{i_{\mathcal{G}}}(\left\\{\delta_{i}\right\\}_{i\in\mathcal{G}}),\;\lambda_{i_{\mathcal{G}}}(\left\\{\omega_{i}\right\\}_{i\in\mathcal{G}}),\;\lambda_{i_{\mathcal{G}}}(\left\\{v_{i}\right\\}_{i\in\mathcal{G}})\\}.$
(19)

Having presented the above node stability quantification method, we now rank
the nodes according to their corresponding LE value. As discussed earlier, a
stable state trajectory has a negative LE. As such, let $\mathrm{S}_{i}$
define the stability index of a node in a power system network. The stability
index of a node is equivalent to the index of the LE from the set of ordered
exponents in ascending order and can be expressed as

$\mathrm{S}_{i}:=\mathrm{index}(\lambda_{i}),\quad\lambda_{i}\in\tilde{\bm{\Lambda}},$
(20)

where
$\tilde{\bm{\Lambda}}=\mathrm{col}\left\\{\lambda_{i}\right\\}_{i\in\mathcal{N}}$
represents a column vector of the ordered LEs $\lambda_{i}$ in ascending
order. The next section introduces the overall stability metric that
quantifies a renewable load’s impact on the power grid as a whole. This metric
indicates the uncertainty propagation onto the power grid from a RER load
injection that is being applied to a single node within the network.

###  IV-B Quantifying stability against renewables and optimally allocating
their locations

The method for assessing the stability of generator and load buses within a
network relies on estimating the LEs and then ranking the stability of the
buses according to the value of the exponents. That is, buses with larger
negative LEs are considered more stable and thus are given a higher stability
index $\mathrm{S}_{i}$. The stability index that is based on the rank, under
the context of renewables allocation, does not assess the impact of such
renewable load injection, along with its uncertainty propagation, on the
overall stability of the network.

With that in mind, we approach allocating RER by considering a LE-based
stability measure that enables studying the impact of renewables injections on
the overall stability of the power network. Meaning that, we want to quantify
the impact of an uncertain renewable injection on the overall stability of the
power system and then choose the nodes that least impact the overall system
stability. Before posing the RER allocation problem within power networks that
are described by a complete NL-DAE representation, we introduce the following
stability measure that is equivalent to computing the sum of LEs along all the
nodes of the power network. The LEs-based measure quantifies the uncertainty
that is propagated with the power grid from a load injection at a network
node.

######  Theorem 1.

The $\mathrm{log\,det}$ of matrix $\tilde{\bm{\Xi}}(\bm{\gamma})$ quantifies
the overall stability of a system. As such, the
$\mathrm{log\,det}(\tilde{\bm{\Xi}}(\bm{\gamma}))$ is said to be equivalent
the Lyapunov spectrum of exponents (13) according to the following relation

$\mathrm{log\,det}(\tilde{\bm{\Xi}}(\bm{\gamma}))\equiv\beta\sum_{i=1}^{N}\lambda_{i},\quad\forall\;i\in\mathcal{N},$
(21)

where $\beta$ is a constant equivalent to $\tfrac{1}{2\mathrm{N}}$. The
constant $\beta$ is due to the definition of LEs being computed over the
trajectory.

###### Proof.

From (16) observe that the LEs are the eigenvalues of the $\mathrm{log}$ of
matrix $\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})$. This matrix is positive
semi-definite, whereby its determinant can be computed as

$\displaystyle\mathrm{det}(\tilde{\bm{\Xi}}(\bm{\gamma}))$
$\displaystyle=\left(\prod_{i=1}^{N}\hat{\lambda}_{i}\right),$ (22)

where $\hat{\lambda}_{i}$ is the $i$-$th$ eigenvalue of matrix
$\tilde{\bm{\Xi}}(\bm{\gamma}))$. Now taking the $\mathrm{log}$ of the
eigenvalues $\hat{\lambda}_{i}$, we obtain according to (16) the $i$-$th$
Lyapunov exponent. As such applying $\mathrm{log}$ to (22) we obtain the
following

$\displaystyle\mathrm{log\,det}(\tilde{\bm{\Xi}}(\bm{\gamma}))$
$\displaystyle=\beta\left(\sum_{i=1}^{N}{\lambda}_{i}\right),$ (23)

this concludes the proof. ∎

The above theorem presents a method for stability quantification after a load
injection at a system node; it suggests that computing (21) is equivalent to
computing the overall stability of all the nodes in a power network after a
renewable load is applied to a single node. As such, this metric enables the
quantification of the impact on stability, through computing the spectrum of
LEs, resulting from the allocation of uncertain renewables at a node within a
power network.

In order to identify the set of nodes that result in minimal uncertainty
propagation, i.e., nodes that are indicative of stability when a renewable
load injection is applied to it, we are required to solve a combinatorial
optimization problem. With that in mind, we now pose the optimal RER
allocation problem based on the overall stability impact, computed based on
(21), of a renewable injection along the buses within a power network. This
combinatorial problem is computationally expensive; it increases in complexity
for large networks. As such, we pose the above problem as a set optimization
problem for which we provide a computationally efficient algorithm to solve.
With that in mind, let $\mathcal{S}$ denote the set of RER to be allocated
within an existing power network. Let the maximum the number of RER be $N$. As
such, the maximum number of RER is limited to the number of buses within a
power network.

The optimal RER allocation problem P1 for the NL-DAE power system (5) can be
written as a set optimization problem by defining the set function
$\mathcal{L}{(\mathcal{S})}:2^{\mathcal{V}}\rightarrow\mathbb{R}$ as the
$\mathrm{log\,det}$ of parametrized matrix
$\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})$, where
$\mathcal{V}:=\\{i\in\mathbb{N}\,\;|\,\;\;0<i\leq N\\}$.

$\displaystyle(\mathrm{\textbf{P1}})\;\;\operatorname*{maximize}_{\mathcal{{\bm{S}}}}\;\;\;$
$\displaystyle\mathcal{L}{(\mathcal{S})}:=\mathrm{log\,det}\left(\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\right),\leavevmode\nobreak\
$ (24a) $\displaystyle\operatorname*{subject\ to}\;\;\;$
$\displaystyle\left\lvert\mathcal{S}\right\rvert=s,\;\;\mathcal{S}\subseteq\mathcal{V},$
(24b)

where $\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\;\in\;\mathbb{R}^{n\times n}$
represents the state-trajectory deformation matrix resulting from allocating
an uncertain renewable load injection at a node defined by set
$\mathcal{{S}}$. As such, $\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})$ can be
defined as

$\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma}):=\sum_{i=1}^{N}\alpha_{i}\tilde{\bm{\Xi}}_{i}(\bm{\gamma}),\quad\alpha_{i}\;\in\;\\{0,1\\}^{\mathcal{S}},$
(25)

where $\alpha_{i}$ is equivalent to zero if set $\mathcal{S}$ has a zero at
the $i$-$th$ node index, and is one otherwise. This means that
$\bm{\alpha}:=\\{\alpha_{i}\\}_{i\in\mathcal{N}}$ represents the allocation of
RER within the power network, where a value of one is given to a node with an
RER injection. Matrix $\tilde{\bm{\Xi}}_{i}(\bm{\gamma})$ denotes the
deformation of state trajectories when an RER load injection is applied to a
single node of index $i\in\mathcal{N}$. We note that based on Theorem 1, the
optimal allocation problem P1 is equivalent to computing the sum of the LEs
across all the nodes after the application of a renewable load injection at a
single perturbed node. As such, at each iteration of the allocation problem
P1, we solve the maximum value for $\mathcal{L}(\mathcal{S})$ by choosing the
node that has the largest sum of LEs. This allocation problem will result in
an optimal set $\mathcal{S}^{*}$ representing the ranked nodes from most
critical to most stable. This is due to the fact that we require at each node
to have a negative LE. The solution provides the operator the nodes that are
indicative of stability; it informs the operator about which nodes in the
system, when allocated a renewable load injection, lead to an unstable or
stable system. Accordingly, to obtain the most stable node, we set
$|\mathcal{S}|=N$ and choose the most stable node to be the last index of set
$\mathcal{S}^{*}$.

Notice that the optimal problem P1 is a combinatorial set optimization
problem, and thus it is computationally expensive. To avoid the complexity of
utilizing global solvers for solving the combinatorial problem, we exploit the
submodular*** A set function
$\mathcal{L}:2^{\mathcal{V}}\rightarrow\mathbb{R}$ is said to be submodular if
and only if for any $\mathcal{A},\mathcal{B}\subseteq\mathcal{V}$ with
$\mathcal{A}\subseteq\mathcal{B}$, it holds true that
$\mathcal{L}(\mathcal{A}\cup\\{s\\})-\mathcal{L}(\mathcal{A})\geq\mathcal{L}(\mathcal{B}\cup\\{s\\})-\mathcal{L}(\mathcal{B})$.
and monotone increasing†††A set function
$\mathcal{L}:2^{V}\rightarrow\mathbb{R}$ is monotone increasing if
$\forall\;A,B\subseteq V$, it holds true that $A\subseteq
B\rightarrow\mathcal{L}(A)\leq\mathcal{L}(B)$ properties of the
$\mathrm{log\,det}$ set function that allows for solving P1 using greedy
algorithms. The following proposition establishes the submodularity of
$\mathrm{log\,det}$ set function.

######  Proposition 3.

The $\mathrm{log\,det}$ of (25) denoted by set function
$\mathcal{L}(\mathcal{S}):2^{\mathcal{V}}\rightarrow\mathbb{R}$ as follows

$\leavevmode\nobreak\
\mathcal{L}(\mathcal{S})=\mathrm{log\,det}\left(\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\right),$
(26)

is, for $\mathcal{S}\subseteq\mathcal{V}$, submodular and monotone increasing.

The complete proof of Proposition 3 is available in Appendix E. The
aforementioned proposition establishes the submodularity of the optimal
allocation problem P1, this allows the use of greedy algorithms to solve the
computationally exhaustive problem in a more efficient method. The greedy
algorithm [59, Algorithm 1] is utilized to solve the aforementioned submodular
optimization problem. The optimality guarantee of the greedy algorithm for
submodular set optimization is detailed in Appendix D. We note that this
optimality guarantee reaches $99\%$ optimality when solving an allocation
problem that has a submodular objective function.

1 input: Network parameters, $s$, $i_{\mathcal{G}}$, $i_{\mathcal{N}}$,
$\mathcal{V}$, $\mathcal{N}$

2 initialize:
$({\tilde{P}}_{\mathrm{R_{i}}}^{0},{\tilde{Q}}_{\mathrm{R_{i}}}^{0})\leftarrow(\bm{0},\bm{0})$,
$i\leftarrow 1$, $k\leftarrow 1$

3 for  _$i=1$ to $N$_ do

4 assign:
$({\tilde{P}}_{\mathrm{R_{i}}}^{0},{\tilde{Q}}_{\mathrm{R_{i}}}^{0})\leftarrow(1+\tfrac{\beta}{100})(\mathrm{P}_{\mathrm{R_{i}}}^{0},\mathrm{Q}_{\mathrm{R_{i}}}^{0})$

5 for  _$k=\;1$ to $\mathrm{N}-1$_ do

6 simulate: (4) by implementing Algorithm 2

7 compute: $\tilde{\bm{\Xi}}_{i}(\bm{\gamma})$ by implementing Algorithm 3

8

9

10initialize: $\mathcal{S}\leftarrow\emptyset$, $i\leftarrow 1$,
${\alpha_{i}}\leftarrow{0}\;\forall\;i\;\in\;\mathcal{N}$

11 for  _$i\leq s$_ do

12 set:
$\mathcal{L}{(\mathcal{S})}=\mathrm{log\,det}\left(\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\right)=\mathrm{log\,det}\left(\sum_{i=1}^{N}\alpha_{i}\tilde{\bm{\Xi}}_{i}(\bm{\gamma})\right)$

13 compute:
$\mathcal{G}_{i}=\mathcal{L}(\mathcal{S}\cup\\{a\\})-\mathcal{L}(\mathcal{S})$,
$\forall a\in\mathcal{V}\setminus\mathcal{S}$

14 assign:
$\mathcal{S}\leftarrow\mathcal{S}\cup\left\\{\mathrm{arg\,max}_{a\in\mathcal{V}\setminus\mathcal{S}}\,\mathcal{G}_{i}\right\\}$

15 update: $i\leftarrow i+1$

16

17output: $\mathcal{S}^{*}$

Algorithm 1 Framework for identifying optimal RER injection locations within a
power network

The proposed problem P1 can be solved efficiently according to the greedy
algorithm. The framework for identifying optimal stable nodes and uncertain
renewable injection locations is outlined in Algorithm 1 which embeds the
greedy algorithm to obtain $\mathcal{S}^{*}$.

##  V Numerical Case Study

In this section, we demonstrate the proposed stable node identification and
RER allocation framework on standard power systems. The goal here is to
identify stable nodes that allow for attenuation of the uncertain and
intermittent loads from RER on the overall stability of a power network. With
that in mind, we attempt to answer the following questions.

  * •

$(\mathrm{Q}1)$ Does the stability index $\mathrm{S}_{i}$ depict different
forms of system instability (frequency, voltage and rotor angle) in power
systems?

  * •

$(\mathrm{Q}2)$ Does solving the optimal RER allocation framework result in
optimal stable nodes that offer the least overall stability implication on the
power network?

  * •

$(\mathrm{Q}3)$ Is the clearing time after inducing a disturbance from an RER
load at an optimal node shorter than disturbances added at other nodes? This
verifies the optimality of the allocation problem.

  * •

$(\mathrm{Q}4)$ Is the overall stability measure numerically equivalent to the
summation of LEs (see Theorem 1)?

###  V-A Implementation of the stability and allocation framework

(a)  (b)  (c)  (d)

Figure 2: Lyapunov spectrum of exponents for (a) $\mathrm{case}$-$\mathrm{9}$
(9-buses) and (b) $\mathrm{case}$-$\mathrm{200}$ (200-buses). The (left)
column depicts the spectrum of LEs computed for all the system states. The
(right) column depicts the LEs for each bus within the network, i.e., the
stability index.

We consider two standard power systems of contrasting size for the assessment
of the proposed framework: (i) $\mathrm{case}$-$\mathrm{9}$: A
$\mathrm{9}$-bus power system with $\mathrm{3}$ synchronous generators—Western
System Coordinating Council (WSCC), and (ii) $\mathrm{case}$-$\mathrm{200}$: A
$\mathrm{200}$-bus power system with $\mathrm{49}$ synchronous
generators—ACTIVSg200-Bus network ”Illinois200 case”. The simulations and
optimization problem are performed in MATLAB R2021b running on a Macbook Pro
having an Apple M1 Pro chip with a 10-core CPU and 16 GB of RAM.

The generator parameters are extracted from power system toolbox (PST) [22]
case file $\mathrm{data3m9b.m}$ for $\mathrm{case}$-$\mathrm{9}$. For
$\mathrm{case}$-$\mathrm{200}$ the generator parameters are chosen based on
the ranges provided in the PST toolbox. Regulation and chest time constants
for the generators are chosen as $R_{\mathrm{D}i}=0.2\;\mathrm{Hz/pu}$ and
$T_{\mathrm{CH}i}=0.2\;\mathrm{sec}$. The synchronous speed is set to
$\omega_{0}=120\pi\;\mathrm{rad/sec}$ and a power base of $100\;\mathrm{MVA}$
is considered for the power system. The steady-state initial conditions for
the power system are generated from solutions of the power flow obtained from
MATPOWER $\mathrm{runpf}$ function.

The discretization constant for the step-size is $h=0.1$ and simulation time-
span $t=30\;\mathrm{sec}$. We simulate a renewable load disturbance by
introducing at $t>0$ a RER load denoted as
$(\mathrm{P}_{\mathrm{R}}^{0},\mathrm{Q}_{\mathrm{R}}^{0})$ using the NR
method as described in Appendix B. In the scope of this work, RERs are modeled
as negative loads that are injected into the power network at a node in set
$\mathcal{N}$. Uncertainty of an RER load is modeled by adding a load
perturbation of magnitude $\beta$. The perturbed renewables load at a node is
then computed as
$(\mathrm{\tilde{P}}_{\mathrm{R}}^{0},\mathrm{\tilde{Q}}_{\mathrm{R}}^{0})=(1+\tfrac{\beta}{100})(\mathrm{P}_{\mathrm{R}}^{0},\mathrm{Q}_{\mathrm{R}}^{0})$.
To account for the uncertainty of the renewable load we vary $\beta$ between
$\\{2\%,20\%\\}$. To verify accuracy of the proposed NL-DAE transformation
readers are referred to [33]. Note that the power generator model’s (27)
response to the transients induced by the uncertain renewable injections is
automatically regulated by the governor response $\bm{T}_{\mathrm{r}}$.

To quantify the overall stability impact on a power network after a renewable
injection is allocated at a bus (whether a generator or load bus; see,
Proposition 2), we solve the optimal problem P1. The horizon window for LEs
computations is chosen as $\mathrm{N}=t/h=300.$ First, we quantify the
parametrized deformation matrix $\tilde{\bm{\Xi}}_{i}(\bm{\gamma})$ (see,
Proposition 1) for each $i\in\mathcal{N}$. Here, a renewable load injection
along with its uncertainty is applied to a node within the network as
described above. To obtain accurate computations of
$\tilde{\bm{\Xi}}_{i}(\bm{\gamma})$ (see, Section III-B), we utilize Algorithm
3. Then, we utilize Theorem 1 to compute the objective function of P1, such
that the
$\mathrm{log\,det}\left(\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\right)$ is
equivalent to computing the sum of LEs of the system after a renewable
injection is applied at a node $i\in\mathcal{N}$. Based on whether a node is a
generator or load bus, we obtain the LEs be referring to Proposition 2. The
greedy algorithm is then adopted to solve P1 for $|\mathcal{S}|=s$. Here we
want to quantify the impact of stability when a renewable load is applied at
each of the network buses, as such we set $s=N$. The solution to the optimal
allocation problem, denoted by $\mathrm{\mathbf{P1}}$, is the set
$\mathcal{S}^{*}$, where $\mathcal{S}^{*}$ represents an ordered set of nodes.
The nodes are ranked based on their ability to minimize uncertainty
propagation in the system when uncertain renewable loads are allocated there.
This detailed framework is summarized in Algorithm 1.

###  V-B Stability index and node ranking for potential RER allocation

As mentioned in the aforementioned sections, the validity of the NL-DAE model
and the discretization is studied in [33]. With that in mind, we now want to
analyze the stability of the nodes through quantifying the LEs of the network
that already has renewable load injections allocated at random buses. The LEs
for $\mathrm{case}$-$\mathrm{9}$ and $\mathrm{case}$-$\mathrm{200}$ computed
using (13), considering all the differential and algebraic states, are
depicted in left column of Figure 2. The states of the systems are summarized
in (2a)–(2b). Notice that the LEs corresponding to each state range from
positive LEs to negative LEs indicating that the system is inundated with
transients from the allocated renewables. The LEs that characterize a bus
stability are then computed for the buses of the two systems, according to
Proposition 2, by computing the parametrized state-deformation matrix; see,
Proposition 1. The corresponding LEs are computed and then ordered to obtain
the stability index of each bus in the system denoted as $\mathrm{S}_{i}$ for
$i\in\mathcal{N}$. The right column of Figure 2 depicts the ordered LEs that
are equivalent to the stability index as discussed in Section IV-A. For
$\mathrm{case}$-$\mathrm{9}$, Bus No. $1,2,3$ are generator buses; their
stability index depends on frequency, voltage or rotor angle. For Bus. No.$3$
stability index is based on frequency stability and that of Bus No. $2$ and
$3$ is computed based on voltage stability. This is a consequence of Property
1 that is discussed in Remark 4. Similarly for $\mathrm{case}$-$\mathrm{200}$,
the stability quantification of the buses can result from instability related
to frequency, voltage or rotor angle as can be seen from Figure 2. These
answers question $\mathrm{Q}1$ posed at the beginning of this section.

(a)  (b)

Figure 3: Nodes from most to least stable for networks
$\mathrm{case}$-$\mathrm{9}$ (left) and $\mathrm{case}$-$\mathrm{200}$
(right). The square nodes represent generator buses while the circle ones are
loads/renewables.

Based on the validity of stability quantification method resulting from the
computation of the system’s spectrum of LEs, we now solve P1. The optimal
uncertain renewable load allocation problem is solved according to Algorithm 1
as discussed in Section V-A. The ranked nodes (buses) denoted by
$\mathcal{S}^{*}$ are depicted in Figure 3. The buses are ranked from most
stable to least stable, thereby indicating which nodes allow for allocating a
renewable load while having the least impact on the overall stability of the
power network. For $\mathrm{case}$-$\mathrm{9}$ the optimal bus location is a
generator Bus No.$2$ this is followed by a load Bus No.$4$ and a generator Bus
No.$1$. Whereas by referring to Figure 2, notice that the top three stable
nodes in order are $\\{1,2,4\\}$. This is due to the following: when computing
the stability index, we quantify the stability of a node based on the LEs
resulting from frequency, voltage, or rotor angle instability. However, for
the allocation problem when solving
$\mathrm{log\,det}\left(\tilde{\bm{\Xi}}(\mathcal{S},\bm{\gamma})\right)$, we
sum all the LEs of all the buses within the network. That is, we are
quantifying the impact of a renewable load injection on the stability of all
the nodes by computing all the LEs of the system after the uncertain load is
added to the power grid. As such, the stability index of the buses can be
different than the stable nodes identified by solving problem P1. The results
for $\mathrm{case}$-$\mathrm{200}$ are also depicted in Figure 3. For brevity,
we have included the ranking of the top $30$ stable nodes. Notice that, mostly
the stable buses are load buses (buses numbered $50$ to $200$). Some generator
buses are indicative of stability towards RER injections such as Buses No.
$\\{7,2,14\\}$.

Now, in order to validate the optimality of the ranked nodes $\mathcal{S}^{*}$
for each of the test systems, we simulate the dynamics by allocating an
uncertain RER at the most stable node and compare to a system dynamics
scenario that has a random RER allocation. The frequency response of the
generators for $\mathrm{case}$-$\mathrm{9}$ and $\mathrm{case}$-$\mathrm{200}$
are shown in Figure 4. Notice that the transients clear faster for the case
when the same renewable load injections are applied to the optimal stable node
(right column) for both test systems. The clearing time for the transients for
$\mathrm{case}$-$\mathrm{9}$ is $0.5$ seconds faster and has a damped
transient as compared to the random node. For $\mathrm{case}$-$\mathrm{200}$,
the transient response is also damped and has around a one second shorter
response. The resulting faster clearing time indicates that uncertainty from
renewable load injections can propagate slower depending on where the RER is
allocated. This also demonstrates the optimality of the stability result
obtained from solving P1, thereby answering questions $\mathrm{Q}2$ and
$\mathrm{Q}3$.

Finally, we numerically verify Theorem 1. In doing so, we evaluate the
$\mathrm{log\,det}$ of matrix $\tilde{\bm{\Xi}}_{i}(\bm{\gamma})$ and compare
it to the sum of the spectrum of LEs computed using (13). For
$\mathrm{case}$-$\mathrm{9}$ and $\mathrm{case}$-$\mathrm{200}$ we obtain an
equivalence relation as presented in Theorem 1. For
$\mathrm{case}$-$\mathrm{9}$ the sum of all LEs is equal to $0.2349$ and for
$\mathrm{case}$-$\mathrm{200}$ the sum of LEs is $2.6334$. The results thus
provide a clear relation between the stability measure used in the optimal
renewable load allocation problem P1 and LEs exponents that quantify system
stability. The equivalence therefore answers question $\mathrm{Q}4$.

(a)  (b)  (c)  (d)

Figure 4: Clearing time for frequency transients induced by allocating an
uncertain renewable load injection at a random bus (left) and at the most
stable bus computed by solving P1.

##  VI Conclusions, limitations and future work

The paper presents a framework for quantifying the stability of a power
network after an uncertain renewable load injection is allocated at a bus in
the power grid. The method is based on computing the spectrum of LEs of a
system; it is based on stability criteria that consider frequency, voltage and
rotor angle stability. The proposed methods allow for quantifying renewables
uncertainty propagation onto the network as a whole. As such, it informs the
operator/utility where to practically allocate renewables while maintaining
the overall stability of the power grid. Future work on this topic will focus
on incorporating renewables-based dynamic models (of solar and wind farms).
This enables modeling the uncertainty from renewable loads and therefore
capturing the full uncertainty propagation within the nodes of the power
system.

## References

  * [1] B. Li and K. Y. Wong, “Optimizing synchronization stability of the Kuramoto model in complex networks and power grids,” Physical Review E, vol. 95, no. 1, pp. 1–8, 2017.
  * [2] H. Bosetti and S. Khan, “Transient Stability in Oscillating Multi-Machine Systems Using Lyapunov Vectors,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 2078–2086, 2018.
  * [3] X. Chen, X. Jin, and Z. Huang, “Data-driven method for identifying the expression of the Lyapunov exponent from random data,” International Journal of Non-Linear Mechanics, vol. 148, no. October 2022, p. 104268, 2023\.
  * [4] Y. J. Isbeih, M. S. E. Moursi, W. Xiao, and E. El-Saadany, “Generator-based threshold for transient stability assessment,” IET Smart Grid, vol. 2, no. 3, pp. 407–419, 2019.
  * [5] C. W. Liu, J. S. Thorp, J. Lu, R. J. Thomas, and H. D. Chiang, “Detection of Transiently Chaotic Swings in Power Systems using Real-Time Phasor Measurements,” IEEE Transactions on Power Systems, vol. 9, no. 3, pp. 1285–1292, 1994.
  * [6] A. Czornik, A. Konyukh, I. Konyukh, M. Niezabitowski, and J. Orwat, “On Lyapunov and Upper Bohl Exponents of Diagonal Discrete Linear Time-Varying Systems,” IEEE Transactions on Automatic Control, vol. 64, no. 12, pp. 5171–5174, 2019.
  * [7] Q. Zong, W. Yao, H. Zhou, H. Zhao, J. Wen, and S. Cheng, “Transient Stability Assessment of Large-Scale Power System Using Predictive Maximal Lyapunov Exponent Approach,” IEEE Transactions on Power Systems, vol. PP, pp. 1–13, 2023.
  * [8] Y. Shen, W. Wu, B. Wang, and S. Sun, “Optimal Allocation of Virtual Inertia and Droop Control for Renewable Energy in Stochastic Look-Ahead Power Dispatch,” IEEE Transactions on Sustainable Energy, vol. 14, no. 3, pp. 1881–1894, 2023.
  * [9] A. R. Sobbouhi and A. Vahedi, “Transient stability prediction of power system; a review on methods, classification and considerations,” Electric Power Systems Research, vol. 190, no. July 2020, p. 106853, 2021.
  * [10] J. Yan, C. C. Liu, and U. Vaidya, “PMU-based monitoring of rotor angle dynamics,” IEEE Transactions on Power Systems, vol. 26, no. 4, pp. 2125–2133, 2011.
  * [11] D. P. Wadduwage, C. Q. Wu, and U. D. Annakkage, “Power system transient stability analysis via the concept of Lyapunov Exponents,” Electric Power Systems Research, vol. 104, pp. 183–192, 2013.
  * [12] S. Dasgupta, M. Paramasivam, U. Vaidya, and V. Ajjarapu, “Real-time monitoring of short-term voltage stability using PMU data,” IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 3702–3711, 2013.
  * [13] S. Wei, M. Yang, J. Qi, J. Wang, S. Ma, and X. Han, “Model-Free MLE Estimation for Online Rotor Angle Stability Assessment with PMU Data,” IEEE Transactions on Power Systems, vol. 33, no. 3, pp. 2463–2476, 2018.
  * [14] L. Escot and J. E. Sandubete, “Estimating Lyapunov exponents on a noisy environment by global and local Jacobian indirect algorithms,” Applied Mathematics and Computation, vol. 436, p. 127498, 2023.
  * [15] S. K. Khaitan and J. D. McCalley, “VANTAGE: A Lyapunov exponents based technique for identification of coherent groups of generators in power systems,” Electric Power Systems Research, vol. 105, pp. 33–38, 2013\.
  * [16] M. Rashidi and E. Farjah, “Lyapunov exponent-based optimal PMU placement approach with application to transient stability assessment,” IET Science, Measurement and Technology, vol. 10, no. 5, pp. 492–497, 2016.
  * [17] N. Safari, C. Y. Chung, and G. C. D. Price, “Novel Multi-Step Short-Term Wind Power Prediction Framework Based on Chaotic Time Series Analysis and Singular Spectrum Analysis,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 590–601, 2017.
  * [18] M. Rashidi and E. Farjah, “LEs based framework for transient instability prediction and mitigation using PMU data,” IET Generation, Transmission and Distribution, vol. 10, no. 14, pp. 3431–3440, 2016.
  * [19] S. Rajeh, Impact of the community structure on the dynamics of complex networks.  PhD thesis, Université Bourgogne Franche-Comté, 2023.
  * [20] K. Yoon, D. Choi, S. H. Lee, and J. W. Park, “Optimal Placement Algorithm of Multiple DGs Based on Model-Free Lyapunov Exponent Estimation,” IEEE Access, vol. 8, pp. 135416–135425, 2020.
  * [21] M. Bazrafshan, N. Gatsis, and E. Dallanese, “Placement and Sizing of Inverter-Based Renewable Systems in Multi-Phase Distribution Networks,” IEEE Transactions on Power Systems, vol. 34, no. 2, pp. 918–930, 2019.
  * [22] P. W. Sauer, M. A. Pai, and J. H. Chow, Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox.  Wiley-IEEE Press, second edi ed., 2017.
  * [23] S. A. Nugroho, A. Taha, N. Gatsis, and J. Zhao, “Observers for Differential Algebraic Equation Models of Power Networks: Jointly Estimating Dynamic and Algebraic States,” IEEE Transactions on Control of Network Systems, vol. 5870, no. c, 2022.
  * [24] B. Guo, O. Karaca, T. Summers, and M. Kamgarpour, “Actuator Placement under Structural Controllability Using Forward and Reverse Greedy Algorithms,” IEEE Transactions on Automatic Control, vol. 66, no. 12, pp. 5845–5860, 2021.
  * [25] V. Kunkel, Peter; Mehrmann, Differential-Algebraic Equations: Analysis and Numerical Solution.  Zurich: European Mathmatical Society, 2006.
  * [26] V. H. Linh and V. Mehrmann, “Lyapunov, Bohl and sacker-sell spectral intervals for differential- algebraic equations,” Journal of Dynamics and Differential Equations, vol. 21, no. 1, pp. 153–194, 2009.
  * [27] S. L. Campbell and C. W. Gear, “The index of general nonlinear DAEs,” Numerische Mathematik, vol. 72, no. 2, pp. 173–196, 1995.
  * [28] Y. Chen and S. Trenn, “The differentiation index of nonlinear differential‐algebraic equations versus the relative degree of nonlinear control systems,” Pamm, vol. 20, no. 1, pp. 2020–2022, 2021.
  * [29] T. Groß, S. Trenn, and A. Wirsen, “Solvability and stability of a power system DAE model,” Systems and Control Letters, vol. 97, pp. 12–17, 2016\.
  * [30] M. L. Crow, Computational methods for electric power systems, third edition.  2015\.
  * [31] K. E. Brenan, S. L. Campbell, S. L. V. Campbell, and L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations.  Society for Industrial and Applied Mathematics, 1996.
  * [32] S. G. Krantz and H. R. Parks, The Implicit Function Theorem: History, Theory, and Applications.  Springer New York, 2013.
  * [33] M. H. Kazma and A. F. Taha, “ODE Transformations of Nonlinear DAE Power Systems,” arXiv: To Appear in IEEE PESGM 2024, no. July, pp. 1–7, 2023\.
  * [34] F. Milano, “Semi-Implicit Formulation of Differential-Algebraic Equations for Transient Stability Analysis,” IEEE Transactions on Power Systems, vol. 31, no. 6, pp. 4534–4543, 2016.
  * [35] C. W. Gear, “Simultaneous Numerical Solution of Differential-Algebraic Equations,” IEEE Transactions on Circuit Theory, vol. 18, no. 1, pp. 89–95, 1971.
  * [36] F. A. Potra, M. Anitescu, B. Gavrea, and J. Trinkle, “A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction,” International Journal for Numerical Methods in Engineering, vol. 66, no. 7, pp. 1079–1124, 2006.
  * [37] F. Milano, ed., Advances in Power System Modelling, Control and Stability Analysis.  London: IET, 2nd editio ed., 2022.
  * [38] M. Leszczyński, P. Perlikowski, T. Burzyński, T. M. Kowalski, and P. Brzeski, “Review of sample-based methods used in an analysis of multistable dynamical systems,” Chaos, vol. 32, no. 8, 2022.
  * [39] Y. Sun and C. Q. Wu, “A radial-basis-function network-based method of estimating Lyapunov exponents from a scalar time series for analyzing nonlinear systems stability,” Nonlinear Dynamics, vol. 70, no. 2, pp. 1689–1708, 2012.
  * [40] B. Hayes and F. Milano, “Viable Computation of the Largest Lyapunov Characteristic Exponent for Power Systems,” Proceedings - 2018 IEEE PES Innovative Smart Grid Technologies Conference Europe, ISGT-Europe 2018, pp. 1–6, 2018.
  * [41] Y. Kawano and J. M. Scherpen, “Empirical differential Gramians for nonlinear model reduction,” Automatica, vol. 127, p. 109534, 2021.
  * [42] M. H. Kazma and A. F. Taha, “Observability for Nonlinear Systems: Connecting Variational Dynamics, Lyapunov Exponents, and Empirical Gramians,” arXiv, 2024.
  * [43] S. A. Nugroho, A. F. Taha, and V. Hoang, “Nonlinear Dynamic Systems Parameterization Using Interval-Based Global Optimization: Computing Lipschitz Constants and Beyond,” IEEE Transactions on Automatic Control, vol. 67, no. 8, pp. 3836–3850, 2022.
  * [44] S. V. Ershov and A. B. Potapov, “On the concept of stationary Lyapunov basis,” Physica D: Nonlinear Phenomena, vol. 118, no. 3-4, pp. 167–198, 1998.
  * [45] S. Balasuriya, “Uncertainty in finite-time lyapunov exponent computations,” Journal of Computational Dynamics, vol. 7, no. 2, pp. 313–337, 2020.
  * [46] F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi, “Characterizing dynamics with covariant lyapunov vectors,” Physical Review Letters, vol. 99, no. 13, pp. 1–4, 2007.
  * [47] L. S. Young, “Mathematical theory of Lyapunov exponents,” Journal of Physics A: Mathematical and Theoretical, vol. 46, no. 25, 2013.
  * [48] L. Barreira, Lyapunov Exponents.  1998\.
  * [49] P. Manneville, “Characterization of Temporal Chaos,” Dissipative Structures and Weak Turbulence, pp. 247–284, 1990.
  * [50] J. Frank and S. Zhuk, “A detectability criterion and data assimilation for nonlinear differential equations,” Nonlinearity, vol. 31, no. 11, pp. 5235–5257, 2018.
  * [51] M. Tranninger, R. Seeber, S. Zhuk, M. Steinberger, and M. Horn, “Detectability Analysis and Observer Design for Linear Time Varying Systems,” IEEE Control Systems Letters, vol. 4, no. 2, pp. 331–336, 2020\.
  * [52] D. Martini, D. Angeli, G. Innocenti, and A. Tesi, “Ruling Out Positive Lyapunov Exponents by Using the Jacobian’s Second Additive Compound Matrix,” IEEE Control Systems Letters, vol. 6, pp. 2924–2928, 2022.
  * [53] L. Dieci and E. S. Van Vleck, “Computation of a few Lyapunov exponents for continuous and discrete dynamical systems,” Applied Numerical Mathematics, vol. 17, no. 3, pp. 275–291, 1995.
  * [54] L. Dieci and C. Elia, “SVD algorithms to approximate spectra of dynamical systems,” Mathematics and Computers in Simulation, vol. 79, no. 4, pp. 1235–1254, 2008.
  * [55] L. Dieci, M. S. Jolly, and E. S. Van Vleck, “Numerical techniques for approximating Lyapunov exponents and their implementation,” Journal of Computational and Nonlinear Dynamics, vol. 6, no. 1, pp. 1–7, 2011.
  * [56] A. Pikovsky and A. Politi, Lyapunov Exponents: A Tool to Explore Complex Dynamics, vol. 70.  2017\.
  * [57] P. Masarati and A. Tamer, “Sensitivity of trajectory stability estimated by Lyapunov characteristic exponents,” Aerospace Science and Technology, vol. 47, pp. 501–510, 2015.
  * [58] M. Balcerzak, A. Dabrowski, B. Blazejczyk-Okolewska, and A. Stefanski, “Determining Lyapunov exponents of non-smooth systems: Perturbation vectors approach,” Mechanical Systems and Signal Processing, vol. 141, p. 106734, 2020.
  * [59] M. H. Kazma, S. A. Nugroho, A. Haber, and A. F. Taha, “State-Robust Observability Measures for Sensor Selection in Nonlinear Dynamic Systems,” 2023 62nd IEEE Conference on Decision and Control (CDC), no. Cdc, pp. 8418–8426, 2023.
  * [60] G. L. Nemhauser, L. A. Wolsey, and M. L. Fisher, “An analysis of approximations for maximizing submodular set functions-I,” Mathematical Programming, vol. 14, no. 1, pp. 265–294, 1978.
  * [61] T. H. Summers, F. L. Cortesi, and J. Lygeros, “On Submodularity and Controllability in Complex Dynamical Networks,” IEEE Transactions on Control of Network Systems, vol. 3, pp. 91–101, mar 2016.

##  Appendix A Power System Model

For a synchronous generator $i\in\mathcal{G}$; its $4$-th order differential
dynamics can be written as

$\displaystyle\dot{\delta_{i}}=\omega_{i}-\omega_{0},$ (27a) $\displaystyle
M_{i}\dot{\omega}_{i}=T_{\mathrm{N}i}-P_{\mathrm{G}i}-D_{i}(\omega_{i}-\omega_{0}),$
(27b) $\displaystyle
T^{{}^{\prime}}_{\mathrm{d0}i}\dot{E^{{}^{\prime}}_{i}}=-\frac{x_{\mathrm{d}i}}{x^{{}^{\prime}}_{\mathrm{d}i}}E^{{}^{\prime}}_{i}+\frac{x_{\mathrm{d}i}-{x_{\mathrm{d}i}^{{}^{\prime}}}}{x^{{}^{\prime}}_{\mathrm{d}i}}v_{i}\cos(\delta_{i}-\theta_{i})+E_{\mathrm{fd}i},$
(27c) $\displaystyle
T_{\mathrm{CH}i}\dot{T}_{\mathrm{N}i}=T_{\mathrm{N}i}-\frac{1}{R_{\mathrm{D}i}}(\omega_{i}-\omega_{0})+T_{\mathrm{r}i},$
(27d)

where the time varying components in (27) are: $\delta_{i}$ the rotor angle
$\mathrm{(rad)}$, $\omega_{i}$ generator rotor speed $\mathrm{(rad/sec)}$,
${E}^{{}^{\prime}}_{i}$ generator transient voltage $\mathrm{(pu)}$,
${T}_{\mathrm{N}i}$ generator mechanical torque $(\mathrm{pu})$. Generator
inputs are: $E_{\mathrm{fd}i}$ generator internal field voltage
$(\mathrm{pu})$, $T_{\mathrm{r}i}$ governor reference signal $(\mathrm{pu})$.
Constants in (27) are: $M_{i}$ is the rotor inertia constant
$(\mathrm{pu}\times\mathrm{sec}^{2})$, $D_{i}$ is the damping coefficient
$(\mathrm{pu}\times\mathrm{sec}^{2})$, $x_{\mathrm{d}i}$ and $x_{\mathrm{q}i}$
are the direct-axis synchronous reactance $(\mathrm{pu})$,
$x^{{}^{\prime}}_{\mathrm{d}i}$ is the direct-axis transient reactance
$(\mathrm{pu})$, $T^{{}^{\prime}}_{\mathrm{d0}i}$ is the direct-axis open-
circuit time constant $\mathrm{(sec)}$, $T_{\mathrm{CH}i}$ is the chest valve
time constant $\mathrm{(sec)}$, $R_{\mathrm{D}i}$ is the speed governor
regulation constant $\mathrm{(Hz/pu)}$, and $\omega_{0}$ is the synchronous
speed $(120\pi\;\mathrm{rad/sec})$.

The algebraic constraints of the power system represent the relation between
the internal states of a synchronous generator, and it’s generated power
$P_{\mathrm{G}i}$ and $Q_{\mathrm{G}i}$, i.e., real and reactive power. The
algebraic constraints of the NL-DAE system can be written as (28) with
$i\in\mathcal{G}$

$\displaystyle P_{\mathrm{G}i}=$
$\displaystyle\tfrac{1}{x^{{}^{\prime}}_{\mathrm{d}i}}E^{{}^{\prime}}_{i}v_{i}\sin(\delta_{i}-\theta_{i})-\tfrac{x_{\mathrm{q}i}-x^{{}^{\prime}}_{\mathrm{d}i}}{2x^{{}^{\prime}}_{\mathrm{d}i}x_{\mathrm{q}i}}v^{2}_{i}\sin(2(\delta_{i}-\theta_{i}))$
(28a)
$\displaystyle\begin{split}Q_{\mathrm{G}i}=&\tfrac{1}{x^{{}^{\prime}}_{\mathrm{d}i}}E^{{}^{\prime}}_{i}v_{i}\cos(\delta_{i}-\theta_{i})-\tfrac{x_{\mathrm{q}i}-x^{{}^{\prime}}_{\mathrm{d}i}}{2x^{{}^{\prime}}_{\mathrm{d}i}x_{\mathrm{q}i}}v^{2}_{i}\\\
&-\tfrac{x_{\mathrm{q}i}-x^{{}^{\prime}}_{\mathrm{d}i}}{2x^{{}^{\prime}}_{\mathrm{d}i}x_{\mathrm{q}i}}v^{2}_{i}\cos(2(\delta_{i}-\theta_{i})).\end{split}$
(28b)

where $\theta_{ij}=\theta_{i}-\theta_{j}$ is the bus angle, $v_{i}$ is the bus
voltage $\mathrm{(pu)}$.

The power balance between the set of generator and load buses with
$i\in\mathcal{G}\cup\mathcal{L}$ can be written as (29) such that,
$N:=|\mathcal{N}|$ is the number of buses within the transmission network
while, $G:=|\mathcal{G}|$, $L:=|\mathcal{L}|$ and $R:=|\mathcal{R}|$ are the
number of generator, load and renewable buses. The power balance in (29)
resembles the power transfer between RER, generators and loads as follows

$\displaystyle P_{\mathrm{G}i}+P_{\mathrm{L}i}+P_{\mathrm{R}i}$
$\displaystyle=\sum_{j=1}^{N}v_{i}v_{j}(G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij}),$
(29a) $\displaystyle Q_{\mathrm{G}i}+Q_{\mathrm{L}i}+Q_{\mathrm{R}i}$
$\displaystyle=\sum_{j=1}^{N}v_{i}v_{j}(G_{ij}\cos\theta_{ij}-B_{ij}\sin\theta_{ij}),$
(29b)

where matrices $G_{ij}$ and $B_{ij}$ denote, respectively, the conductance and
susceptance between bus $i$ and $j$.

##  Appendix B Newton-Raphson algorithm

The NR iterative method ensures state convergence for each time-step $k$ by
evaluating the Jacobian of the nonlinear dynamics for each $k$ under a NR
iteration $i$. The Jacobian is used to compute increment
$\Delta\bm{x}_{k}^{(i)}$. The computed increment is then used to update the
state vector to the next time-step as
$\bm{x}^{(i+1)}_{k}=\bm{x}_{k}^{(i)}+\Delta\bm{x}_{k}^{(i)}$ for each NR
iteration $i$. Once a convergence criterion is satisfied,
$||\Delta\bm{x}_{k}^{(i)}||_{2}\leq\varepsilon$, time-step $k$ advances to
$k+1$ and iteratively proceeds for the whole time-span $t$. Here $\varepsilon$
is a small and positive convergence criterion. With that in mind, the
iteration NR increment $\Delta\bm{x}_{k}^{(i)}$ can be written as

$\leavevmode\nobreak\
\Delta\bm{x}^{(i)}_{k}=\left[\bm{A}_{g}(\bm{z}^{(i)}_{k})\right]^{-1}\begin{bmatrix}\bm{\varphi}(\bm{z}^{(i)}_{k})\end{bmatrix},$
(30)

where
$\bm{\varphi}(\bm{z}^{(i)}_{k},\bm{z}_{k-1}):=\bm{\varphi}(\bm{z}^{(i)}_{k})$
is the discrete-time nonlinear model (3) in implicit form such that
$\bm{\varphi}(\bm{z}^{(i)}_{k})=0\;\forall\;k$. The Jacobian
$\bm{A}_{g}(\bm{z}^{(i)}_{k}):=\begin{bmatrix}\tfrac{\partial\bm{\varphi}(\bm{z}^{(i)}_{k})}{\partial\bm{x}_{k}}\end{bmatrix}\in\mathbb{R}^{n\times
n}$ is defined as

$\leavevmode\nobreak\
\bm{A}_{g}=\begin{bmatrix}\bm{I}_{n_{d}}-\tilde{h}\bm{F}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1})&-\tilde{h}\bm{F}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1})\\\
-\tilde{h}\bm{G}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1})&\bm{I}_{n_{a}}-\tilde{h}\bm{G}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1})\end{bmatrix},$
(31)

where matrices
$\bm{F}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1}):=\bm{F}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k})+\bm{F}_{\bm{x}_{d}}(\bm{z}_{k-1})\in\mathbb{R}^{n_{d}\times
n_{d}}$ and
$\bm{F}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1}):=\bm{F}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k})+\bm{F}_{\bm{x}_{a}}(\bm{z}_{k-1})$
$\in\mathbb{R}^{n_{d}\times n_{a}}$ are the Jacobians of (4) with respect to
$\bm{x}_{d}$ and $\bm{x}_{a}$. Matrix $\bm{I}_{n_{d}}$ is an identity matrix
of dimension similar to $\bm{F}_{\bm{x}_{d}}(\cdot)$. Matrices
$\bm{G}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1}):=\bm{G}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k})+\bm{G}_{\bm{x}_{d}}(\bm{z}_{k-1})\in\mathbb{R}^{n_{a}\times
n_{d}}$ and
$\bm{G}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k},\bm{z}_{k-1}):=\bm{G}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k})+\bm{G}_{\bm{x}_{a}}(\bm{z}_{k-1})\in\mathbb{R}^{n_{a}\times
n_{a}}$ are the Jacobians with respect to $\bm{x}_{d}$ and $\bm{x}_{a}$.
Matrix $\bm{I}_{n_{a}}$ is an identity matrix of dimension similar to
$\bm{G}_{\bm{x}_{a}}(\cdot)$. The NR iterative method is summarized in
Algorithm 2.

1 Input: $\bm{x}_{0}$, $\epsilon$, $\max$ iteration

2 Output:
$\bm{x}^{*}=\bm{x}_{k+1}\;\forall\;k\;\in\;\\{1,2,\ldots,\mathrm{N}-1\\}$

3 Initialize: $i\leftarrow 0$, $\bm{x}^{(i)}\leftarrow\bm{x}_{0}$

4 forall  _$k\;\in\;\\{1,2,\ldots,\mathrm{N}-1\\}$_ do

5 while  _$||\Delta\bm{x}_{k}^{(i)}||_{2}$ $>\epsilon$ and $i<\max$ iteration_
do

6 compute: Jacobian $\bm{A}_{g}(\bm{z}^{(i)})$ as in (31)

7 solve for: $\Delta\bm{x}^{(i)}$ as in (30)

8 update: $\bm{x}^{(i+1)}\leftarrow\bm{x}^{(i)}+\Delta\bm{x}^{(i)}$

9 update: error $\leftarrow\|\Delta\bm{x}^{(i)}\|_{2}$

10 if  _error $\leq\epsilon$_ then

11 return $\bm{x}^{(i+1)}=x^{*}$

12

13 else

14 $i\leftarrow i+1$

15

16

17

Algorithm 2 Newton-Raphson Algorithm

##  Appendix C Discrete-QR Method

The discrete QR-factorization for the variational transition matrix (11) can
be obtained according to the methods developed in [54]. For brevity, we
summarize the implementation of the discrete-QR method in the following
algorithm.

1 Input: Initial factorization $\bm{\Phi}^{k}_{0}=\bm{Q}_{k}\bm{R}^{k}_{0}$

2 Output: QR factorization $\;\forall\;k\;\in\;\\{1,2,\ldots,\mathrm{N}-1\\}$

3 Initialize: $\bm{Q}_{0}$, $\bm{R}^{0}_{0}$

4 forall  _$k\;\in\;\\{1,2,\ldots,\mathrm{N}-1\\}$_ do

5 compute: $\bm{\Phi}_{0}^{k}(\bm{x}_{0})$
$\bm{\Phi}_{0}^{k}=\left(\bm{I}_{n}+\tfrac{\partial\tilde{\bm{f}}(\bm{x}_{k},\bm{x}_{k-1})}{\partial\bm{x}_{k}}\right)\tfrac{\partial\bm{x}_{k}}{\partial\bm{x}_{0}},\quad\bm{\Phi}_{0}^{0}(\bm{x}_{0})=\bm{I}_{n}$

6 compute QR factorization:
$\bm{\Phi}^{k}_{0}\bm{Q}_{k-1}=\bm{Q}_{k}\bm{R}^{k}_{k-1}$

7 record and update: $\bm{Q}_{k-1}\leftarrow\bm{Q}_{k}$

8 record and update:$\bm{R}^{k-1}_{k-2}\leftarrow\bm{R}^{k}_{k-1}$

9

Algorithm 3 Discrete-QR Factorization Algorithm

##  Appendix D Submodular set function maximization

For any set function, denoted by $\mathcal{L}(\cdot)$, that is submodular and
monotone increasing then the maximization of the set function computed using
the greedy algorithm offers a theoretical worst-case bound according to the
following theorem.

######  Theorem 2.

([60]) Let $\mathcal{L}:2^{V}\rightarrow\mathbb{R}$ be a submodular and
monotone increasing set function, $\mathcal{L}^{*}$ be the optimal solution of
P1 and $\mathcal{L}^{*}_{\mathcal{S}}$ be the solution computed using the
greedy algorithm. Then, the following performance bound holds true

$\displaystyle\mathcal{L}^{*}_{\mathcal{S}}-\mathcal{L}(\emptyset)\geq\left(1-\frac{1}{e}\right)\left(\mathcal{L}^{*}-f(\emptyset)\right),\quad\text{with}\;\;\;\mathcal{L}(\emptyset)=0,$

where $e\approx 2.71828$.

We note that the above bound is theoretical, and generally a greedy approach
performs better in practice. In practice when considering a submodular set
maximization problem, it has been shown the performance accuracy achieves a
$99\%$ guarantee; see [61] and the many references that cite this work.

##  Appendix E Proof of Proposition 3

The following is the proof for Proposition 3.

###### Proof.

Let $\mathcal{L}_{s}:2^{V-\\{s\\}}\rightarrow\mathbb{R}$ denote a derived set
function is defined as follows

$\displaystyle\mathcal{L}_{s}(\mathcal{S})$
$\displaystyle=\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}\left({\mathcal{S}\cup\\{s\\}}\right)\right)-\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}(\mathcal{S})\right),$
$\displaystyle=\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}(\mathcal{S})+\tilde{\bm{\Xi}}(\\{s\\})\right)-\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}(\mathcal{S})\right).$

We first show $\mathcal{L}(\mathcal{S})$ that is monotone decreasing for any
$s\in V$. Let $\mathcal{A}\subseteq\mathcal{B}\subseteq\mathcal{V}-\\{s\\}$,
and let
$\tilde{\bm{\Xi}}(\bm{\alpha})=\tilde{\bm{\Xi}}(\mathcal{A})+\bm{\alpha}\left(\tilde{\bm{\Xi}}(\mathcal{B})-\tilde{\bm{\Xi}}(\mathcal{A})\right)$
for $\bm{\alpha}\in[0,1]$. Then for

$\displaystyle\tilde{\mathcal{L}_{s}}(\tilde{\bm{\Xi}}(\bm{\alpha}))=\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}(\bm{\alpha})+\tilde{\bm{\Xi}}(\mathcal{S})\right)-\mathrm{log}\,\mathrm{det}\left(\tilde{\bm{\Xi}}(\bm{\alpha})\right),$

we obtain the following

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\bm{\alpha}}\tilde{\mathcal{L}_{s}}(\tilde{\bm{\Xi}}(\bm{\alpha}))$
$\displaystyle=\mathrm{trace}\Big{[}\left(\left(\tilde{\bm{\Xi}}(\bm{\alpha})+\tilde{\bm{\Xi}}(\mathcal{S})\right)^{-1}-\tilde{\bm{\Xi}}(\bm{\alpha})^{-1}\right)\left(\tilde{\bm{\Xi}}(\mathcal{B})-\tilde{\bm{\Xi}}(\mathcal{A})\right)\Big{]}$
$\displaystyle\leq 0.$

Such that

$\displaystyle\left(\left(\tilde{\bm{\Xi}}(\bm{\alpha})+\tilde{\bm{\Xi}}(\mathcal{S})\right)^{-1}-\;\;\tilde{\bm{\Xi}}(\bm{\alpha})^{-1}\right)^{-1}\preceq
0,$

and

$\displaystyle\left(\tilde{\bm{\Xi}}(\mathcal{B})-\tilde{\bm{\Xi}}(\mathcal{A})\right)\succeq
0,$

then the above inequality holds. Thus, we have $\mathcal{L}_{s}$ is monotone
decreasing, and $\mathcal{L}(\mathcal{S})$ is submodular. On such note, the
above proposition holds true. ∎