#  ODE Transformations of Nonlinear DAE Power Systems*

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>
*The discrepancy between this version of the paper and the PES General Meeting 2024 version (available on IEEE Explore) is due to the imposed page limit for all conference papers. This version includes additional comments as well as an easier paper layout for better readability.

###### Abstract

Dynamic power system models are instrumental in real-time stability,
monitoring, and control. Such models are traditionally posed as systems of
nonlinear differential algebraic equations (DAEs): the dynamical part models
generator transients and the algebraic one captures network power flow. While
the literature on control and monitoring for ordinary differential equation
(ODE) models of power systems is indeed rich, that on DAE systems is not. DAE
system theory is less understood in the context of power system dynamics. To
that end, this paper presents two new mathematical transformations for
nonlinear DAE models that yield nonlinear ODE models while retaining the
complete nonlinear DAE structure and algebraic variables. Such transformations
make (more accurate) power system DAE models more amenable to a host of
control and state estimation algorithms designed for ODE dynamical systems. We
showcase that the proposed models are effective, simple, and computationally
scalable.

###### Index Terms:

Time-domain simulation, transient stability analysis, nonlinear descriptor
models, power systems.

##  I Introduction

Power systems monitoring, state estimation, control, and transient stability
analysis are all reliant on high-fidelity models of multi-machine power
systems. In power grids, transient stability analysis determines how the power
system maintains synchronicity under time-varying conditions and large
uncertainties from load disturbances [1, 2]. Such analysis relies on time-
domain simulations of the system that is expressed as a set of differential
algebraic equations (DAEs) [3, 4]. The differential-algebraic nature couples
the system dynamics with power flow constraints, thus resulting in a more
accurate model. Nonlinear DAEs are an extreme case of stiff dynamical systems
[5], meaning that the system has time constants that span several orders of
magnitude; in particular, the algebraic constraints exhibit null time
constants.

In general, nonlinear DAEs are solved using implicit discretization schemes
[6]. Multi-step methods offer stable and efficient schemes when dealing with
nonlinear DAEs [7]. Such discrete-time modeling methods include: backward
differential formulas (BDF) [8, 9], backward Euler (BE) method [9], and
trapezoidal implicit (TI) method [10, 9]. Simulating discrete-time models
requires an integrative time-step algorithm [11], and the solvability of power
system DAEs under implicit integration methods is well-established. The
Newton-Raphson (NR) method [12, 9] is generally implemented within power
system simulation packages to solve discretized DAEs [13]. Despite such well-
developed time-domain numerical solutions, from a systems’ theory perspective,
the literature on nonlinear DAE power networks is limited—unlike that of
ordinary differential equations (ODE models) [14].

Common modeling for systems’ control and estimation is based on an ODE
formulation; this is due to the aforementioned limitation on system theory for
DAEs. Typically, the formulation of ODE systems from DAE models is performed
by either neglecting the algebraic constraints or by formulating a decoupled
modeling approach [15]. The simplified models potentially limit the transient
stability simulations and, ultimately, the estimation and control performance.
Time-domain simulations resulting from the full DAE models of a power network
can, for instance, give an accurate depiction of the dynamics under
topological changes triggered by faults and be modeled to include uncertain
loads from renewable energy resources.

The limitation on model fidelity from a control and estimation perspective is
expressed in the form of the following research question: How do we extend
existing systems control theory, developed for ODE dynamical systems, to
accurately apply it to DAE models? Descriptor systems—arising from DAEs
models—appear in numerous applications, with a few examples being chemical,
electrical and mechanical systems. For such reason, there is a rise in
interest towards translating control and stability theory to the analysis of
descriptor systems [16].

Recent studies—see, [17, 18, 19, 16, 15, 20] and reference therein—present
literature on developing the state estimation and control theory of DAE
systems, in particular that of linear DAEs. However, in this paper, we aim to
address the limitations on state estimation, control, and transient stability
analysis of power systems by giving a new perspective on DAE to ODE system
modeling. Therefore, we instead attempt to address the posed research question
in the form of: Is there a methodology to accurately restructure the DAE
system into an ODE model without loss of information and therefore exploit
existing ODE systems control theoretic?

To that end, in this paper we introduce two simple yet effective methods to
transform nonlinear DAE power system models into ODEs. The idea is to
formulate ODE-structured representations of the network dynamics while
depicting the full nonlinear DAE structure along with the algebraic
constraints. These transformations allow the utilization of the rich
literature on control and estimation of ODE models of power systems.The main
contributions of this paper are:

  * •

We present two mathematical transformations that retain the complete nonlinear
DAE dynamics while achieving a nonlinear ODE-structure. The first
transformation is based on applying the implicit function theorem (IFT) to
reformulate the algebraic constraints into an ODE model (Section II-A). The
second transformation is approximation-based; it yields an effective
approximation towards ODE-structured algebraic constraints (Section II-B).

  * •

