PINNSim: A Simulator for Power System Dynamics based on Physics-Informed
  Neural Networks
The dynamic behaviour of a power system can be described by a system of
differential-algebraic equations. Time-domain simulations are used to simulate
the evolution of these dynamics. They often require the use of small time step
sizes and therefore become computationally expensive. To accelerate these
simulations, we propose a simulator - PINNSim - that allows to take
significantly larger time steps. It is based on Physics-Informed Neural
Networks (PINNs) for the solution of the dynamics of single components in the
power system. To resolve their interaction we employ a scalable root-finding
algorithm. We demonstrate PINNSim on a 9-bus system and show the increased time
step size compared to a trapezoidal integration rule. We discuss key
characteristics of PINNSim and important steps for developing PINNSim into a
fully fledged simulator. As such, it could offer the opportunity for
significantly increasing time step sizes and thereby accelerating time-domain
simulations.
Jochen Stiasny, Baosen Zhang, Spyros Chatzivasileiadis
Solving Differential-Algebraic Equations in Power Systems Dynamics with Neural Networks and Spatial Decomposition

###### Abstract

The dynamics of the power system are described by a system of differential-algebraic equations. Time-domain simulations are used to understand the evolution of the system dynamics. These simulations can be computationally expensive due to the stiffness of the system which requires the use of finely discretized time-steps. By increasing the allowable time-step size, we aim to accelerate such simulations. In this paper, we use the observation that even though the individual components are described using both algebraic and differential equations, their coupling only involves algebraic equations. Following this observation, we use Neural Networks (NNs) to approximate the components' state evolution, leading to fast, accurate, and numerically stable approximators, which enable larger time-steps. To account for effects of the network on the components and vice-versa, the NNs take the temporal evolution of the coupling algebraic variables as an input for their prediction. We initially estimate this temporal evolution and then update it in an iterative fashion using the Newton-Raphson algorithm. The involved Jacobian matrix is calculated with Automatic Differentiation and its size depends only on the network size but not on the component dynamics. We demonstrate this NN-based simulator on the IEEE 9-bus test case with 3 generators.

## I Introduction

Solving Differential-Algebraic Equations (DAEs) constitutes one of the indispensable computational tasks for analyzing power system dynamics, e.g., for the analysis of transient stability phenomena. The size and complexity of the power system make solving these DAEs computationally demanding. In addition, the ongoing transition of the power system towards more distributed and heavily inverter-based generation increases this burden even further. Unsurprisingly, the acceleration of DAE simulations has therefore been a research topic throughout many years. The fundamental problem setup as described in [1] has undergone little change in the past decades. The main lines of work, as [2] reports, have been on 1) simplifying models, 2) decomposing the problem (spatially, temporally, and across the numerical schemes; see also [3]) to allow for parallelization, and 3) using pre-computed analytical approximations for faster online evaluation. Additionally, the use of parallel processing offers potential for acceleration [4].

In recent years, using Machine Learning (ML) to solve differential equations has emerged as an active area of research. ML-based methods can be evaluated in real-time and have the capacity to approximate complicated functions, including ones that usually cannot be expressed in closed form. Physics-Informed Neural Networks (PINNs), as reviewed in [5], are used for approximating the solution of partial differential equation. The ideas of domain decomposition are used in [6, 7] and can extend beyond PINNs.

In the area of power system dynamics, works in [8, 9, 10] have used ML-based methods to predict transient behaviors of these systems. They effectively provide approximate DAE-solvers. However, it is unclear whether these approaches are scalable with respect to the system size and different network configurations without becoming prohibitively expensive to train. A potential alternative to address this problem of scalability is the development of hybrid methods, that combine ML-based approaches with conventional numerical DAE-solution schemes to leverage their respective advantages. In [11], the authors train a NN to approximate the unknown dynamics of a system component which can then be integrated into a DAE-solver. But this approach does not fundamentally overcome the scaling challenges in solving DAEs.

In order to enable larger time-step sizes and hence to accelerate the simulation, we focus in this paper on the use of NNs as approximators for the state evolution of single components; we thereby essentially apply a spatial decomposition to the system. In contrast to other explicit approximators for the state evolution of components, NNs are fast, highly accurate, and numerically stable. However, it is nontrivial to model the interactions of these components in the network. The solution approach, we present, is structurally similar the Semi-Analytical Solution (SAS)-based method in [12].

