# PINNSim: A Simulator for Power System Dynamics
based on Physics-Informed Neural Networks

Jochen Stiasny, Spyros Chatzivasileiadis  Division for Power and Energy
Systems
Technical University of Denmark
Kgs. Lyngby, Denmark
{jbest<EMAIL_ADDRESS>Baosen Zhang  Electrical and Computer Engineering
University of Washington
Seattle, WA 98195
<EMAIL_ADDRESS>
###### Abstract

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.

###### Index Terms:

Dynamical systems, Differential-algebraic equations, Physics-Informed Neural
Networks, Time-domain simulation

Submitted to the 23rd Power Systems Computation Conference (PSCC 2024). This
work was supported by the ERC Project VeriPhIED, funded by the European
Research Council, Grant Agreement No: 949899.

##  I Introduction

Differential equations form a centre piece in the modelling of the dynamic
behaviour of power systems. They provide very widely applicable component and
system models, however, their solution requires numerical integration methods.
These tools for Time-Domain Simulations (TDSs) constitute a workhorse of power
system analysis and can reliably deliver accurate solutions. The used
algorithms are scalable, versatile, and applicable without case-specific
modifications, but they incur a high computational burden [1]. Efforts around
improved solvers and parallelisation can reduce this burden, see [2, 3, 4] for
an overview. However, a fundamental problem is the small time step size that
is required by the schemes to ensure accurate results and numerical stability.
Hence many time steps and function evaluations are necessary.

The proposition of Neural Network (NN)-based approaches [5], and in particular
Physics-Informed Neural Networks (PINNs) [6], has very different
characteristics. In the training process, PINNs learn the solution to the
provided differential equations and subsequently provide fast and sufficiently
accurate solutions over long time intervals – addressing two main issues of
conventional numerical integration methods. PINNs have been adopted to power
system dynamics in [7] to predict single machine dynamics. In [8], the related
methodology of operator learning is also applied to single machine dynamics.
The applicability of learning-based methods to multi-machine systems is
explored in [9] and with alternative architectures in [10, 11]. The speed
advantage persists and the accuracy remains sufficiently high for the trained
setup. However, the common challenge lies in achieving generalisation, i.e.,
being able to apply the same NN for an unseen task without retraining. Such
flexibility would make NN-based TDS approaches comparable with conventional
numerical integration methods. However, in the case of multi-machine systems,
this requires the generalisation to a very large variety of setups. As a
consequence, the dimensionality of the associated learning tasks increases
dramatically with the system size and thereby becomes very or even too costly
to learn [12].

To nonetheless harness the advantages of PINNs, i.e., fast solutions for large
time steps, we propose PINNSim. This algorithm allows a modular integration of
PINNs into TDSs for multi-machine systems. We learn the solution of all
dynamic components in the system separately, i.e., with entirely independent
PINNs, and before executing any simulation. Thereby, the learning tasks are of
relatively low dimensionality and hence remain tractable. To perform a TDS, we
solely need to determine the interaction between the components. To this end,
we utilise a root-finding algorithm similar to conventional TDS algorithms; at
this stage, no learning is required. As a result, the scalability and
flexibility of PINNSim become decoupled from the difficulties of learning
high-dimensional problems. At the same time, we still can benefit from fast
and accurate solutions for larger time step sizes due to the use of PINNs.

In Section II, we present the conceptual motivation behind PINNSim. In Section
III, we introduce the power system formulation and the methodology of PINNSim.
Section IV describes the case study and in Section V, we discuss the results
by highlighting key characteristics of PINNSim. Section VI discusses the steps
for developing PINNSim from a proof of concept to a fully fledged simulator.
Section VII concludes.

##  II Concept

To facilitate the presentation of the methodology in Section III, we want to
first describe the conceptual idea behind PINNSim. It originates from the
problem that we face when solving Differential-Algebraic Equations (DAEs), a
form of differential equations where a set of algebraic equations constraints
the differential equations. It formulates as

$\displaystyle\frac{d}{dt}\bm{x}$
$\displaystyle=\bm{f}\left(\bm{x},\bm{y}\right)$ (1a) $\displaystyle\bm{0}$
$\displaystyle=\bm{g}\left(\bm{x},\bm{y}\right)$ (1b)

and we refer to $\bm{x}(t)$ as the differential variables and to $\bm{y}(t)$
as the algebraic variables. The update function
$\bm{f}\left(\bm{x},\bm{y}\right)$ and the algebraic relationship
$\bm{g}\left(\bm{x},\bm{y}\right)$ govern the dynamics of the system111For
notational clarity, we formulate an autonomous, unforced system. The
conceptual idea can also accommodate non-autonomous and forced systems.. Our
interest focuses on the particular form of index-1 DAEs or semi-analytical
DAEs [13]. This form implies that $\bm{g}$ can be differentiated once with
respect to time $t$ which is possible when
$\frac{\partial\bm{g}{}}{\partial\bm{y}}$ is non-singular. Based on 1, we
could describe the temporal evolution of $\bm{x}$ and $\bm{y}$ as

$\displaystyle\bm{x}(t)$
$\displaystyle=\bm{x}_{0}+\int_{t_{0}}^{t}\bm{f}\left(\bm{x},\bm{y}\right)d\tau$
(2a) $\displaystyle\bm{y}(t)$ $\displaystyle=\bm{g}^{\prime}(\bm{x})$ (2b)

where $\bm{x}_{0}=\bm{x}(t_{0})$ represents the initial condition and
$\bm{g}^{\prime}(\bm{x})$ describes the solution to the algebraic equations
given $\bm{x}$. However, usually no analytical expression describes the
evolution $\bm{x}(t)$ and $\bm{y}(t)$ for a given $\bm{x}_{0}$. Hence, we
revert to numerical integration methods to obtain an approximate solution.

To resolve the non-trivial integration operation and the implicit relationship
$\bm{g}$, numerical schemes often restrict the functional form of $\bm{x}$ to
a certain approximation $\hat{\bm{x}}$. For instance, Runge-Kutta (RK) methods
assume a polynomial form of
$\hat{\bm{x}}{}(t)=\bm{x}_{0}+a_{1}(t-t_{0})+a_{2}(t-t_{0})^{2}+\dots$ as they
match the Taylor expansion up to a certain degree by construction. The
different RK schemes prescribe the order of the scheme and the computation of
the coefficients. When algebraic variables are present, they have to be
interfaced with the approximation of the differential variable to incorporate
their interaction with each other. Simultaneous and partitioned integration
methods describe such routines [1] and are used in many variations. The
accuracy of these constructions is dependent on the order of the used
integration scheme and potential interface errors. These considerations and
aspects of numerical stability limit the usable time step size.

By choosing a different functional form for $\hat{\bm{x}}$, like Fourier-
series based or around the Adomian decomposition [14], the practical time step
size might be increased. Due to their high flexibility of their functional
form, PINNs, and NNs in general as suggested in [5], can allow significantly
larger time steps. In fact, we can even choose a functional form of
$\hat{\bm{y}}{}(t)$ and approximate
$\hat{\bm{x}}{}=\operatorname{PINN}(t,\hat{\bm{y}}{})$ in dependence of it.
However, this great approximation flexibility of PINNs comes with the
challenge of generalising well across the entire domain of interest, i.e.,
being accurate over the entire domain and not only on the training dataset.
When the dimensionality of $\bm{x}$ and $\bm{y}$ in 1 increases, the training
of a PINN to approximate a wide range of solutions becomes increasingly
difficult and eventually intractable. With PINNSim, we avoid this problem by
exploiting the structure of the power system specific DAEs. The structure
allows a decomposition of 1 into multiple smaller sub-problems which remain
tractable from a learning perspective. At the same time, the use of PINNs
allows for large and accurate time steps of these sub-problems. To simulate
the entire system, we need to align the sub-problems’ solutions by enforcing
the algebraic relationship $\bm{g}\left(\bm{x}{},\bm{y}{}\right)$.

##  III Methodology

This section describes the problem formulation and its decomposition for
PINNSim in Section III-A. We then describe the main elements of the algorithm
in Sections III-B, III-C and III-D and lastly the entire algorithm in Section
III-E.

###  III-A Problem formulation and solution approach

The key relation for power system dynamics is described by the current
balance, i.e., Kirchhoff’s current law, and it needs to hold at all times $t$
and for each of the $n$ buses in the network. We distinguish between the
complex current injections stemming from the network
$\bar{\bm{i}}^{N}{}\in\mathbb{C}^{n}$ and from connected components
$\bar{\bm{i}}^{C}{}\in\mathbb{C}^{n}$. The former can be described by the
algebraic relationship (neglecting electro-magnetic transients)

$\displaystyle\bar{\bm{i}}^{N}{}=\bar{\bm{Y}}\bar{\bm{v}}.$ (3)

where $\bar{\bm{v}}\in\mathbb{C}^{n}$ represents the complex voltages at the
buses and $\bar{\bm{Y}}\in\mathbb{C}^{n\times n}$ the complex admittance
matrix. The component current injections $\bar{i}_{i}^{C}$ at bus $i$ can stem
from a dynamic component, i.e., their behaviour is governed by differential
equations,

$\displaystyle\bar{i}_{i}^{C}$ $\displaystyle=h_{i}(\bm{x}_{i},\bar{v}_{i})$
(4a) $\displaystyle\frac{d}{dt}\bm{x}_{i}$
$\displaystyle=\bm{f}_{i}\left(\bm{x}_{i},\bar{v}_{i},\bm{u}_{i}\right)$ (4b)

where the current injection depends on the state vector
$\bm{x}_{i}\in\mathbb{R}^{p}$ and the voltage $\bar{v}_{i}$ of the bus $i$ to
which the component is connected222For notational ease, we assume that
component $i$ is connected to bus $i$. If component $i$ was connected to bus
$j$, $\bar{v}_{i}$ would be replaced by $\bar{v}_{j}$. If a component is
connected to multiple buses, all of the corresponding voltages will be
included in $h_{i}(\cdot)$ and $\bm{f}_{i}\left(\cdot\right)$ in 4b and 4a..
The update function $\bm{f}{}_{i}$ can furthermore depend on control inputs
$\bm{u}_{i}$. For static components the current injections becomes a function
of the local voltage $\bar{i}_{i}^{C}=h_{i}(\bar{v}_{i})$. If no component is
connected to bus $i$, then $\bar{i}_{i}^{C}=0$. If multiple components were
connected to the same bus, the current injections would be summed up.

Figure 1: Structure of the DAEs that govern the power system dynamics with
current injection $\bar{\bm{i}}^{C}$ from the components, e.g., generators and
loads, and current flows in the network $\bar{\bm{i}}^{N}$. Adapted from [15].

Figure 1 schematically depicts this structure with $m$ dynamic components, $l$
static component, and $n$ buses. When we formulate the current balance

$\displaystyle\bm{0}=\bar{\bm{i}}^{C}{}(\bm{x}(t),\bar{\bm{v}}(t))-\bar{\bm{i}}^{N}{}(\bar{\bm{v}}(t)),$
(5)

and consider the system’s structure, we note that the component currents
$\bar{\bm{i}}^{C}$ have only a local dependency, i.e., to their corresponding
state vectors $\bm{x}_{i}$ and the local voltage $\bar{v}_{i}$. Their
interaction arises from the coupling through the network currents
$\bar{\bm{i}}^{N}$, which in turn only depends on all complex voltages
$\bar{\bm{v}}$.

PINNSim exploits this structure as follows: We assume a voltage evolution
$\hat{\bar{\bm{v}}}{}(t)$ – we denote all approximations with the hat symbol
$\hat{}$ – based on which we evaluate the resulting network currents
$\hat{\bar{\bm{i}}}^{N}$ and the component currents $\hat{\bar{\bm{i}}}^{C}$.
For the latter, we furthermore require the state evolution $\hat{\bm{x}}_{i}$
of each component; we approximate the necessary temporal integration with
PINNs. When comparing the resulting currents

$\displaystyle\hat{\bar{\bm{i}}}^{C}{}\left(\hat{\bm{x}}(t),\hat{\bar{\bm{v}}}(t)\right)-\hat{\bar{\bm{i}}}^{N}{}\left(\hat{\bar{\bm{v}}}(t)\right)$
(6)

the result will not equal $\bm{0}$ by default, i.e., the current balance 5 is
not obeyed. However, the mismatch indicates the quality of the approximation.
By adjusting the assumed voltage evolution $\hat{\bar{\bm{v}}}{}(t)$, we
iteratively reduce this mismatch. The resulting voltage
$\hat{\bar{\bm{v}}}{}(t)$ and state $\hat{\bm{x}}{}(t)$ evolutions constitute
the solution of PINNSim to the system of DAEs.

###  III-B Parametrisation of the voltage evolution $\hat{\bar{v}}_{i}(t)$

We describe the evolution of the complex voltage in polar form
$\bar{v}_{i}(t)=V_{i}(t)e^{j\theta_{i}(t)}$. To approximate $\bar{v}_{i}(t)$,
we express the voltage magnitude $V_{i}$ and the voltage angle $\theta_{i}$ as
a power-series with respect to time up to power $r$ as suggested in [16]

$\displaystyle\hat{\bar{v}}_{i}(t)$
$\displaystyle=\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)}.$
(7)