We show that it is viable to model and perform both continuous and discrete
time-domain simulations for transient stability analysis on the proposed
systems. We also illustrate that the resulting Hessian matrices computations
that arise when simulating the discretized IFT transformed system can be
approximated using finite difference approximations (Section III).

  * •

We assess the validity of the proposed models on standard power networks under
transient time-domain simulations subject to load disturbances (Section IV).

We note that, in this paper, we do not focus on extending DAE systems
stability and control theory, or on formulating new discretization techniques
for such models, but rather on introducing the aforementioned transformations
onto the nonlinear DAE power system. The transformations’ significance is
accentuated by the fact that we obtain nonlinear ODE structured models that
retain the complete nonlinear DAE structure along with the algebraic
variables. Utilizing the ODE structured models, the rich literature on control
and monitoring of ODE dynamical systems can be exploited.

##  II Proposed Nonlinear DAE Transformations

In this paper, we consider a general nonlinear model of multi-machine power
system dynamics $(\mathcal{N},\mathcal{E})$, where
$\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of transmission
lines, $\mathcal{N}=\mathcal{G}\cup\mathcal{L}$ is the set of all buses in the
network, while $\mathcal{G}$ and $\mathcal{L}$ are the set of generator and
load buses, respectively. The model represents both the generator dynamics and
algebraic constraints. Readers can refer to [13, Ch. 7] for the full
description of the power network utilized within this work. Note that the
theoretical developments herein still apply to any semi-explicit nonlinear DAE
model of a power system written as

$\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 generator are defined as
$\bm{x}_{d}:=\bm{x}_{d}(t)\in\mathbb{R}^{n_{d}}$, the algebraic states as
$\bm{x}_{a}:=\bm{x}_{a}(t)\in\mathbb{R}^{n_{a}}$ and the system input as
$\bm{u}:=\bm{u}(t)\in\mathbb{R}^{n_{u}}$. Functions
$\bm{f}(\cdot):\mathbb{R}^{n_{d}}\times\mathbb{R}^{n_{a}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{d}}$
and
$\bm{g}(\cdot):\mathbb{R}^{n_{d}}\times\mathbb{R}^{n_{a}}\rightarrow\mathbb{R}^{n_{a}}$
are nonlinear and define the system dynamics and power flows.

The existence of a solution for nonlinear DAEs can be determined by proving
that a DAE is strangeness-free [21, Hypothesis 4.2], i.e., the strangeness
index is equal to zero. This index is a generalization of the differentiation
index of DAEs. Refer to [21] for the detailed hypothesis that defines the
strangeness index of nonlinear DAEs.

######  Definition 1.

The differentiation index [22, 21, 23] of descriptor system refers to the
number of differentiations required to obtain ODEs using algebraic
manipulations.

For linear DAE systems, a set of differential algebraic equations is of index
one if and only if it is regular. Regularity is an important property for
linearized DAEs; it is a condition for the existence of a consistent unique
solution for every initial condition [14].

######  Definition 2.

Regularity of a linearized DAE (1) around an initial state can be
characterized by matrix pair $(\bm{E},\bm{A})$, such that it is regular if and
only if $\mathrm{det}(s\bm{E}-\bm{A})\neq 0$ for $s\in\mathbb{C}$.

The linearized representation around an operating point of the power system
DAE (1) can be rewritten as

$\leavevmode\nobreak\ \bm{E}\dot{\bm{x}}=\bm{A}\bm{x}+\bm{B}\bm{u},$ (2)

where $\bm{E}\in\mathbb{R}^{n_{d}+n_{a}}$ represents the singular mass matrix
that has ones on its diagonal entry for the differential equations and zeros
for the algebraic equations. The constant state-space matrices are defined as
$\bm{A}\in\mathbb{R}^{n_{d}+n_{a}}$ and $\bm{B}\in\mathbb{R}^{n_{u}}$.

In [15], it is shown that the linearized representation of power system (1)
has a differentiation index of one and is regular. Such condition guarantees
that for each consistent initial condition, a unique solution exists. Under
the linearized model, this condition holds true if and only if there exist
paths whereby every load is connected to a generator bus. To prove a solution
is unique for nonlinear DAEs is rather complex and still considered an open
problem [14]. We note here that proving the linearized dynamic system (2) is
regular and of index one is a prerequisite for providing evidence that the
nonlinear system can be strangeness-free. The differentiation index,
regularity and strangeness-free property of a nonlinear DAE system ensure the
solvability of the system. Considering the aforementioned results from [15]
and the complexity of proving such condition for nonlinear DAEs, the following
assumption holds true within this paper.

######  Assumption 1.

The DAE (1) is strangeness-free, of differentiation index one and regular.
Thus, under any consistent initial conditions a unique solution exists and the
partial derivatives of $\bm{f}(\cdot)$ and $\bm{g}(\cdot)$ with respect to
$\bm{x}_{d}$ and $\bm{x}_{a}$ are non-singular [21].