The SAS-based methods make use of time-power series [12] or Adomian Decomposition [13] to approximate the state evolution. An essential part in these approaches is to model the interactions between various components (e.g., how generators interact with each other through power flow along the lines). Since the components are coupled with each other through algebraic equations, one way to "anticipate" their interaction is to predict or estimate the evolution of the algebraic variables as a function of time. For example, finding an approximate power series expansion in the time variable can lead to larger time-steps and therefore accelerate the simulations [2, 12]. Hence the quality of these approximations is of critical importance.

In this paper, we propose a NN-based simulator for power system dynamics. The simulator consists of a NN for each
dynamic system component and a network model that governs their interaction. To model these interactions, the NNs use the temporal evolution of the network algebraic variables as an input for their prediction. As the temporal evolution of the network algebraic variables is initially unknown - when it is known, we "solved" the problem - we first estimate it by means of a power series with respect to time. We then formulate an optimization problem to update this estimate by matching it with the prediction of the NNs.

Section II introduces the problem formulation and Section III describes the proposed solution approach. Section IV introduces the case study and Section V the results. Section VI discusses the results and Section VII concludes.

## II Problem formulation

We consider an autonomous, semi-explicit system of DAEs

\[\frac{d}{dt}\mathbf{x} =\mathbf{f}(\mathbf{x}(t),\mathbf{y}(t),\mathbf{u}) \tag{1a}\] \[\mathbf{0} =\mathbf{g}(\mathbf{x}(t),\mathbf{y}(t)) \tag{1b}\]

with differential variables \(\mathbf{x}(t)\in\mathbb{R}^{m}\), algebraic variables \(\mathbf{y}(t)\in\mathbb{R}^{n}\), control inputs \(\mathbf{u}\in\mathbb{R}^{l}\), and the functions \(\mathbf{f}:\mathbb{R}^{m+n+l+1}\mapsto\mathbb{R}^{m}\) and \(\mathbf{g}:\mathbb{R}^{m+n+1}\mapsto\mathbb{R}^{n}\). We furthermore assume that \(\mathbf{g}\) is continuously differentiable with respect to \(\mathbf{y}\), hence, this system of DAEs is of Hessenberg index-1 form [14]. Given the initial state \(\mathbf{x}_{0}=\mathbf{x}(t_{0})\) of the system, our goal is to evaluate the trajectory of the dynamic variables \(\mathbf{x}(t)\) and of the algebraic variables \(\mathbf{y}(t)\) for the time interval \(t\in[t_{0},t_{0}+\Delta t]\).

The DAE in (1) is nontrivial to solve and there has been a great deal of work dedicated to solving them in many applications [15, 16]. In this paper, we consider solving (1) in the context of power systems, where \(\mathbf{y}(t)\) describes bus terminal voltages and \(\mathbf{x}(t)\) includes variables related to dynamics within the power system components. For example, \(\mathbf{x}(t)\) includes the internal states of the generator (see, [17] and the references within). A particular challenge about the DAEs that describe power systems dynamics is that they can be stiff, and a number of techniques have been developed to address this challenge [1, 18, 19, 20]. These techniques, however, tend to be computationally intensive and are difficult to use in real-time or when a large number of simulations need to be performed in a short time (e.g., for systems with high level of renewable resources).

To overcome these computational challenges, we leverage the underlying sparsity in power systems. Namely, we use the fact that the differential variables associated with each bus, evolve based only on the variables at that bus, except for a coupling constraint on the algebraic variables.1 This coupling constraint is the current balance, derived from Kirchhoff's and Ohm's laws, and constitutes (1b):

Footnote 1: In this paper, we do not consider the very fast electromagnetic transients on the transmission lines.

\[\mathbf{0}=\bar{Y}_{N}\bar{\mathbf{V}}(t)-\bar{\mathbf{I}}(\mathbf{x}(t),\bar{\mathbf{V}}(t)). \tag{2}\]

\(\bar{\mathbf{V}}\in\mathbb{C}^{n}\) contains the complex voltages at \(n\) buses in the system and \(\bar{Y}\in\mathbb{C}^{n\times n}\) defines the complex admittance matrix of the network. The complex current injections from the \(n\) buses, i.e., from components such as synchronous generators, converters, or loads, are collected in the vector \(\bar{\mathbf{I}}\in\mathbb{C}^{n}\). We assume \(m\) components to be connected to the system with the differential states \(\mathbf{x}_{i}\in\mathbb{R}^{p}\) and control inputs \(\mathbf{u}_{i}\in\mathbb{R}^{q}\) for the \(i\)-th component.2 The \(i\)-th component is connected at the \(i\)-th bus and therefore dependent only on the complex voltage \(\bar{V}_{i}\in\mathbb{C}^{1}\). Each component is governed by a set of differential equations \(\mathbf{f}_{i}:\mathbb{R}^{p+2+1+q}\mapsto\mathbb{R}^{p}\) and a function \(h_{i}:\mathbb{R}^{p+2}\mapsto\mathbb{C}^{1}\) that determines the current injection \(\bar{I}_{i}\)