The coefficients $V_{0,i},V_{1,i},\dots,V_{r,i}$ and
$\theta_{0,i},\theta_{1,i},\dots,\theta_{r,i}$ form the parameters which we
will later on update to improve the approximation. The vector $\bm{\Xi}{}_{i}$
collects all parameters at bus $i$

$\displaystyle\bm{\Xi}{}_{i}=\begin{bmatrix}V_{0,i}&\theta_{0,i}&\ldots&V_{r,i},\theta_{r,i}\end{bmatrix},\quad\bm{\Xi}{}_{i}\in\mathbb{R}^{2(r+1)}.$
(8)

This parametisation is repeated for all $n$ buses in the system. We
concatenate all $\bm{\Xi}{}_{i}$ in the vector
$\bm{\Xi}{}\in\mathbb{R}^{2(r+1)n}$.

###  III-C Solving component dynamics with PINNs

The exact solution for the evolution of the differential variables
$\bm{x}_{i}(t)$ of a single dynamic component $i$ can be obtained by
integration of 4b

$\displaystyle\bm{x}_{i}(t)$
$\displaystyle=\bm{x}_{0,i}+\int_{t_{0}}^{t}\bm{f}_{i}\left(\bm{x}_{i}(\tau),\bar{v}_{i}(\tau),\bm{u}_{i}\right)d\tau.$
(9)