The aforementioned assumption is mild and holds true for the test cases
considered in the numerical studies section. Under such conditions, the two
transformations that reformulate (1) into nonlinear ODEs can be posed.

######  Remark 1.

The time-domain numerical solvability of the DAE system (1) implies that a
unique solution for different inputs and consistent initial conditions exists
[24].

Hence, the time-domain simulations presented in Section IV validate Assumption
1 for the power system represented as (1) and therefore provide explicit proof
that the nonlinear power system DAE is strangeness-free, of differentiation
index one and regular.

###  II-A IFT-Based Nonlinear DAE Model: $\mathrm{ODE}$-$\mathrm{DAE}$

The first transformation relies on applying the IFT [25, Theorem 3.3.1] to
resolve the algebraic constraints into ODEs $(\mathrm{ODE}$-$\mathrm{DAE})$.
This method entails differentiating the algebraic constraints in (1b) with
respect to time variable $t$.

######  Proposition 1.

Consider the nonlinear descriptor system (1) and that Assumption 1 holds true
for any initial condition and control input, then the descriptor system can be
restructured into a nonlinear ODE system by applying the IFT [25, Theorem
3.3.1] to resolve the algebraic constraints into ODEs. As such, the nonlinear
descriptor 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=-(\bm{G}_{\bm{x}_{a}})^{-1}\bm{G}_{\bm{x}_{d}}\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u})=\tilde{\bm{g}}(\bm{x}_{d},\bm{x}_{a},\bm{u}).\vspace{-0.1cm}$
(3b)

###### Proof.

Under Assumption 1, the partial derivative
$\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{a}}$ is non-
singular, refer to [14, 15]. Then, implicitly differentiating the algebraic
constraints (1b) using the implicit function theorem we obtain the following

$\leavevmode\nobreak\
0=\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{d}}\frac{\partial\bm{x}_{d}}{\partial
t}+\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{a}}\frac{\partial\bm{x}_{a}}{\partial
t}.$ (4)

We define the Jacobian matrices of the implicit algebraic constraints with
respect to $\bm{x}_{a}$ and $\bm{x}_{d}$ as
$\bm{G}_{\bm{x}_{a}}:=\frac{\partial\bm{g}(\bm{x}_{d},\bm{x}_{a})}{\partial\bm{x}_{a}}$
$\in\mathbb{R}^{n_{d}\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}}$. The time-derivative
$\frac{\partial\bm{x}_{d}}{\partial
t}:=\dot{\bm{x}}_{d}=\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u})$, refer to (1a). As
such, the algebraic constraints (1b) can be rewritten as (3b). ∎

The resulting $\mathrm{ODE}$-$\mathrm{DAE}$ (3b) representing the algebraic
constraints (1b) is now a function of (1a) and therefore depends on control
input $\bm{u}$. Nevertheless (1b) can be explicitly formulated and then
differentiated with respect to $t$, however this requires the complex task of
explicitly reconstructing the algebraic equations—which is rather difficult to
perform. In contrast, the IFT method allows us to utilize the semi-implicit
system (1) along with all the dynamic and algebraic relationships that are
inherent to the system while also being represented as a set of ODEs.

It is important to note that this index reduction approach can be applied to
nonlinear DAEs of differentiation index $n$. In other words, this theoretical
method can also be applied for higher index DAEs, however it requires several
rounds of differentiation to reformulate the algebraic constraints into a set
of ODE equations. The application on DAE of index greater than one is outside
the scope of this paper and power systems in general—power systems are
typically of index one.

###  II-B Approximate Nonlinear DAE Model: $\mathrm{Approx}$-$\mathrm{DAE}$

The second transformation formulates an approximate DAE model
$(\mathrm{Approx}$-$\mathrm{DAE})$ that is based on introducing a positive
scalar term, denoted by $\mu$, to the DAE system at the algebraic constraint
equations level. In particular, the left-hand side of (1b) is replaced by
$\mu\dot{\bm{x}}_{a}$. The rationale behind utilizing such simplistic
alternative as compared to the first transformation is evident at the
discrete-time modeling level; it offers an alternative to dealing with the
Hessian matrix computations that arise under the $\mathrm{ODE}$-$\mathrm{DAE}$
approach—refer to section III. Under such transformation, $\mu>0$ is defined
as a relatively small term that simulates the system’s dynamics while
satisfying the power flow constraints. As a result, the descriptor system (1)
can be rewritten as

$\displaystyle\dot{\bm{x}}_{d}$
$\displaystyle=\bm{f}(\bm{x}_{d},\bm{x}_{a},\bm{u}),$ (5a)
$\displaystyle\mu\dot{\bm{x}}_{a}$
$\displaystyle=\bm{g}(\bm{x}_{d},\bm{x}_{a})+\bm{O}(\mu),$ (5b)

where the approximation error $\bm{O}(\mu)$ is of order $\mu$, such that as
$\mu\rightarrow 0$ the error becomes null.

######  Remark 2.