Footnote 2: All the results in the paper hold if \(p\) is different for each bus \(i\), but we avoid this notation to reduce clutter in the variables.

\[\frac{d}{dt}\mathbf{x}_{i} =\mathbf{f}_{i}(\mathbf{x}_{i}(t),\bar{V}_{i}(t),\mathbf{u}_{i}), i=1,\ldots,m \tag{3a}\] \[\bar{I}_{i} =h_{i}(\mathbf{x}_{i}(t),\bar{V}_{i}(t)), i=1,\ldots,m \tag{3b}\]

The concatenation of all differential variables is of dimension \(m\cdot p\), i.e., \(\mathbf{x}\in\mathbb{R}^{m\cdot p}\).

## III Methodology

We approach the solution of (2) and (3) by finding an explicit expression for \(\bar{\mathbf{V}}(t)\) for \(t\in[t_{0},t_{0}+\Delta t]\). Given \(\bar{\mathbf{V}}(t)\), the solution of (3a) can be computed quickly, as we now only need to solve \(m\) much smaller and simpler problems instead of a single large and difficult one. We start by initially estimating the temporal evolution of the algebraic variables \(\hat{\hat{\mathbf{V}}}(t)\) (Section III-A). Using this estimate, we employ a NN for each component \(i\) to predict the trajectory \(\hat{\mathbf{x}}_{i}(t)(t,\mathbf{x}_{0},\hat{\bar{\mathbf{V}}}(t))\) (Section III-B). We then can evaluate (2) to assess the quality of the guess \(\hat{\bar{\mathbf{V}}}(t)\) and update it iteratively based on the mismatch of the currents in (2) (Section III-C). As soon as this update scheme converged, we "solved" the problem for the present time-step and can repeat the procedure for the next time-step. We are thereby able to simulate trajectories of arbitrary length (Section III-D).

This approach allows us, firstly, to incorporate NNs, which are a powerful explicit function approximator, into solution algorithms for DAEs. NNs can approximate \(\hat{\mathbf{x}}(t)\) well for large time-steps of complicated functions, they do not experience numerical instability, and they are fast to evaluate. The second benefit relates to the fact, that the proposed solution algorithm is indifferent to the dimensionality of the differential variable vector, i.e., \(\dim(\mathbf{x})\). This addresses a major challenge for solving DAEs in very large power systems as \(\dim(\mathbf{x})\) rapidly grows when many complicated components are present.

### _Parametrization of the voltage estimate \(\hat{\hat{\mathbf{V}}}(t)\)_

Based on the initial state \(\mathbf{x}_{0}\), we can compute the initial voltages \(\bar{\mathbf{V}}\) that satisfy the current balance. We now take an initial estimate for the evolution of the complex voltage in its polar form at each bus for the time interval \([t_{0},t_{0}+\Delta t]\). We express this estimate as a power-series with respect to
time for the voltage magnitude \(V_{i}\) and the voltage angle \(\theta_{i}\)

\[\bar{V}_{i}(t) =V_{i}(t)e^{j\theta_{i}(t)} \tag{4}\] \[\approx\left(\sum_{k=0}^{r}V_{k,i}(t-t_{0})^{k}\right)e^{j\left( \sum_{k=0}^{r}\theta_{k,i}(t-t_{0})^{k}\right)} \tag{5}\]

up to power \(r\). The coefficients \(V_{0,i},V_{1,i},\ldots,V_{r,i}\) and \(\theta_{0,i},\theta_{1,i},\ldots,\theta_{r,i}\) form the parameters which we will later on update to improve the initial estimate. We subsequently refer to this estimate as \(\hat{V}_{i}(t,\mathbf{\Xi}_{i})\). The vector \(\mathbf{\Xi}_{i}\) collects all parameters at bus \(i\)

\[\mathbf{\Xi}_{i}=\begin{bmatrix}V_{0,i}&\theta_{0,i}&\ldots&V_{r,i}, \theta_{r,i}\end{bmatrix}\in\mathbb{R}^{1\times 2r}. \tag{6}\]

This parametization is repeated for all \(n\) buses in the system. We collect all \(\mathbf{\Xi}_{i}\) in the vector \(\mathbf{\Xi}\in\mathbb{R}^{2rn}\).

### _Explicit approximation of \(\hat{\mathbf{x}}_{i}(t)\) by a Neural Network (NN)_