As there usually exists no explicit analytical solution to 9, we need to
approximate the solution. As motivated earlier, using PINNs allows us to
accurately approximate 9 over large time steps $\Delta t=t-t_{0}$ while being
fast to evaluate. The approximation will be dependent on the time step size
$\Delta t$, the initial state $\bm{x}_{0,i}$, potentially the control input
$\bm{u}_{i}$, and the voltage profile $\bar{v}_{i}$. As we cannot train the
PINN with an arbitrary $\bar{v}_{i}$, we restrict the learning to the same
form as in 7, i.e., $\bar{v}_{i}\approx\hat{\bar{v}}_{i}(\Delta
t,\bm{\Xi}_{i})$. Hence, the approximation of 9 will depend also on
$\bm{\Xi}_{i}$

$\displaystyle\hat{\bm{x}}_{i}(t)=\operatorname{PINN}(\Delta
t,\bm{x}_{0,i},\bm{\Xi}_{i},\bm{u}_{i}).$ (10)

The functional form of the PINN is a standard feed-forward NN with $K$ hidden
layers parameterised by the weight matrices $\bm{W}^{k}$ and bias vectors
$\bm{b}^{k}$ and a non-linear activation function $\sigma$

$\displaystyle\bm{z}_{0}$ $\displaystyle=[\Delta
t,\bm{x}_{0,i},\bm{\Xi}_{i},\bm{u}_{i}]$ (11a) $\displaystyle\bm{z}_{k+1}$
$\displaystyle=\sigma\left(\bm{W}^{k+1}\bm{z}_{k}+\bm{b}^{k+1}\right),\,\forall
k=0,1,\ldots,K-1$ (11b) $\displaystyle\hat{\bm{x}}_{i}$
$\displaystyle=\bm{x}_{0,i}+\Delta t(\bm{W}^{K}\bm{z}_{K}+\bm{b}^{K}).$ (11c)