The proposed $\mathrm{Approx}$-$\mathrm{DAE}$ model is not the model reduction
technique commonly referred to as the singular perturbation technique [26, 27,
28, 11]. Singular perturbation is applied to dynamical systems that have
multiple time-scale dynamics, including electrical systems, to reduce the
number of simulated states; it results in a simplified model representation.

TABLE I: Discretized Dynamics of the Algebraic Constraints $\mathrm{DAE}$ | $\bm{0}=\begin{cases}-\tilde{h}\bm{g}(\bm{x}_{k})&\text{for BDF}\\\ -\tilde{h}(\bm{g}(\bm{x}_{k})+\bm{g}(\bm{x}_{k-1}))&\text{for TI}\end{cases}$
---|---
$\mathrm{ODE}$-$\mathrm{DAE}$ | $\bm{0}=\begin{cases}\bm{x}_{a,k}-\Sigma_{s=1}^{k_{g}}\alpha_{s}\bm{x}_{a,k-s}-\tilde{h}\tilde{\bm{g}}(\bm{z}_{k})&\text{for BDF}\\\ \bm{x}_{a,k}-\bm{x}_{a,k-1}-\tilde{h}(\tilde{\bm{g}}(\bm{z}_{k})+\tilde{\bm{g}}(\bm{z}_{k-1}))&\text{for TI}\end{cases}$
$\mathrm{Approx}$-$\mathrm{DAE}$ | $\bm{0}=\begin{cases}\mu\bm{x}_{a,k}-\Sigma_{s=1}^{k_{g}}\alpha_{s}\mu\bm{x}_{a,k-s}-\tilde{h}\bm{g}(\bm{x}_{k})&\text{for BDF}\\\ \mu\bm{x}_{a,k}-\mu\bm{x}_{a,k-1}-\tilde{h}(\bm{g}(\bm{x}_{k})+\bm{g}(\bm{x}_{k-1}))&\text{for TI}\end{cases}$

Herein, we utilize the small positive scalar $\mu$ to reduce to stiffness of
the DAE model, i.e, transforming the null time-scale constants for the
existing algebraic constraints to become of order $\mu$. The introduced
dynamic $\mu\dot{\bm{x}}_{a}$ is an arbitrary modification to the system’s
behavior and depends highly on the choice of $\mu$. The value of $\mu$
contributes to the stiffness of the system, i.e., the order of $\mu$ defines
the time-scale for the algebraic equations that are modeled as dynamic
equations. Such that as $\mu\rightarrow 0$, the time-scale becomes null, and
therefore the system reverts to having algebraic equations that are extremely
stiff.

##  III Time-domain Modeling for Transient Stability Analysis

The choice of discretization method must rely upon the system’s stiffness,
desired accuracy, and the performance of computation resources. Stiff
dynamical systems can be identified from time constants on local subsystems
that have contrasting magnitudes by a large margin. Nonlinear dynamic power
systems under transient conditions are in practice modeled as discrete-time
state-space models and are solved using numerical methods [24]. Nonlinear DAEs
exhibit stiff dynamics and are solved using implicit discretization methods
that offer stable computational methods as compared with explicit methods [6].

In this work we approach discretizing the nonlinear DAE system and the
proposed transformations using BDF, BE, and TI discrete-time modeling methods.
Specifically, we investigate how the two transformations perform when embedded
within vintage discretization methods. BDF depends on its discretization index
denoted by $k_{g}$, which is stable for $2\leq k_{g}\leq 5$. We note that for
$k_{g}=1$, BDF renders into the BE method. For succinctness, we refer to BE as
BDF while designating the order of $k_{g}=1$. We start by discretizing the
differential dynamics (1a) that are common between the presented models. To
that end, the discrete-time representation of the generator dynamics
represented in (1a) can be written as

$\leavevmode\nobreak\
\bm{0}=\begin{cases}\bm{x}_{d,k}-\Sigma_{s=1}^{k_{g}}\alpha_{s}\bm{x}_{d,k-s}-\tilde{h}\bm{f}(\bm{z}_{k})&\text{for
BDF,}\\\
\bm{x}_{d,k}-\bm{x}_{d,k-1}-\tilde{h}(\bm{f}(\bm{z}_{k})+\bm{f}(\bm{z}_{k-1}))&\text{for
TI,}\end{cases}\\\ $ (6)

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

The discretization constants $\beta$ and $\alpha_{s}$ for BDF method depend on
the order of index $k_{g}$ and are calculated as

$\leavevmode\nobreak\
\beta=\Big{(}\sum_{s=1}^{k_{g}}\frac{1}{s}\Big{)}^{-1}\;,\;\;\;\alpha_{s}=(-1)^{(s-1)}\beta\sum_{j=s}^{k_{g}}\frac{1}{j}\begin{pmatrix}j\\\
s\end{pmatrix}.$ (7)