The exact solution for the evolution of the differential variables \(\mathbf{x}_{i}(t)\) of machine \(i\) can be obtained by integration of (3a)

\[\mathbf{x}_{i}(t)=\mathbf{x}_{0,i}+\int_{t_{0}}^{t}\mathbf{f}_{i}\left(\mathbf{x}_{i}(\tau), \bar{V}_{i}(\tau),\mathbf{u}_{i}\right)d\tau. \tag{7}\]

As there usually exists no explicit analytical solution to (7), we seek to approximate the solution. In a first step, we replace \(\bar{V}_{i}(\tau)\) by the estimate for the bus voltage \(\hat{\bar{V}}_{i}(\tau)\)

\[\mathbf{x}_{i}(t)\approx\mathbf{x}_{0,i}+\int_{t_{0}}^{t}\mathbf{f}_{i}\left( \mathbf{x}_{i}(\tau),\hat{\bar{V}}_{i}(\tau),\mathbf{u}_{i}\right)d\tau. \tag{8}\]

Because an explicit analytical solution is usually still not attainable, we now approximate (8) by a NN to obtain

\[\hat{\mathbf{x}}_{i}^{NN}(t)=NN(t,\mathbf{x}_{0,i},\hat{\bar{V}}_{i}(t), \mathbf{u}_{i}). \tag{9}\]

As we have parameterized the estimate \(\hat{\bar{V}}_{i}(t)\) in terms of the parameters \(\mathbf{\Xi}_{i}\), the input \(\mathbf{z}_{0}\) to the NN becomes

\[\mathbf{z}_{0}=[t,\mathbf{x}_{0,i},\mathbf{\Xi}_{i},\mathbf{u}_{i}],\quad\mathbf{z}_{0}\in\mathbb{ R}^{1+p+2r+q} \tag{10}\]

We use a multi-layer perceptron NN which is a sequence of linear transformations parameterized by weight matrices \(\mathbf{W}_{k}\) and bias vectors \(\mathbf{b}_{k}\) and the application of a non-linear function \(\sigma\) in each of the \(K\) hidden layers

\[\mathbf{z}_{k+1} =\sigma\left(\mathbf{W}_{k+1}\mathbf{z}_{k}+\mathbf{b}_{k+1}\right),\,\forall k =0,1,\ldots,K-1 \tag{11a}\] \[\hat{\mathbf{x}}_{i}^{NN} =\mathbf{x}_{0,i}+(\mathbf{W}_{K}\mathbf{z}_{K}+\mathbf{b}_{K}). \tag{11b}\]

In our simulation, we use \(\sigma(\cdot)=\tanh(\cdot)\) and three layers, i.e., \(K=3\). As a baseline, we consider an explicit Runge-Kutta (RK)-based approximation scheme for (8)

\[\hat{\mathbf{x}}_{i}^{RK}(\Delta t,\mathbf{x}_{0,i},\mathbf{\Xi}_{i})=\mathbf{x}_ {0,i}+\Delta t\sum_{j=1}^{\nu}b_{j}\mathbf{k}^{(j)} \tag{12}\] \[\mathbf{k}^{(j)}=\mathbf{f}\left(t_{0}+c_{j}\Delta t,\mathbf{x}_{0,i}+\Delta t \sum_{l=1}^{\nu}a_{jl}\mathbf{k}^{(l)},\mathbf{\Xi}_{i},\mathbf{u}_{i}\right), \tag{13}\]

where \(a_{jl}\) is strictly lower triangular. Other choices such as time-power series [12] or Adomian decomposition methods [13] are equally possible. In either case the state evaluation is explicitly dependent on \(\mathbf{\Xi}_{i}\) which together with \(t\) accounts for the evolution of the network algebraic variables.

### _Update scheme for voltage parametrization_

Independent of the approximation method used for \(\hat{\mathbf{x}}_{i}\), we can now compute each component's current injection \(\hat{\bar{I}}_{i}(\hat{\mathbf{x}}_{i},\hat{\bar{V}}_{i})\). We then define the current mismatch \(\mathbf{\varepsilon}_{I}^{(k)}\in\mathbb{R}^{2n}\) in iteration \(k\)

\[\mathbf{\varepsilon}_{I}^{(k)}:=\begin{bmatrix}\Re\left(\bar{Y}_{N}\hat{\mathbf{V}}(t, \mathbf{\Xi}^{(k)})-\hat{\bar{\mathbf{I}}}(t,\mathbf{x}_{0},\mathbf{\Xi}^{(k)})\right)\\ \Im\left(\bar{Y}_{N}\hat{\mathbf{V}}(t,\mathbf{\Xi}^{(k)})-\hat{\bar{\mathbf{I}}}(t,\mathbf{x} _{0},\mathbf{\Xi}^{(k)})\right)\end{bmatrix}. \tag{14}\]