We adjust the last layer to enforce that the initial condition $\bm{x}_{0,i}$
is met if $\Delta t=0$, i.e., $t=t_{0}$. Thereby we ensure numerical
consistency and can improve inaccuracies related to recursive application of
the PINN [5]. To train a PINN, we optimise

$\displaystyle\min_{\bm{W}^{1},\bm{b}^{1},\dots,\bm{W}^{K},\bm{b}^{K}}\quad\mathcal{L}_{x}+\alpha\mathcal{L}_{c},$
(12)

where $\mathcal{L}_{x}$ evaluates the prediction error based on a dataset
$\mathcal{D}$ of simulated data points

$\displaystyle\mathcal{L}_{x}=\frac{1}{\left\lvert\mathcal{D}\right\rvert}\sum_{j=1}^{\left\lvert\mathcal{D}\right\rvert}\left\lVert\bm{x}_{i}^{(j)}-\hat{\bm{x}}_{i}^{(j)}\right\rVert_{2}^{2}$
(13)

and a physics-based loss term $\mathcal{L}_{c}$ evaluated on collocation
points which do not require any simulation

$\displaystyle\mathcal{L}_{c}=$
$\displaystyle\frac{1}{\left\lvert\mathcal{D}\right\rvert}\sum_{j=1}^{\left\lvert\mathcal{D}_{c}\right\rvert}\left\lVert\frac{d}{dt}\hat{\bm{x}}_{i}^{(j)}-\bm{f}\left(\hat{\bm{x}}_{i}^{(j)},\bar{v}_{i}^{(j)},\bm{u}_{i}^{(j)}\right)\right\rVert_{2}^{2}.$
(14)

The two loss terms are weighted with the hyperparameter $\alpha$. For a more
detailed explanation of PINNs, we refer to [9]. Based on the evolution of the
state $\hat{\bm{x}}_{i}(t)$ and the voltage $\hat{\bar{v}}_{i}(t)$, we can
calculate the current injection $\hat{\bar{i}}_{i}^{C}$ by evaluating

$\displaystyle\hat{\bar{i}}_{i}^{C}(t)=h_{i}(\hat{\bm{x}}_{i}(t),\hat{\bar{v}}_{i}(t)).$
(15)

For each of the $m$ dynamic components, we will train a separate PINN and all
training can be can be performed in advance of executing PINNSim.

###  III-D Updating the voltage profile

The approximations $\hat{\bar{v}}_{i}(t)$ and $\hat{\bm{x}}_{i}(t)$ are both
continuous functions and therefore the same holds for
$\hat{\bar{i}}_{i}^{C}(t)$ and $\hat{\bar{i}}_{i}^{N}(t)$. The current balance
5 will only be fulfilled, if
$\hat{\bar{i}}_{i}^{C}(t)=\hat{\bar{i}}_{i}^{N}(t)$ for the entire time step
and at all buses. We express this notion by requiring that the norm
$\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert$ shall be
0. This norm

$\displaystyle\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert=\sqrt{\langle\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N},\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\rangle}$
(16)

is induced by the inner product $\langle\bar{a},\bar{b}\rangle$ between two
complex functions $\bar{a}(t),\bar{b}(t)$ over the interval
$[t_{0},t_{0}+\Delta t]$

$\displaystyle\langle\bar{a},\bar{b}\rangle$
$\displaystyle=\int_{t_{0}}^{t_{0}+\Delta t}\bar{a}(t)\,\bar{b}^{*}(t)\,dt.$
(17)

As $\hat{\bar{i}}_{i}^{C}$ and $\hat{\bar{i}}_{i}^{N}$ depend on the
parametrised voltages $\hat{\bar{\bm{v}}}{}(t,\bm{\Xi})$, we subsequently aim
to find a parametrisation $\bm{\Xi}$ that minimises
$\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert$. To this
end, we formulate the following optimisation problem for the entire system as
the sum of the norms (we square 16 to avoid the calculation of the square
root)

$\displaystyle\min_{\bm{\Xi}}\sum_{i=1}^{n}\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert^{2}.$
(18)

The above expression involves the integration

$\displaystyle\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert^{2}$
$\displaystyle=\int_{t_{0}}^{t_{0}+\Delta
t}\Re\left(\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right)^{2}+\Im\left(\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right)^{2}dt$

which we approximate with the Midpoint rule, i.e., we split the interval
$[t_{0},t_{0}+\Delta t]$ into $s$ equally sized intervals of width
$\frac{\Delta t}{s}$ and sum the function value at the middle of these
intervals

$\displaystyle\begin{split}\left\lVert\hat{\bar{i}}_{i}^{C}-\hat{\bar{i}}_{i}^{N}\right\rVert^{2}\approx\frac{\Delta
t}{s}\sum_{j=1}^{s}&\Re\left(\hat{\bar{i}}_{i}^{C}(t_{j})-\hat{\bar{i}}_{i}^{N}(t_{j})\right)^{2}\\\
&+\Im\left(\hat{\bar{i}}_{i}^{C}(t_{j})-\hat{\bar{i}}_{i}^{N}(t_{j})\right)^{2}.\end{split}$
(19)

Thereby, we can approximate the optimisation in 18 as

$\displaystyle\min_{\bm{\Xi}}\frac{\Delta t}{s}\bm{\rho}{}^{\top}\bm{\rho}{}$
(20)