The implicit discrete-time representation of the algebraic constraints for the
descriptor system (1) and the proposed models are summarized in Table I.
Solving discrete-time models requires an integrative time-step algorithm [11].
The Newton-Raphson (NR) method [12, 9] is generally implemented to solve the
implicit discrete-time nonlinear descriptor dynamics.

TABLE II: Variations in the Jacobian of ODE-Transformed Systems $\mathrm{ODE}$-$\mathrm{DAE}$ | $\begin{bmatrix}-\tilde{h}\widetilde{\bm{G}}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k})&\bm{I}_{n_{a}}-\tilde{h}\widetilde{\bm{G}}_{\bm{x}_{a}}(\bm{z}^{(i)}_{k})\end{bmatrix}$
---|---
$\mathrm{Approx}$-$\mathrm{DAE}$ | $\begin{bmatrix}-\tilde{h}\bm{G}_{\bm{x}_{d}}(\bm{x}^{(i)}_{k})&\mu\bm{I}_{n_{a}}-\tilde{h}\bm{G}_{\bm{x}_{a}}(\bm{x}^{(i)}_{k})\end{bmatrix}$

To implement the NR method, the Jacobian pertaining to the nonlinear dynamics
is evaluated. At each time step $k$, the increment $\Delta\bm{x}_{k}^{(i)}$
defined as (8) is evaluated and used to update the state variables
$\bm{x}^{(i+1)}_{k}=\bm{x}_{k}^{(i)}+\Delta\bm{x}_{k}^{(i)}$ under the NR
iteration $i$ until a convergence criterion is satisfied. Once NR iteration
converges, time step $k$ advances until the dynamics over time span $t$ is
simulated. The iteration increment $\Delta\bm{x}^{(i)}_{k}$ can be written as

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

where
$\bm{z}^{(i)}_{k}:=[\bm{x}^{(i)}_{d,k},\bm{x}^{(i)}_{a,k},\bm{u}^{(i)}_{k}]$,
and $\bm{\phi}(\bm{z}^{(i)}_{k})$ denotes the discretized system dynamics,
represented in (6) and Table I, under the NR iteration index $i$. The Jacobian
of $\bm{\phi}(\bm{z}^{(i)}_{k})$ is defined as
$\bm{A}_{g}(\bm{z}^{(i)}_{k}):=\begin{bmatrix}\tfrac{\partial\bm{\phi}(\bm{z}^{(i)}_{k})}{\partial\bm{x}_{k}}\end{bmatrix}$
and can be written as

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

where matrices $\bm{F}_{\bm{x}_{d}}(\cdot)\in\mathbb{R}^{n_{d}\times n_{d}}$
and $\bm{F}_{\bm{x}_{a}}(\cdot)$ $\in\mathbb{R}^{n_{d}\times n_{a}}$ are the
Jacobians of (6) 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}}(\cdot)$ and
$\bm{G}_{\bm{x}_{a}}(\cdot)$ retain the same definition as before, however
under the discretized dynamics presented in Table I.

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

Figure 1: Transient differential and algebraic state trajectories under load
disturbance: (a) $\mathrm{case}$-$\mathrm{9}$ ($\alpha_{L}=2\%$) and (b)
$\mathrm{case}$-$\mathrm{200}$ ($\alpha_{L}=15\%$).

(a)  (b)

Figure 2: RMSE on dynamic and algebraic state trajectories for the
$\mathrm{Approx}$-$\mathrm{DAE}$ model while changing $\mu$-value.

Regarding the transformed discrete-time models, the Jacobian of the
differential equations remains the same, however, that of the algebraic
constraints summarized in Table I differs. The Jacobian of the algebraic
equations of the transformed systems are represented in Table II. For the
$\mathrm{Approx}$-$\mathrm{DAE}$ system, the only difference is with the
addition of a $\mu\bm{I}_{n_{a}}$ perturbation term, where $\bm{I}_{n_{a}}$ is
an identity matrix of dimension similar to $\bm{G}_{\bm{x}_{a}}(\cdot)$.

(a)  (b)  (c)  (d)  (e)  (f)

Figure 3: Error norm on discrete-time domain state trajectories—under load
disturbance—for the proposed model simulated under different discretization
methods: (a) $\mathrm{case}$-$\mathrm{9}$ ($\alpha_{L}=2\%$) and (b)
$\mathrm{case}$-$\mathrm{200}$ ($\alpha_{L}=15\%$).

For the $\mathrm{ODE}$-$\mathrm{DAE}$ system, matrices
$\widetilde{\bm{G}}_{\bm{x}_{d}}(\cdot):=\frac{\partial\tilde{\bm{g}}(\bm{x}_{d},\bm{x}_{a},\bm{u})}{\partial\bm{x}_{d}}\in\mathbb{R}^{n_{a}\times
n_{d}}$ and
$\widetilde{\bm{G}}_{\bm{x}_{a}}(\cdot):=\frac{\partial\tilde{\bm{g}}(\bm{x}_{d},\bm{x}_{a},\bm{u})}{\partial\bm{x}_{a}}\in\mathbb{R}^{n_{a}\times
n_{a}}$ define the Jacobian of the algebraic equations Table I, row 2 . For
brevity, we define Jacobian $\tilde{\bm{G}}_{\bm{x}_{d}}(\cdot)$ in (10) and
leave $\widetilde{\bm{G}}_{\bm{x}_{a}}(\cdot)$ for the reader to infer.