We deliberately used the notation \(\hat{\bar{\mathbf{V}}}(t,\mathbf{\Xi}^{(k)})\) and \(\hat{\bar{\mathbf{I}}}(t,\mathbf{x}_{0},\mathbf{\Xi}^{(k)})\) to make clear that \(\mathbf{\varepsilon}_{I}^{(k)}\) is a function dependent only on \(t\), \(\mathbf{x}_{0}\), and \(\mathbf{\Xi}_{i}^{(k)}\).

Figure 1 illustrates the real part of the current injection of the component \(\bar{I}_{1}\) and the respective network current \((\bar{Y}_{N}\hat{\bar{\mathbf{V}}})_{1}\) at bus 1. The current mismatch \(\Re(\mathbf{\varepsilon}_{I,1}^{(k)})\) corresponds to area between the two curves. By updating the parameters \(\mathbf{\Xi}\) we reduce the mismatch as seen in each of the columns in Fig. 1. If \(\mathbf{\varepsilon}_{I}=\mathbf{0}\) across the entire time-step \([t_{0},t_{0}+\Delta t]\) and \(\hat{\mathbf{x}}\) accurately solves (3a), then we found a solution to the system of DAEs. In practice, we evaluate (14) on \(s\) collocation points within the time interval, \(\mathbf{T}=[t_{1},\ldots,t_{s}]\), which we can freely choose. The two columns in Fig. 1 show two possible choices for \(s=2\) and \(s=5\). Formally, the objective in order to solve the system of DAEs becomes

\[\min_{\mathbf{\Xi}}\left\|\begin{bmatrix}\mathbf{\varepsilon}_{I}(t_{1}, \mathbf{x}_{0},\mathbf{\Xi})\\ \vdots\\ \mathbf{\varepsilon}_{I}(t_{s},\mathbf{x}_{0},\mathbf{\Xi})\end{bmatrix}\right\|_{2} \tag{15}\]

and we solve this optimization problem by applying an iterative procedure, summarized by Algorithm 1, to determine \(\mathbf{\Xi}\) based on the Newton-Raphson algorithm. The update value \(\Delta\mathbf{\Xi}^{(k)}\) in Step 8 of Algorithm 1 is determined by the least

Fig. 1: Current with \(s\) collocation points at iteration 0 and 5.
square problem

\[\min_{\Delta\mathbf{\Xi}^{(k)}}\ \left\|A\Delta\mathbf{\Xi}^{(k)}-B\right\|_{2} \tag{16a}\] \[A=\begin{bmatrix}\nicefrac{{\partial\mathbf{\varepsilon}_{I}(t_{1},\mathbf{x }_{0},\mathbf{\Xi}^{(k)})}}{{\partial\mathbf{\Xi}^{(k)}}}\end{bmatrix}B=\begin{bmatrix} \mathbf{\varepsilon}_{I}(t_{1},\mathbf{x}_{0},\mathbf{\Xi}^{(k)})\\ \vdots\\ \mathbf{\varepsilon}_{I}(t_{s},\mathbf{x}_{0},\mathbf{\Xi}^{(k)})\end{bmatrix} \tag{16b}\]

which we solve for repeatedly until we reach the defined tolerance \(\Delta\mathbf{\Xi}^{\max}\) or the maximum number of iterations \(k^{\max}\). The problem in (16) can be under-determined, over-determined, or unique, primarily depending on \(s\) and the number of voltage parameters \(\mathbf{\Xi}\). As we are free to choose the number of collocation points, we impose \(2\cdot n\cdot s\geq|\mathbf{\Xi}|\) which rules out under-determined settings unless the Jacobian matrices are not full rank3.

Footnote 3: The Jacobian of the network current \(\hat{\mathbf{V}}\hat{\mathbf{V}}\) is actually deficient by \(1\) rank due to the symmetry around the voltage angles, however, as long as the Jacobian of the component’s current injection matrix is non-zero, i.e., there is some dependency of the voltage parameters, the problem is well posed.