where $\bm{\rho}$ collects all summands in a vector

$\displaystyle\begin{split}\bm{\rho}{}&=\begin{bmatrix}\Re\left(\hat{\bar{i}}_{1}^{C}(t_{1})-\hat{\bar{i}}_{1}^{N}(t_{1})\right)\\\
\Im\left(\hat{\bar{i}}_{1}^{C}(t_{1})-\hat{\bar{i}}_{1}^{N}(t_{1})\right)\\\
\vdots\\\
\Re\left(\hat{\bar{i}}_{n}^{C}(t_{s})-\hat{\bar{i}}_{n}^{N}(t_{s})\right)\\\
\Im\left(\hat{\bar{i}}_{n}^{C}(t_{s})-\hat{\bar{i}}_{n}^{N}(t_{s})\right)\end{bmatrix},\quad\bm{\rho}{}\in\mathbb{R}^{2ns}\end{split}$
(21)

We solve the above non-linear least square problem by iteratively updating the
parameters $\bm{\Xi}^{(k+1)}=\bm{\Xi}^{(k)}+\Delta\bm{\Xi}$. To determine
$\Delta\bm{\Xi}$, we compute the Jacobian

$\displaystyle\bm{J}$
$\displaystyle=\frac{\partial\bm{\rho}{}}{\partial\bm{\Xi}}\Bigr{|}_{\bm{\Xi}^{(k)}},\quad\bm{J}\in\mathbb{R}^{2ns\times
2(r+1)n}$ (22)

at the values $\bm{\Xi}^{(k)}$ and solve the linear problem

$\displaystyle\left(\bm{J}^{\top}\bm{J}\right)\Delta\bm{\Xi}=-\bm{J}^{\top}\bm{\rho}{}.$
(23)

To calculate $\bm{J}$ in 22, we use Automatic Differentiation (AD) [17] as all
computations for $\hat{\bar{\bm{i}}}^{C}$ and $\hat{\bar{\bm{i}}}^{N}$ are
expressed as explicit functions. While the size of $\bm{J}$ can become large,
it has a very sparse structure that can be exploited in its construction and
when solving 23 – herein lies the key to the scalability of PINNSim as it
closely resembles the structure of conventional integration schemes such as
the trapezoidal method.

###  III-E PINNSim: The full time-stepping simulator

Algorithm 1 shows the integration of the previous sections into the full
PINNSim algorithm that is used for the computation of a time step.

Algorithm 1 PINNSim - single time step

Require: $\bm{x}_{0}$, $t_{0}$, $t_{\max}$, $\Delta t$, $s$, $r$,
$\Delta\bm{\Xi}^{\max}$, $k^{\max}$
Initialise: $k=0,\bm{\Xi}^{(0)},\Delta\bm{\Xi}^{(0)}$

1:while $\Delta\bm{\Xi}^{(k)}>\Delta\bm{\Xi}^{\max}$ and $k<k^{\max}$ do

2: for component $i=1,\dots,m$ & query points $j=1,..,s$ do

3: Predict state with PINN
$\hat{\bm{x}}_{i}(t_{j},\bm{x}_{0,i},\bm{\Xi}_{i}^{(k)})$

4: Calculate component injections
$\hat{\bar{i}}_{i}^{C}=h_{i}(t_{j},\hat{\bm{x}}_{i},\bm{\Xi}_{i}^{(k)})$

5: Calculate contribution of $\hat{\bar{i}}_{i}^{C}$ to $\bm{J}$ end

6: for query points $j=1,\dots,s$ do

7: Calculate network injections $\hat{\bar{\bm{i}}}^{N}(t_{j},\bm{\Xi}^{(k)})$

8: Calculate contribution of $\hat{\bar{\bm{i}}}^{N}$ to $\bm{J}$ end

9: Assemble $\bm{J}$ and solve
$\left(\bm{J}^{\top}\bm{J}\right)\Delta\bm{\Xi}^{(k+1)}=-\bm{J}\bm{\rho}$

10: Update iteration $\bm{\Xi}^{(k+1)}=\bm{\Xi}^{(k)}+\Delta\bm{\Xi}^{(k)}$,
$k=k+1$ end

11:return Trajectory across time step
$t,\hat{\bm{x}}{}(t),\hat{\bar{\bm{v}}}{}(t)$ with $\bm{\Xi}^{(\text{final})}$

Its accuracy is dependent on the approximation quality of $\hat{\bm{x}}{}(t)$
and $\hat{\bar{\bm{v}}}(t)$ and on the tolerance settings for $\Delta\bm{\Xi}$
and the maximum number of iterations $k^{\max}$. Therefore, requirements on
the resulting tolerance limit the suitable time step size $\Delta t$. By
repeatedly applying Algorithm 1, we obtain a time-stepping scheme that then
allows the simulation of dynamics beyond $\Delta t$.

##  IV Case study

This section describes the modelling of the power system dynamics and the
implementation of PINNSim.

###  IV-A Power system modelling

As an example for a dynamic component, we consider a two-axis generator model
as modelled in [15].