$\leavevmode\nobreak\
\begin{split}\widetilde{\bm{G}}_{\bm{x}_{d}}(\bm{z}^{(i)}_{k}):=-(&-\bm{G}_{\bm{x}_{a}}^{-1}\tfrac{\partial\bm{G}_{\bm{x}_{a}}}{\partial\bm{x}_{d}}\bm{G}_{\bm{x}_{a}}^{-1}\bm{G}_{\bm{x}_{d}}\bm{f}(\cdot)\\\
&+\bm{G}_{\bm{x}_{a}}^{-1}\tfrac{\partial\bm{G}_{\bm{x}_{d}}}{\partial\bm{x}_{d}}\bm{f}(\cdot)+\bm{G}_{\bm{x}_{a}}^{-1}\bm{G}_{\bm{x}_{d}}\bm{F}_{\bm{x}_{d}}).\end{split}$
(10)

We note here that due to the existence of Jacobian matrices within the state-
space formulation (3b) of the IFT model, Hessian matrices
$\bm{G}_{\bm{x}_{a,d}}:=\tfrac{\partial\bm{G}_{\bm{x}_{a}}}{\partial\bm{x}_{d}}\in\mathbb{R}^{n_{a}\times
n_{d}}$ and
$\bm{G}_{\bm{x}_{d,d}}:=\tfrac{\partial\bm{G}_{\bm{x}_{d}}}{\partial\bm{x}_{d}}\in\mathbb{R}^{n_{a}\times
n_{d}}$ appear under the NR iteration. To avoid the computationally convoluted
evaluation of Hessian matrices, we refer to finite difference (FD)
approximation techniques [29, 30] to approximate the matrices. FD
approximations are based on Taylor series expansions of the system around a
consistent operating point. For the approximations herein, we refer to the
central difference approximation [30] which yields a lower error than the
other difference methods. As such $\bm{G}_{\bm{x}_{a,d}}$ can be approximated
as (11) with difference variable $m\rightarrow 0$.

$\leavevmode\nobreak\
\bm{G}_{\bm{x}_{a,d}}=\tfrac{\bm{G}_{\bm{x}_{a}}(\bm{x}^{(i)}_{d,k-m},\bm{x}^{(i)}_{a,k})-\bm{G}_{\bm{x}_{a}}(\bm{x}^{(i)}_{d,k+m},\bm{x}^{(i)}_{a,k})}{2m}+\bm{O}(m^{2}),$
(11)

where $\bm{G}_{\bm{x}_{a,d}}$ represents the central difference approximation
[30], which yields a truncation error $\bm{O}(m^{2})$.

##  IV Case Studies

The main objective of transient time-domain simulations is to trace the
system’s trajectory—after a disturbance—towards equilibrium [11]. Accordingly,
we investigate the viability and effectiveness of performing both continuous-
and discrete-time transient stability simulations of the proposed systems
following a load disturbance. Two networks of contrasting sizes are chosen:
Western System Coordinating Council (WSCC) 9-Bus network
($\mathrm{case}$-$\mathrm{9}$), and ACTIVSg200-Bus network
($\mathrm{case}$-$\mathrm{200}$).

The transient time-domain simulations 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 baseline model herein refers to the nonlinear DAE (1) and is
simulated using MATLAB ODE/DAE solver $\mathrm{ode15s}$ and $\mathrm{ode15i}$.
The settings chosen for the solvers are: $(i)$ absolute tolerance as $1\times
10^{-06}$, $(ii)$ relative tolerance as $1\times 10^{-05}$ and $(iii)$ maximum
step size equal to $0.001$. For the transient discrete-time domain
simulations, the discretization step size is set to $h=0.1$ and the NR
algorithm parameters are: $(i)$ absolute tolerance on $\mathcal{L}_{2}$–norm
of iteration convergence as $10^{-2}$ and $(ii)$ maximum iterations as $10$.

Generator parameters are extracted from the power system toolbox (PST).
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 steady state initial conditions for
the power system are generated from the power flow solution obtained using
MATPOWER. 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.

Starting from the initial steady state conditions, a load disturbance at $t>0$
on initial load $(\mathrm{P}_{\mathrm{L}}^{0},\mathrm{Q}_{\mathrm{L}}^{0})$ is
introduced. The perturbed magnitude under a load disturbance $(\alpha_{L})$ is
computed as
$(\mathrm{\tilde{P}}_{\mathrm{L}}^{0},\mathrm{\tilde{Q}}_{\mathrm{L}}^{0})=(1+\tfrac{\alpha_{L}}{100})(\mathrm{P}_{\mathrm{L}}^{0},\mathrm{Q}_{\mathrm{L}}^{0})$.
The load disturbance $(\alpha_{L})$ for $\mathrm{case}$-$\mathrm{9}$ is chosen
within the range of $\\{1\%,4\%\\}$ with respect to the original loads; for
test $\mathrm{case}$-$\mathrm{200}$ the range is between $\\{5\%,20\%\\}$.