```
0:\(\mathbf{x}_{0}\), \(T\)
0:\(\mathbf{\Xi}^{(0)}\), \(\Delta\mathbf{\Xi}^{\max}\), \(k^{\max}\)
1:while\(\Delta\mathbf{\Xi}^{(k)}>\Delta\mathbf{\Xi}^{\max}\) and \(k<k^{\max}\)do
2:for component \(i=1,\ldots,m\)do
3: Calculate current injections \(\bar{I}_{i}^{(k)}=h_{i}(\mathbf{T},\mathbf{x}_{0,i},\mathbf{\Xi}_{i}^{(k)})\)
4: Calculate Jacobian \(\nicefrac{{\partial\bar{I}_{i}^{(k)}}}{{\partial\mathbf{\Xi}_{i}^{(k)}}}\)
5:end
6: Calculate network injections \(\bar{Y}_{N}\hat{\mathbf{V}}(\mathbf{T},\mathbf{\Xi}^{(k)})\)
7: Calculate network Jacobian \(\nicefrac{{\partial Y_{N}\hat{\mathbf{V}}(\mathbf{T},\mathbf{\Xi}^{(k)})}}{{\partial\mathbf{ \Xi}^{(k)}}}\)
8: Calculate current balance error \(\varepsilon_{I}^{(k)}\)
9: Calculate voltage parameter update \(\Delta\mathbf{\Xi}^{(k)}\)
10: Update iteration \(\mathbf{\Xi}^{(k+1)}=\mathbf{\Xi}^{(k)}+\Delta\mathbf{\Xi}^{(k)}\), \(k=k+1\)
11:end
12:return\(\mathbf{\Xi}\)
```

**Algorithm 1** Voltage parametrization update

Usually, \(\nicefrac{{\partial\hat{\mathbf{I}}}}{{\partial\mathbf{\Xi}^{(k)}}}\) is not straightforward to derive, however, the use of Automatic Differentiation (AD) allows an efficient computation of the derivatives [21]. AD constructs a computational graph for the calculation of \(\hat{\mathbf{x}}\) and \(\hat{\bar{\mathbf{I}}}\) and then applies the chain rule to obtain the the derivatives \(\nicefrac{{\partial\hat{\mathbf{x}}}}{{\partial\mathbf{\Xi}^{(k)}}}\) and \(\nicefrac{{\partial\hat{\mathbf{I}}}}{{\partial\mathbf{\Xi}^{(k)}}}\). As we restricted \(\hat{\mathbf{x}}(t,\mathbf{x}_{0},\mathbf{\Xi})\) to be an explicit integration scheme for (8), the associated computational graph is automatically constructed and of relatively small size which results in fast calculations. The calculation furthermore does not directly depend on the dimension of \(\mathbf{x}\).4

Footnote 4: \(\dim(\mathbf{x})\) might affect the necessary size of the computational graph.

### _Multi-step simulator for power system dynamics_

The sections above yield a state trajectory for a single time-step and its accuracy is dependent on the approximation quality of \(\hat{\mathbf{x}}(t)\) and \(\hat{\mathbf{V}}(t)\). Requirements on the resulting tolerance therefore limit the suitable time-step size \(\Delta t\). By repeatedly applying the above algorithm, we obtain a multi-step scheme that then allows the simulation of dynamics beyond \(\Delta t\). Algorithm 2 summarizes the simplest form of such a multi-step simulator in which we evaluate the trajectory on the points \(T_{\text{eval}}\). We show in Section V that Algorithm 2 yields accurate trajectories and that the use of NNs allows larger time-steps, hence faster simulations, while being more accurate than the RK-based approximators.

```
0:\(\mathbf{x}_{0}\), \(t_{0}\), \(t_{\max}\), \(\Delta t\)
0:\(\mathbf{x}^{\prime}_{0}=\mathbf{x}_{0}\), \(t^{\prime}=t_{0}\), \(T_{\text{eval}}\)
1:while\(t^{\prime}\leq t_{\max}\)do
2: Define collocation points \(\mathbf{T}(t^{\prime},\Delta t)\)
3: Calculate voltage parametrization \(\mathbf{\Xi}(\mathbf{T},\mathbf{x}^{\prime}_{0})\)
4: Store trajectory \([T_{\text{eval}},\hat{\mathbf{x}}(\mathbf{x}^{\prime}_{0},T_{\text{eval}},\mathbf{\Xi})]\)
5: Update \(\mathbf{x}^{\prime}_{0}\leftarrow\hat{\mathbf{x}}(\mathbf{x}^{\prime}_{0},\Delta t,\mathbf{ \Xi})\) and \(t^{\prime}\gets t^{\prime}+\Delta t\)
6:end
7:return Trajectory
```

**Algorithm 2** Multi-step simulator

## IV Case Study

This section describes the modeling of the power system dynamics and the implementation of the different elements of the simulator.

### _Power system modeling_