$\displaystyle\begin{bmatrix}\scriptstyle T^{\prime}_{do}\\\ \scriptstyle
T^{\prime}_{qo}\\\ \scriptstyle 1\\\ \scriptstyle 2H\end{bmatrix}\frac{d}{dt}$
$\displaystyle\begin{bmatrix}\scriptstyle E^{\prime}_{q}\\\ \scriptstyle
E^{\prime}_{d}\\\ \scriptstyle\delta\\\
\scriptstyle\Delta\omega\end{bmatrix}=\begin{bmatrix}\scriptstyle-E^{\prime}_{q}-(X_{d}-X^{\prime}_{d})I_{d}+E_{fd}\\\
\scriptstyle-E^{\prime}_{d}+(X_{q}-X^{\prime}_{q})I_{q}\\\
\scriptstyle\omega_{s}\Delta\omega\\\ \scriptstyle
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}$
(24a) $\displaystyle\begin{bmatrix}I_{d}\\\ I_{q}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}R_{s}&-X^{\prime}_{q}\\\
X^{\prime}_{d}&R_{s}\end{bmatrix}^{-1}\begin{bmatrix}E^{\prime}_{d}-V\sin{(\delta-\theta)}\\\
E^{\prime}_{q}-V\cos{(\delta-\theta)}\end{bmatrix}$ (24b)
$\displaystyle\bar{i}^{C}$
$\displaystyle=(I_{D}+jI_{Q})=(I_{d}+jI_{q})e^{j(\delta-\pi/2)}.$ (24c)