The continuous-time transient state trajectories representing algebraic state
$\theta$ (bus angle) and differential state $\omega$ (generator synchronous
rotor speed) for the DAE model and the proposed transformations are depicted
in Fig. 1. The proposed models yield accurate transient state trajectory
simulations as compared to the baseline DAE model for both test cases.

Before moving forward, we investigate the choice of $\mu$ on the accuracy of
the $\mathrm{Approx}$-$\mathrm{DAE}$ system as compared with the baseline
model. We calculate the root mean square error $(\mathrm{RMSE})$ of the
transformed systems over time period $t$ as
$\mathrm{RMSE}:=\sqrt{\tfrac{\sum_{k=1}^{t}\bm{e}_{k}^{2}}{t}}$, where
$\bm{e}_{k}:=\left\lvert\bm{\hat{x}}_{k}-\bm{x}_{k}\right\rvert$ is the
difference between the states of the baseline DAE model (1) denoted as
$\bm{x}_{k}$ and the states of proposed systems as $\bm{\hat{x}}_{k}$.

Fig. 2 presents RMSE for the $\mathrm{Approx}$-$\mathrm{DAE}$ model when
varying the choice of $\mu$ between $\\{1\times 10^{-04},1\times 10^{-08}\\}$
under continuous and discrete time-domain simulations. It is illustrated that
for both test cases, the choice of $\mu$ affects the error on state
trajectories. For TI and BE discretization method, the choice of $\mu$ does
not alter the RMSE error. Such methods are single step discretization methods
which are inherently less accurate than the multi-step method, BDF, that is
typically used for DAE systems. The RMSE for the continuous and BDF case is
lower and is influenced by the choice of $\mu$ until it becomes asymptotic
after a certain value of $\mu$. For the remainder of this work $\mu$ is chosen
to be $1\times 10^{-06}$.

(a)  (b)

Figure 4: RMSE on transient dynamic and algebraic states under varying load
disturbances.

To assess the accuracy of the discrete time-domain simulations, we compute the
error norm $\zeta:=\left\lVert\hat{\bm{x}}_{k}-{\bm{x}}_{k}\right\rVert_{2}$
between the baseline states $\bm{x}_{k}$ at time and the states resulting from
the proposed models $\bm{\hat{x}}_{k}$. Fig. 3 illustrates the error resulting
from transient discrete time-domain simulation under the different
discretization methods. The proposed models accurately depict the transient
states; it can be noted that BDF method outperforms BE and TI for the DAE
model and the proposed systems. To further investigate the applicability of
the methods, we vary the load disturbance magnitude $\alpha_{L}$ under the
aforementioned ranges. The RMSE for $\mathrm{case}$-$\mathrm{9}$ and
$\mathrm{case}$-$\mathrm{200}$ under such load disturbances are depicted in
Fig. 4.

The performance of the transformed systems yields similar results to that of
the DAE system under the continuous and discrete-time models while simulating
the states under different transient conditions. We note that the performance
under BDF discretization of the $\mathrm{ODE}$-$\mathrm{DAE}$ model, suggests
that the BDF order $k_{g}$ might require changing to account for the altered
stiffness in the dynamics. Refer to Fig. 5 for the computational complexity
resulting from simulating the proposed transformations relative to the
nonlinear DAE system. The DAE and $\mathrm{Approx}$-$\mathrm{DAE}$ models
require the same computational effort, however that of the
$\mathrm{ODE}$-$\mathrm{DAE}$ model shows around a $100\%$ to $200\%$ increase
in the computational complexity under discrete-time techniques. This is
primarily due to the Hessian approximations (10), whereby such an increase is
not evident at the continuous time-domain level. An alternative Taylor series
approximation for (10) can be chosen that requires less computational effort.
Having provided the above results, the validity and accuracy of the proposed
models are demonstrated, therefore proving their applicability for transient
time-domain simulations of multi-machine power systems.

##  V Summary and Future Work

This paper presents transformations that result in solutions to nonlinear DAE
power systems modeled as nonlinear ODEs. The validity and accuracy of the two
methods have been investigated. The proposed transformations yield accurate
depictions of the transient state-space models and can therefore be used in
various feedback control or state estimation algorithms—our future work on
this topic.

(a)  (b)

Figure 5: Computational time for the different time-domain methods and power
system models.