As an example for a dynamic component, we consider a two-axis generator model as modeled in [19]

\[\begin{bmatrix}T^{\prime}_{d\sigma}\\ T^{\prime}_{q\sigma}\\ 1\\ 2H\end{bmatrix}\ \frac{d}{dt}\begin{bmatrix}E^{\prime}_{q}\\ E^{\prime}_{d}\\ \delta\\ \Delta\omega\end{bmatrix}=\begin{bmatrix}-E^{\prime}_{q}-(X_{d}-X^{\prime}_{d})I_{d }+E_{fd}\\ -E^{\prime}_{d}+(X_{q}-X^{\prime}_{q})I_{q}\\ \omega_{\sigma}\Delta\omega\\ P_{m}-E^{\prime}_{d}I_{d}-E^{\prime}_{q}I_{q}-(X^{\prime}_{q}-X^{\prime}_{d})I_{ d}I_{q}-D\Delta\omega\end{bmatrix} \tag{17a}\] \[\begin{bmatrix}I_{d}\\ I_{q}\end{bmatrix}=\begin{bmatrix}R_{s}&-X^{\prime}_{q}\\ X^{\prime}_{d}&R_{s}\end{bmatrix}^{-1}\begin{bmatrix}E^{\prime}_{q}-V\sin{(\delta -\theta)}\\ E^{\prime}_{q}-V\cos{(\delta-\theta)}\end{bmatrix}\] (17b) \[\bar{I}=(I_{D}+jI_{Q})=(I_{d}+jI_{q})e^{j(\delta-\pi/2)}. \tag{17c}\]

(17a) corresponds to (3a) and (17b) and (17c) to (3b) For this study, we simplify above model to a classical model by setting the reactances to \(X^{\prime}_{q}=X^{\prime}_{d}\) and \(X_{q}=X^{\prime}_{d}\) and then finding the integral manifold such that the internal voltages \(E^{\prime}_{q}\) and \(E^{\prime}_{d}\) remain constant at \(E^{\prime}_{q0}\) and \(E^{\prime}_{d0}=0\). For more details we refer to [19]. Now, the rotor angle \(\delta\) and the frequency deviation \(\Delta\omega\) form the state \(\mathbf{x}_{i}\), the magnitude \(V\) and angle \(\theta\) of the terminal voltage form \(\mathbf{y}_{i}\), and the mechanical power \(P_{m}\) and the excitation voltage \(E_{fd}\) form the control input \(\mathbf{u}_{i}\). More detailed models could include higher order electro-mechanical modes, governor dynamics for \(P_{m}\) and exciter dynamics \(E_{fd}\). Similarly, inverter-based resources or voltage dependent loads could be included. The results presented in Section V-A are based on the parameters of generator 1 in Table I.

To illustrate the full multi-step simulator, we consider the IEEE 9-bus system described in [19, pp. 164-167] with three generators (all modeled as classical model as above) with parameters from Table I. The initial conditions \(\mathbf{x}_{0}\) and control inputs \(\mathbf{u}\) are determined from assuming an equilibrium state for the used load flow case. To perturb the system, we reduce the mechanical power \(P_{m}\) of generator 2 to 80% of its initial
value and then observe the trajectory for \(2.5\,\mathrm{s}\) as shown in Fig. 3.

### _State predictors_

We test three state predictors, namely a NN with 3 layers and 32 neurons per layer, and two RK-based predictors. They are the Forward-Euler (FE), and 4th-order Runge-Kutta (RK4) schemes. In all cases the complex voltage is approximated by linear functions, i.e., \(r=1\) in (5) resulting in \(|\mathbf{\Xi}|=4\).

### _Implementation and NN training_

A central functionality we use is AD, therefore, the entire simulator is implemented in PyTorch [22]. It furthermore provides the ML framework to train the NN models. For each generator we used 2500 simulated points, that we split into training (64%), validation (16%) and testing (20%) datasets. Each generator surrogate model is trained for a range of possible linear voltage parametrizations, initial conditions, and prediction time-steps, which results in a seven dimensional input space. The training function includes a physics regularization term with evaluates how well the prediction aligns with the governing equations (3a), known as PINNs. We train, using a L-BFGS optimizer, for 2000 epochs in about 10 minutes. The published code base [23] provides a detailed overview of the dataset creation setup and the training procedure. The timing of the simulations is conducted on a AMD Ryzen 7 PRO (1.9 GHz, 8 cores) and 16GB RAM. The numerical simulations of the system of DAEs are implemented using Assimulo [24].

## V Results

We begin by briefly comparing the performance of the three predictors for \(\hat{\mathbf{x}}_{i}\) (FE, RK4, NN) before testing Algorithm 2.