24a corresponds to 4b and 24b and 24c to 4a. For this study, we simplify the
model above to a classical machine 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 [15]. Now, the rotor angle $\delta$ and the
frequency deviation $\Delta\omega$ form the state $\bm{x}_{i}$, the magnitude
$V$ and angle $\theta$ of the terminal voltage form $\bar{v}_{i}$, and the
mechanical power $P_{m}$ and the excitation voltage $E_{fd}$ form the control
input $\bm{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.

To illustrate PINNSim, we consider the IEEE 9-bus system described in [15, pp.
164–167] with three generators (all modelled as classical machines as above)
with parameters from Table I. The initial conditions $\bm{x}_{0,i}$ and
control inputs $\bm{u}$ are determined from assuming an equilibrium state for
the load flow case in [15]. To perturb the system, we reduce the mechanical
power $P_{m}$ of generator 1 to 50% of its initial value and then observe the
resulting trajectory.

TABLE I: Generator parameters and set points (all in $\mathrm{p.u.}$) Gen. | $H$ | $D$ | $X_{d}$ | $X^{\prime}_{d}$ | $R_{s}$ | $P_{m}$ | $E_{fd}$
---|---|---|---|---|---|---|---
1 | 23.64 | 2.364 | 0.146 | 0.0608 | 0.0 | 0.71 | 1.08
2 | 6.4 | 1.28 | 0.8958 | 0.1969 | 0.0 | 1.612 | 1.32
3 | 3.01 | 0.903 | 1.3125 | 0.1813 | 0.0 | 0.859 | 1.04

###  IV-B Implementation and NN training

The entire simulator is implemented in PyTorch [18] as we utilise AD for the
computation of the Jacobian $\bm{J}$ and it enables the training of the PINN
models. For each generator we used 2500 simulated points in the training
dataset $\mathcal{D}$ and 5000 collocation points. Each PINN consists of 2
hidden layers with 32 neurons, uses a $\tanh{}$ activation function, and is
trained for a range of possible voltage parametrisations $\bm{\Xi}_{i}$ with
$r=2$, initial conditions $\bm{x}_{0,i}$, and prediction time steps $\Delta
t\in[0,0.3]$\mathrm{s}$$, which results in a nine-dimensional input space. We
train, using a L-BFGS optimiser, for 2000 epochs. The published code base [19]
provides a detailed overview of the dataset creation setup and the training
procedure. The numerical simulations of the system of DAEs that serve as
ground truth are implemented using Assimulo [20]. As a comparison to PINNSim,
we implement a trapezoidal integration scheme, see [21].

##  V Results

The potential for accelerating TDSs with PINNSim relies on increasing the time
step size and thereby reducing the total number of time steps compared to
established methods such as the trapezoidal method. We demonstrate this
behaviour in Section V-A on the test case. In Section V-B, we then illustrate
what allows PINNSim these larger time steps and describe in Section V-C how
the involved computational cost compare.

###  V-A Accuracy of PINNSim for a full TDS

First, we consider a simulation of the described test case over
$2.5\text{\,}\mathrm{s}$. Figure 2 shows in the upper panels the resulting
trajectory of the frequency deviation at machine 2, i.e., $\Delta\omega_{2}$.
While PINNSim and the trapezoidal rule accurately capture the dynamics for a
time step size of $\Delta t=$0.05\text{\,}\mathrm{s}$$, the trapezoidal rule
fails to track the evolution for a larger time step size of $\Delta
t=$0.25\text{\,}\mathrm{s}$$. In contrast, PINNSim captures the state
evolution accurately.

We test this behaviour for more time step sizes and plot the maximum error
over this time interval in the lower panel in Fig. 2. The trapezoidal rule
performs well for very small time step sizes, but the errors quickly become
very large - for reference, the overall variation of $\Delta\omega_{2}$ is
around $0.8\text{\,}\mathrm{Hz}$. PINNSim leads to an almost constant error
characteristics, mostly with less than $0.01\text{\,}\mathrm{Hz}$ error. The
slightly increasing error towards smaller time step sizes arises due to
accumulating errors. The larger errors for $\Delta t>$0.3\text{\,}\mathrm{s}$$
are expected, as we trained the PINNs only up this time step size.

Figure 2: The upper panels show the simulated trajectories of machine 2 with a
trapezoidal and PINNSim time stepping schemes for two time step sizes. The low
panel shows the maximum error for $\Delta\omega_{2}$ over the trajectory for a
range of time step sizes.

###  V-B Accuracy of PINNSim on a single time step

In the following, we consider a single time step and show the influencing
factors on the accuracy. To obtain a range of system states, we simulate the
test case for $10\text{\,}\mathrm{s}$ with an accurate solver and then extract
200 instances, i.e., every $0.05\text{\,}\mathrm{s}$.

####  V-B1 Voltage parametrisation

We use these 200 instances as initial values to test the performance of
PINNSim on a single time step for different values of $\Delta t$, the voltage
profile order $r$ and the number of query points $s$.

Figure 3: Maximum error for single time step predictions of different length
$\Delta t$ starting from 200 different initial conditions. Comparison between
the trapezoidal rule and PINNSim with different voltage scheme orders ($r$)
and number of query points $s$. Larger $r$ and $s$, both improve the accuracy
of PINNSim.

We evaluate these time steps with respect to the maximum error of
$\Delta\omega_{2}$ at the end of the time step. Figure 3 reports the results.
A main observation is that PINNSim with $r=1$ (linear voltage evolution) and
$s=r+1$ has a comparable relationship between accuracy and time step size as
the trapezoidal rule. However, when using PINNSim with $r=2$, the time step
size can be significantly increased without incurring much larger errors. For
very small time step sizes PINNSim offers no benefits over the trapezoidal
rule as the error decreases slower when reducing $\Delta t$. We can improve
this characteristic for PINNSim by evaluating 19 on more points, i.e., choose
a larger value for $s$, shown as the dashed lines. The improved performance,
in particular for $r=1$, originates from the better approximation of 18.
However, for large time steps, here for around $\Delta
t>$0.2\text{\,}\mathrm{s}$$, we additionally require $r=2$ for accurate
predictions; the voltage evolution becomes too non-linear, hence, the linear
voltage approximation ($r=1$) is insufficient. For $r=2$, the limitation in
time step size arises as training domain of the PINNs was restricted to
$\Delta t\leq$0.3\text{\,}\mathrm{s}$$.

####  V-B2 PINN accuracy

The performance of PINNSim relies in large parts on how well each PINN
approximates the integration in 9 for a given voltage profile
$\hat{\bar{v}}_{i}(t,\bm{\Xi}_{i})$. Figure 4 illustrates the results on a
test dataset of 4000 points for generator 2. Both panels show the same
results, i.e., the maximum and median error on the test dataset, but on a
logarithmic and linear x-axis. The logarithmic axis clearly shows that for
small time step sizes, the error of PINNs decreases as we included a time
dependency in the final layer 11c. The linear x-axis gives a more intuitive
understanding of how much further PINNs can predict the dynamics accurately.
Only when reaching the limit of the training domain, the accuracy
deteriorates. Revisiting the error plots in Figs. 2 and 3, we can find that it
is the particular error characteristic of PINNs that allows the large time
steps of PINNSim and why it outperforms the trapezoidal method.

Figure 4: Error characteristics of the PINN for machine 2 on a test dataset.
Both plots show the same results but on a logarithmic x-axis (left) and a
linear x-axis (right) to highlight the accuracy of PINNs over large time
steps.

###  V-C Computation cost of a time step

The larger time step sizes of PINNSim come at increasd computational cost
compared to the trapezoidal rule. Providing a fair comparison based on
optimised implementations will be future work, instead, we will here focus on
the main cost of PINNSim from a perspective of scalability.

In Algorithm 1, the evaluation of the inner loop (lines 2-8) requires the
calculation of network and component current injections
$\hat{\bar{\bm{i}}}^{N}$, $\hat{\bar{\bm{i}}}^{C}$ and their sensitivity to
$\bm{\Xi}$. For the network currents, this requires merely matrix-vector
products and for the component currents the evaluations of the PINNs. A single
evaluation (pass) of a PINN requires on the order of
$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$, e.g., to compute the current
injection, and the sensitivity calculations requires only two additional
passes thanks to AD. As all PINNs could be evaluated in parallel, these
computations can easily be scaled to large system. The expensive computation
arises in the solution of the linear system of equations in line 9 of
Algorithm 1, but its sparsity and known structure can be exploited. A closely
related calculation is required for the trapezoidal rule.

The overall cost then scales linearly with the number of iterations $k$ per
time step, i.e., line 1-10 in Algorithm 1, until convergence is reached. We
show in Fig. 5 the value of the objective 19 over the iterations for two time
step sizes and different values of $r$. In the left panel, the least square
problem is fully determined ($s=r+1$), hence the objective value continues to
decrease. In the right panel, the five additional query points render the
problem over-determined, hence the algorithm converges at a non-zero objective
value. In either case, we observe that larger time step sizes require more
iterations, while the order of the voltage profiles primarily affects the
magnitude of the objective value and less its convergence.

Figure 5: Convergence of the objective value in 19 for different number of
query points $s$, voltage profile orders $r$ and time step sizes $\Delta t$.

As stated in [9, 12], the cost of training the PINNs has to be considered as
well. Improving the efficiency of the training, however, can be decoupled from
the analysis of PINNSim.

##  VI Discussion

The presented results shall serve as a proof of concept of PINNSim as a novel
time stepping simulator. We demonstrated that PINNSim outperforms the
trapezoidal rule on the metric of the allowable time step size. The following
describes four aspects that are necessary to turn PINNSim from a proof of
concept into a fully fledged simulator for power system dynamics and to
realise the potential acceleration of TDSs.

  1. 1.

Speed: Numerical methods require a high level of optimisation in the
implementation to become competitive. For PINNSim this concerns primarily the
optimisation of the calculation of the residual $\bm{\rho}$, the Jacobian
$\bm{J}$, and the solution of the sparse linear system in 23. The focus should
therefore lie on controlling memory allocations, utilising parallelisation,
and exploiting sparsity patterns for these computations.

  2. 2.

Accuracy: The accuracy of PINNSim hinges around accurately learned PINNs. We
envision that this process can be highly standardised, so that PINNs can be
trained reliably to high accuracy and with desirable error characteristics for
a wide range of dynamic components.

  3. 3.

Scale: The analysis of the computational complexity of PINNSim suggests its
scalability but it remains to be shown in practice. Furthermore, we need to
investigate the accuracy and convergence properties of the voltage update
scheme for larger systems.

  4. 4.

Applications: In this work, we have considered the dynamics of classical
machine models under a set-point change. Future work should include more
detailed dynamical models and apply PINNSim to short-circuits, topology
changes, and other numerically more demanding setups.

##  VII Conclusion

With PINNSim, we introduced a novel approach for time-domain simulations that
allows the integration of PINNs. As a result we achieve accurate results with
significantly larger time steps than with the common trapezoidal method.
Hence, we require fewer time steps which could result in a significant
acceleration of time-domain simulations in power systems. By design, we only
require the training of PINNs for single dynamical components and thereby
enable the scalability of PINNSim to large systems. Additionally, all
calculations at run-time are highly scalable and parallelisable, necessary
requirements to develop PINNSim into a powerful simulator.

## References

  * [1] B. Stott, “Power system dynamic response calculations,” _Proceedings of the IEEE_ , vol. 67, no. 2, pp. 219–241, Feb. 1979.
  * [2] 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.
  * [3] G. Gurrala, A. Dimitrovski, S. Pannala, S. Simunovic, and M. Starke, “Parareal in Time for Fast Power System Dynamic Simulations,” _IEEE Transactions on Power Systems_ , vol. 31, no. 3, pp. 1820–1830, May 2016.
  * [4] P. Aristidou, “Time-domain simulation of large electric power systems using domain-decomposition and parallel processing methods,” Ph.D. dissertation, Université de Liège, Liège, Belgium, 2015.
  * [5] I. Lagaris, A. Likas, and D. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” _IEEE Transactions on Neural Networks_ , vol. 9, no. 5, pp. 987–1000, Sep. 1998.
  * [6] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, no. C, Nov. 2018.
  * [7] G. S. Misyris, A. Venzke, and S. Chatzivasileiadis, “Physics-Informed Neural Networks for Power Systems,” in _2020 IEEE Power & Energy Society General Meeting (PESGM)_. Montreal, QC, Canada: IEEE, 2020, pp. 1–5.
  * [8] C. Moya, G. Lin, T. Zhao, and M. Yue, “On Approximating the Dynamic Response of Synchronous Generators via Operator Learning: A Step Towards Building Deep Operator-based Power Grid Simulators,” Jan. 2023, arXiv:2301.12538 [cs.LG].
  * [9] J. Stiasny and S. Chatzivasileiadis, “Physics-informed neural networks for time-domain simulations: Accuracy, computational cost, and flexibility,” _Electric Power Systems Research_ , vol. 224, p. 109748, Nov. 2023.
  * [10] W. Cui, W. Yang, and B. Zhang, “Predicting Power System Dynamics and Transients: A Frequency Domain Approach,” Nov. 2021, arXiv: 2111.01103 [eess.SY].
  * [11] C. Moya and G. Lin, “DAE-PINN: a physics-informed neural network model for simulating differential algebraic equations with application to power networks,” _Neural Computing and Applications_ , vol. 35, no. 5, pp. 3789–3804, Feb. 2023.
  * [12] J. Stiasny, “Physics-Informed Neural Networks for Power System Dynamics,” Ph.D. dissertation, Technical University of Denmark, 2023.
  * [13] K. E. Brenan, S. L. Campbell, and L. R. Petzold, _Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations_. Society for Industrial and Applied Mathematics, Jan. 1995.
  * [14] G. Adomian, “A review of the decomposition method in applied mathematics,” _Journal of Mathematical Analysis and Applications_ , vol. 135, no. 2, pp. 501–544, Nov. 1988.
  * [15] P. W. Sauer and M. A. Pai, _Power system dynamics and stability_. Upper Saddle River, N.J: Prentice Hall, 1998.
  * [16] B. Wang, N. Duan, and K. Sun, “A Time–Power Series-Based Semi-Analytical Approach for Power System Simulation,” _IEEE Transactions on Power Systems_ , vol. 34, no. 2, pp. 841–851, Mar. 2019.
  * [17] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, “Automatic Differentiation in Machine Learning: a Survey,” _Journal of Machine Learning Research_ , vol. 18, no. 153, pp. 1–43, 2018.
  * [18] A. Paszke _et al._ , “PyTorch: An imperative style, high-performance deep learning library,” in _Advances in neural information processing systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019.
  * [19] J. Stiasny, “Publicly available implementation,” 2023. [Online]. Available: https://github.com/jbesty
  * [20] C. Andersson, C. Führer, and J. Åkesson, “Assimulo: A unified framework for ODE solvers,” _Mathematics and Computers in Simulation_ , vol. 116, pp. 26–43, Oct. 2015.
  * [21] F. Milano, _Power System Modelling and Scripting_ , ser. Power Systems. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.