## References

  * [1] L. Aolaritei, D. Lee, T. L. Vu, and K. Turitsyn, “A Robustness Measure of Transient Stability under Operational Constraints in Power Systems,” _IEEE Control Systems Letters_ , vol. 2, no. 4, pp. 803–808, 2018.
  * [2] W. Cui and B. Zhang, “Equilibrium-Independent Stability Analysis for Distribution Systems with Lossy Transmission Lines,” _IEEE Control Systems Letters_ , vol. 6, pp. 3349–3354, 2022.
  * [3] 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.
  * [4] ——, “On Current and Power Injection Models for Angle and Voltage Stability Analysis of Power Systems,” _IEEE Transactions on Power Systems_ , vol. 31, no. 3, pp. 2503–2504, 2016.
  * [5] 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.
  * [6] J. Y. Astic, M. Jerosolimski, and A. Bihain, “The mixed adams - BDF variable step size algorithm to simulate transient and long term phenomena in power systems,” _IEEE Transactions on Power Systems_ , vol. 9, no. 2, pp. 929–935, 1994.
  * [7] R. K. Brayton, F. G. Gustavson, and G. D. Hachtel, “Differential-Algebraic Systems Using,” vol. 60, no. 1, pp. 98–108, 1972.
  * [8] C. W. Gear, “Simultaneous Numerical Solution of Differential-Algebraic Equations,” _IEEE Transactions on Circuit Theory_ , vol. 18, no. 1, pp. 89–95, 1971.
  * [9] F. Milano, Ed., _Advances in Power System Modelling, Control and Stability Analysis_ , 2nd ed. London: IET, 2022.
  * [10] 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.
  * [11] J. D. Lara, R. Henriquez-Auba, D. Ramasubramanian, S. Dhople, D. S. Callaway, and S. Sanders, “Revisiting Power Systems Time-domain Simulation Methods and Models,” _arXiv_ , pp. 1–11, 2023.
  * [12] A. J. Conejo and L. Baringo, _Power System operation_. Springer, 2018, vol. 14, no. 54.
  * [13] P. W. Sauer, M. A. Pai, and J. H. Chow, _Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox_ , second edi ed. Wiley-IEEE Press, 2017.
  * [14] 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\.
  * [15] 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.
  * [16] J. Jaiswal and N. K. Tomar, “Existence Conditions for ODE Functional Observer Design of Descriptor Systems,” _IEEE Control Systems Letters_ , vol. 6, no. 2, pp. 355–360, 2022.
  * [17] P. Kunkel and V. Mehrmann, “Optimal control for unstructured nonlinear differential-algebraic equations of arbitrary index,” _Mathematics of Control, Signals, and Systems_ , vol. 20, no. 3, pp. 227–269, 2008.
  * [18] M. Althoff and B. H. Krogh, “Reachability analysis of nonlinear differential-algebraic systems,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 2, pp. 371–383, 2014.
  * [19] Y. Liu and K. Sun, “Solving Power System Differential Algebraic Equations Using Differential Transformation,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 2289–2299, may 2020.
  * [20] M. Nadeem and A. F. Taha, “Robust Dynamic State Estimation of Multi-Machine Power Networks with Solar Farms and Dynamics Loads,” _Proceedings of the IEEE Conference on Decision and Control_ , vol. 2022-Decem, no. Cdc, pp. 7174–7179, 2022.
  * [21] P. Kunkel and V. Mehrmann, _Differential-Algebraic Equations: analysis and Numerical Solution_. EMS Press, 2005\.
  * [22] S. L. Campbell and C. W. Gear, “The index of general nonlinear DAEs,” _Numerische Mathematik_ , vol. 72, no. 2, pp. 173–196, 1995.
  * [23] 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.
  * [24] M. L. Crow, _Computational methods for electric power systems, third edition_ , 2015.
  * [25] S. G. Krantz and H. R. Parks, _The Implicit Function Theorem: History, Theory, and Applications_. Springer New York, 2013.
  * [26] J. H. Chow, J. R. Winkelman, M. A. Pai, and P. W. Sauer, “Singular perturbation analysis of large-scale power systems,” _International Journal of Electrical Power and Energy Systems_ , vol. 12, no. 2, pp. 117–126, 1990.
  * [27] N. H. McClamroch, “On the Connection Between Nonlinear Differential-Algebraic Equations and Singularly Perturbed Control Systems in Nonstandard Form,” _IEEE Transactions on Automatic Control_ , vol. 39, no. 5, pp. 1079–1084, 1994.
  * [28] F. Shen, P. Ju, M. Shahidehpour, Z. Li, C. Wang, and X. Shi, “Singular Perturbation for the Dynamic Modeling of Integrated Energy Systems,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 1718–1728, 2020\.
  * [29] R. E. Mickens, _Nonstandard Finite Difference Models of Differential Equations_. World Scientific Publishing, 1994.
  * [30] S. P. Venkateshan and P. Swaminathan, “Chapter 8 - Numerical Differentiation,” in _Computational Methods in Engineering_ , S. P. Venkateshan and P. Swaminathan, Eds. Boston: Academic Press, 2014, pp. 295–315.