### _Predictor comparison_

The approximator \(\hat{\mathbf{x}}_{i}\) should have two properties: being fast and being accurate for large time-steps \(\Delta t\). Figure 2 shows the two properties. The left panel displays the error in the differential variable \(\Delta\omega\) across 200 points of random initial conditions \(\mathbf{x}_{0}\) and voltage parameterizations \(\mathbf{\Xi}_{i}\). In terms of accuracy, the NN performs well and only for smaller time-steps (\(\Delta t<0.05\mathrm{s}\)), the RK4 approximation becomes more accurate. The RK-schemes exhibit the expected dependency of the time-step - the local truncation error should follow \(\mathcal{O}(\Delta t)\) and \(\mathcal{O}(\Delta t^{4})\) for FE and RK4. The error of the NN in contrast is near independent of \(\Delta t\). At the same time, the NN is the fastest approximator. While the FE approximator is similarly fast, its poor accuracy makes it undesirable and while more accurate, the RK4 scheme has the drawback of comparably long run-time.

### _Simulator results_

We now focus on the performance of the simulator in Algorithm 2. Figure 3 displays the prediction for \(\Delta t=0.05\,\mathrm{s}\) using NN-based approximators \(\hat{\mathbf{x}}_{i}^{NN}\).

The results correspond to the first row in Table II where we report the maximum absolute errors of \(V\) and \(\Delta\omega\) along the trajectory and the run-time. Table II shows further results for the RK4-based simulator, for time-steps of \(\Delta t=0.1\,\mathrm{s}\) and \(\Delta t=0.15\,\mathrm{s}\) and for collocation points at \(\mathbf{T}=[0.3,0.7]\Delta t\) and \(\mathbf{T}=[0.1,0.3,0.5,0.7,0.9]\Delta t\), i.e., \(s=2\) and \(s=5\). The NN-based simulator is consistently faster and more accurate than the RK-based simulator, except for the case of \(\Delta t=0.05\,\mathrm{s}\). These results confirm the predictor characteristics observed in Section V-A. The simulation run-time scales approximately inversely proportional with \(\Delta t\). Deviations can arise due to varying numbers of iteration in Algorithm 1, however, in the reported cases, we observe usually 5-7 iterations. Increasing

Fig. 3: Prediction of a state trajectory (\(\Delta\omega_{i}\)) and an algebraic variable \(V_{i}\) for time-step size \(\Delta t=0.05\,\mathrm{s}\) with a NN-based simulator. The predictions (dashed lines) coincide with the ground truth (gray lines).

Fig. 2: Predictor characteristics: (left) prediction error of \(\Delta\omega\) for 200 predictions with random \(\mathbf{x}_{0}\) and \(\mathbf{\Xi}\), (right) run-time per point. These results show that NN constitute an accurate and fast predictor.
the number of collocation points \(s\) results in a small increase in run-time in all cases. In terms of accuracy, we observe that more collocation points can lead to better accuracy, when the overall solution quality is good. However, for too large time-steps, here \(\Delta t=0.15\,\mathrm{s}\), the effect might reverse. The choice of the location of the collocation points, i.e., \(\mathbf{T}\), also matters.

## VI Discussion

The above results illustrate the conceptual benefits of a simulator leveraging NN-based predictors for the dynamic components of a power system. The upside of the approach should become more clear when integrating higher order models that would usually require very small time-steps for accurate resolution. A second direction of improvements lies in the implementation. All dynamic elements should be pre-compiled and distributed on the available cores to speed up the calculation of the residual and the Jacobian of the residual as this is the area where most computational time is spend. Furthermore, borrowing from classical numerical methods, dishonest Newton-Raphson-schemes and the exploitation of the structure in the Jacobians should be improved. The latter is due to the sparsity in the network equations and the dependence of the components current injection only on its terminal voltage. Another important aspect is to investigate the relation between the number of collocation points, the resulting current balance errors and overall errors. Analytical analysis of the scheme will be necessary to build confidence in the accuracy but also to identify limitations of the approach. Ultimately, larger test cases need to be studied.

## VII Conclusion

This work presented a simulation approach that centers around finding iteratively an approximation of the evolution of the algebraic variables in the power system DAEs. The approximation of the dynamic state evolutions by NNs, instead of classical explicit numerical integration schemes, allows larger time-steps to be realized while being fast to execute. This work aimed at providing a proof of concept, it is foreseeable that future work on this method shares many typical questions with established DAE solvers, hence, by applying various existing techniques the computational performance and scalability of the approach should improve significantly.
