Scientific Machine Learning for Power System Dynamic

Simulation

by

Matthew. A. Bossart

B.E., Dartmouth College, 2015

M.S., University of Colorado Boulder, 2023

A thesis submitted to the

Faculty of the Graduate School of the

University of Colorado in partial fulfillment

of the requirements for the degree of

Doctor of Philosophy

Department of Electrical, Computer, and Energy Engineering

2024

Committee Members:

Bri-Mathias Hodge, Chair

Dragan Maksimović

Kyri Baker

Eduardo Corona

Duncan Callaway
Bossart, Matthew. A. (Ph.D., Electrical, Computer, and Energy Engineering)

Scientific Machine Learning for Power System Dynamic Simulation

Thesis directed by Prof. Bri-Mathias Hodge

It is imperative to decarbonize the power system as quickly as possible to respond to the threat

of climate change. In pursuit of this goal, power systems around the world are undergoing a period

of dramatic change. Large, centralized, fossil-fuel based generators are being replaced by distributed

renewables such as wind and solar photovoltaics. While thermal generators interface to the grid via

synchronous machines, the fastest growing renewables are connected to the grid via power electronic

inverters. This change alone—from synchronous machines to inverters—fundamentally alters the

dynamic behavior of the power system at the sub-second to minutes time scales. Furthermore, this

change impacts the tools that are used to study and operate the power system at fast time scales.

The status quo for studying large-signal dynamic behavior is numerical simulation. In this area,

simplified formulations that have been used historically cannot reliably capture the true behavior

of an inverter-dominated system. On the other hand, more accurate formulations were designed

for small studies and cannot scale to the largest systems. The power systems community needs

additional tools and methods for simulating the power system that can scale computationally while

capturing the behavior of an inverter-dominated system.

This thesis proposes methods for incorporating machine learning or data-driven models into

the existing scientific computing workflow (i.e. scientific machine learning). This hybrid approach

takes advantage of the strengths of each paradigm; the expressiveness and speed of modern machine

learning architectures is combined with the interpretability of models built from first principles.

This thesis addresses practical barriers to incorporating data-driven models in order to improve the

trade-off between computational burden and model fidelity. Chapter 3 describes a novel data-driven

model architecture based on deep equilibrium layers and neural ordinary differential equations

(ODEs). While neural ODEs can be naturally combined with the differential equations describ-
 iii

ing traditional power system models, the proposed model additionally satisfies the initialization

requirements of a power system device within the simulation workflow. Unlike the neural ODE

model, which is used in conjunction with a numerical solver, other proposed data-driven models

predict directly the state values at a future point in time. Chapter 4 proposes a method for incorpo-

rating such models through the inclusion of delays in the system model. One of the key additional

challenges in adopting any data-driven model is the need for training prior to deployment. Chapter

5 proposes a tool for power system simulation that is end-to-end differentiable with respect to pa-

rameters. The result is the ability to train the models proposed in previous chapters directly within

a simulation framework, thereby simplifying the process of training and using data-driven models,

regardless of the model architecture. Using data-driven models in power systems simulations is

promising but not fully explored; this thesis develops the tools and methodologies for incorporating

such data-driven models and sets the stage for continued development of models so that they can

be adopted into industrial practice.

 Dedication

For my parents.
 v

Acknowledgements

First, thank you to my advisor Professor Bri-Mathias Hodge for the unwavering support.

It was a true comfort to know that you always had my best interest in mind both personally

and professionally. Thank you to Professor Dragan Maksimović for accepting me as a student

and mentoring me in my first year and to the rest of my thesis committee—Kyri Baker, Eduardo

Corona, Duncan Callaway, and Spyros Chatzivasileiadis—for their time and effort. Thank you to

my colleagues from the Hodge lab, particularly Marija Marković, Wallace Kenyon, Amir Sajadi,

Muhy Eddin Za’tar, Marena Trujillo, and Fiona Majeau for the friendship, support, and fruitful

conversations. Thank you to the entire ECEE department, especially Patty Sue Femrite, Lori

Mehan and Audrey Cook; my interactions with you brightened my day on many occasions.

Thank you to my collaborators at UC Berkeley: Jose Daniel Lara, Rodrigo Henriquez-Auba,

and Ciaran Roberts. Working with you on the AGM project was a formative experience of my PhD.

One of the highlights of my PhD was my research visit to ETH. Thank you to Professor Gabriella

Hug and the students of the Power Systems Laboratory for welcoming me to your lab, to Professor

Spyros Chatzivasileiadis for advising me, and to the Marconi family for graciously hosting me in

your home for four months. I will always cherish my time in Zurich.

Thank you to my family for your continued love and support. To my parents, for supporting

my education and instilling in me the idea that a career can and should be enjoyable.

Lastly, an enormous thank you to Sean. I’m so thankful to have had you by my side for this

entire journey, you provided the support through the ups and downs and words of encouragement

that I often needed.

 vi

Contents

Chapter

1 Introduction 1

1.1 A Changing Power System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Power System Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Transient Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4 Research Impact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.4.1 Contingency Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.4.2 Interconnection Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4.3 Validation of Novel Controls . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2 Preliminaries and Literature Review 13

2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.1 Power System Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.2 Scientific Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1.3 Machine Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.1.4 The Intersections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.2 State of Practice for Power System Transient Simulations . . . . . . . . . . . . . . . 19

2.2.1 Waveform Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.2.2 Positive Sequence Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.3 [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
 vii

2.4 Methodology Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.5 Numerical Test Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.6 Literature Review and Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.6.1 Chapter 3: Relevant Literature . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.6.2 Chapter 4: Relevant Literature . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.6.3 Chapter 5: Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2.7 Summary of Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3 Acceleration of Power System Dynamic Simulations using a Deep Equilibrium Layer and

Neural ODE Surrogate 44

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.2 Surrogate Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.2.1 Power Systems Time Domain Simulation Procedure . . . . . . . . . . . . . . 49

3.2.2 Network Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.2.3 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.2.4 Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.3 Training Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.3.1 Data Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.3.2 Loss Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.3.3 Training Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.4.1 Computational Set Up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.4.2 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4 Integrating Data-Driven Dynamic Power System Models with Delay Differential Equations 72

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.1.1 Scope and Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

 viii

4.1.2 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

4.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.2.1 Power System Time-Domain Simulation . . . . . . . . . . . . . . . . . . . . . 77

4.2.2 System Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.2.3 Integrating Surrogate Models with Delay Equations . . . . . . . . . . . . . . 79

4.2.4 Solving Delay Differential Equations . . . . . . . . . . . . . . . . . . . . . . . 81

4.2.5 Computational Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.3.1 Test System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.3.2 Single Surrogate Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

4.3.3 Multiple Surrogate Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

4.3.4 Characterizing Computational Trade-offs . . . . . . . . . . . . . . . . . . . . 92

4.4 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5 Differentiable Programming for Transient Power System Simulations 98

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.2 Background and Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

5.3 Differentiating Power System Simulations . . . . . . . . . . . . . . . . . . . . . . . . 103

5.3.1 Nonlinear Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

5.3.2 Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.4 Software Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

5.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.5.1 Single Machine Infinite Bus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.5.2 Enhancing Understanding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

5.5.3 Parameter Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

5.5.4 Training Data-driven Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.5.5 Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

 ix

5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6 Conclusions and Future Work 122

6.1 Summary of Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

6.2 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

6.2.1 Robust Training of Data-driven Models for Power System Dynamic Simulations126

6.2.2 Capturing Uncertainty in Power System Transient Models . . . . . . . . . . . 127

6.2.3 Benchmarking Data-driven Models for Contingency Analysis . . . . . . . . . 128

6.2.4 Physics Informed Neural Networks for Transient Simulations . . . . . . . . . 129

Bibliography 130
 x

Tables

Table

3.1 Summary of dynamic models used in the 36-bus case study system . . . . . . . . . . 61

3.2 Summary of tunable hyperparameters considered for data-driven and physics-based

surrogates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.3 Comparison of model complexity in terms of number of states and number of pa-

rameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.1 Summary of differences between proposed direct-solution surrogate models . . . . . . 80

 xi

Figures

Figure

1.1 Global energy consumption by sector . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Global electric vehicle share by category . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Global electricity generation by source . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4 Overview of power system changes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.5 Model error vs. computational burden for power system transient simulations . . . . 7

1.6 Power system dynamic timescales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1 Venn diagram: the key topics of this thesis . . . . . . . . . . . . . . . . . . . . . . . 14

2.2 Categorization of ways to study a system . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3 [Link] inverter meta-model . . . . . . . . . . . . . . . . . . 24

2.4 [Link] software dependencies . . . . . . . . . . . . . . . . . 25

2.5 Methodology overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.6 A categorization of machine learning models: time-stepping vs. direct-solution models 33

3.1 Depiction of integrating data-driven surrogate models for computational benefits . . 46

3.2 Categorization of machine learning applications to power systems . . . . . . . . . . . 47

3.3 Overview of Chapter 3 methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.4 Block diagram of the neural ordinary differential equation and deep equilibrium layer

based surrogate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.5 Illustrative simulation results for the propose surrogate, both before and after training 55
 xii

3.6 Overview of the computational setup for generating results including software de-

pendencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.7 Diagram of the 36-bus case study system along with the two surrogate models . . . 60

3.8 The voltage and current trajectories of the training dataset . . . . . . . . . . . . . . 62

3.9 Mean average error comparison for various permutations of the physics-based surrogate 65

3.10 A sample trajectory result from the test dataset . . . . . . . . . . . . . . . . . . . . 67

3.11 Simulation time comparison for the full-order simulation and surrogate simulations . 68

3.12 Mean average error comparison for the data-driven and physics-based surrogate models 68

3.13 Mean average error comparison for an out-of-sample test dataset . . . . . . . . . . . 69

4.1 Replacing transient simulations with machine learning vs. integrating machine learn-

ing models into transient simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.2 Time-stepping models vs. direct-solution models . . . . . . . . . . . . . . . . . . . . 75

4.3 The assumed system structure for incorporating data-driven surrogate models . . . . 78

4.4 Splitting the system at a single bus of common coupling . . . . . . . . . . . . . . . . 79

4.5 Overview of the computational setup for generating results including software de-

pendencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.6 Diagram of the 2-area case study system . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.7 Depiction of replacing a full area with a single direct-solution surrogate . . . . . . . 88

4.8 Single surrogate: Comparison of current trajectories for the ground truth systems

and various surrogates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

4.9 Work-precision diagram for a delay differential equation system . . . . . . . . . . . . 90

4.10 Depiction of replacing individual generators with direct-solution surrogates . . . . . 91

4.11 Multiple surrogates: Comparison of current trajectories for the ground truth system

and various surrogates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

4.12 Comparison of simulation time between differential equation and delay differential

equation systems for different reductions in the number of states . . . . . . . . . . . 93

 xiii

4.13 Comparison of simulation time between differential equation and delay differential

equation systems for different lengths of the delay for the terminal data surrogate . . 94

4.14 Comparison of simulation time between differential equation and delay differential

equation systems for different number of parameters in the terminal data surrogate . 95

5.1 Block diagram depicting the steps of a power system transient simulation . . . . . . 102

5.2 The software dependencies for the forward numerical solution of a power system

simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

5.3 The software dependencies for the gradient computation of a power system simulation110

5.4 Code sample demonstrating the computation of the gradient of a simulation . . . . . 111

5.5 Illustrative example of the gradient with respect to various parameters over time for

a one machine infinite bus system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

5.6 Time-domain voltage trajectories and gradients for a maximum voltage deviation

objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

5.7 Time-domain frequency trajectories and gradients for a maximum frequency devia-

tion objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

5.8 Results of a parameter tuning problem where the inclusion of the initialization in

the gradient calculation is critical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

5.9 Ground truth and modified systems for a case study of incorporating machine learn-

ing surrogate models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.10 Comparison of the interfacing voltage and current trajectories for the ground truth

model, data-driven surrogate model, and common power system surrogates . . . . . 117

5.11 Mean average error comparison of various surrogates . . . . . . . . . . . . . . . . . . 118

5.12 Benchmarking results depicting the forward simulation time and gradient computa-

tion time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

 xiv

List of Acronyms and Abbreviations

AD automatic differentiation

API application programming interface

DAE differential algebraic equation

DDE delay differential equation

DE differential equation

DEQ deep equilibrium

DER distributed energy resource

DP differentiable programming

EMT electromagnetic transient

FFNN feed-forward neural network

IBR inverter-based resource

IVP initial value problem

LSTM long short-term memory

ML machine learning

NZE net-zero emissions

ODE ordinary differential equation

PINN physics-informed neural network

[Link] [Link]

PV photovoltaic

RMS root mean square

 xv

SciML scientific machine learning

TSA trajectory sensitivity analysis

 Chapter 1

Introduction

1.1 A Changing Power System

To limit the negative impacts of climate change, it is imperative to decarbonize the global

energy sector as quickly as possible. Corporations and governments throughout the world have put

forth common goals to reach net-zero carbon emissions by 2050 [1] as climate science suggests that

this target is roughly compatible with limiting warming to 1.5 degrees Celsius above pre-industrial

levels [2]. To accomplish such a rapid transition from a fossil-fuel based system to a renewable

system will require dramatic changes across sectors. In [3], the International Energy Agency lays

out a macro-level roadmap for transitioning all sectors to net-zero emissions (NZE) by 2050.

One promising approach to widespread decarbonization is to decarbonize the electric power

system while simultaneously electrifying other sectors. This is a popular strategy given that the

electricity sector has access to established and cost-effective technologies for zero-carbon genera-

tion. Figure 1.1 shows the final global energy consumption across industry, transport, and building

sectors in the years leading up to 2050 for the NZE scenario. The common trend in each sec-

tor is that electricity makes up a larger share of energy consumption as part of the pathway to

overall decarbonization. From the power system perspective, electrifying other industries means

incorporating additional load into the system. One example of this trend is the adoption of electric

vehicles. Figure 1.2 shows the increasing share of electric vehicles across different vehicle categories.

While smaller vehicles are more easily electrified, in the NZE scenarios all vehicle categories are

majority battery electric by 2050. Other industries such as aviation, shipping, and certain indus-
 2
Figure 1.1: Final global energy consumption across industry, transport, and building sectors
leading up to 2050 for the NZE scenario.1

Figure 1.2: Global electric vehicle share by vehicle category in the years leading up to 2050 for
the NZE scenario.1

trial processes are more difficult to electrify, hence the transition will be delayed; however, progress

is being made even on these more difficult sectors. In summary, for many industries, converting to

electricity is the most straightforward path to decarbonization.

1
Image: IEA 2021; NetZero by 2050 A Roadmap for the Global Energy Sector, [Link]
zero-by-2050, License: CC BY 4.0.
 3
Figure 1.3: Global electricity generation by source leading up to 2050 for the NZE scenario. Wind
and solar PV are expected to dominate the generation mix by 2050.2

While converting other industries to electricity amounts to the addition of new loads, de-

carbonizing the power system itself requires replacing fossil fuel-based generation with renewable

sources of generation. There are many different carbon-free sources of electricity generation—

including hydroelectric, nuclear, geothermal, wind, and solar photovoltaic (PV). Figure 1.3 shows

the breakdown of total electricity generation by source for the years leading up to 2050 in the

NZE scenario. While the future power system will undoubtedly include a diverse generation fleet,

wind and solar PV stand out as well-established and affordable options with the potential to scale.

Both technologies have experienced rapid decreases in cost to the point that they are the most

economical source of new generation in many areas of the world [4]. The global average levelized

cost of energy for onshore wind generation and solar PV decreased by 69% and 89% respectively

between 2010 and 2022 [5]. While there are significant locational differences in the quality of wind

and solar, the resource adequacy is not as restrictive as other technologies such as hydro-electric or

geothermal.

The adoption of new technologies for both load and generation leads to fundamental changes
2
Image: IEA 2021; NetZero by 2050 A Roadmap for the Global Energy Sector, [Link]
zero-by-2050, License: CC BY 4.0.
 4

in the system itself. A simplified summary of these changes are depicted in Figure 1.4. Three

general trends can be observed in this transition. First, power system generation is becoming more

distributed as large centralized generation facilities are replaced with many smaller distributed

energy resources (DERs). Second, the power system is increasingly inverter-based as many of the

fastest growing technologies are interfaced to the grid through power electronic converters. Finally,

with the addition of DERs and so called “smart-grid” technologies, the potential exists for the

coordination of many devices and more active control of the system generation and load.

Figure 1.4: An illustrative overview of how technology additions manifest in the future power
system. Image adapted from [6] and reproduced with permission. (©2020 NREL).

Present Future

Large centralized generation Distributed energy resources

Machine-based Inverter-based
Passive load Active load

1.2 Power System Stability

The fundamental role of the electric power system is to reliably deliver electricity from gen-

eration to load. The power system is designed to exist in a state of quasi-equilibrium wherein the

entire system remains synchronized and the physical states of the system remain close to their

nominal values. The power system is said to be stable if it can maintain this steady state during

normal operating conditions and can return to such a state when subjected to credible disturbances.

To capture the wide range of dynamic phenomena that are observed in power systems and impact
 5

stability, it is useful to categorize the concept of stability further, as in [7]. For example, rotor angle

stability refers to the ability of synchronous machines throughout the system to remain synchro-

nized, voltage stability refers to the ability of the system to maintain voltage levels close to nominal

values, and frequency stability refers to maintaining the system frequency close to nominal by bal-

ancing the generation and consumption of active power [7]. While the exact timescales associated

with specific phenomena vary, power system dynamic stability generally refers to dynamics at the

sub-second to minute timescales. Problems that can arise within any of these categories must be

addressed in order to maintain power system reliability.

In order to anticipate and mitigate any stability concerns before they appear in the real power

system, it is important to have reliable methods for studying dynamic stability in power systems.

Of course, the concept of dynamic stability is not unique to power system, but exists throughout

the physical sciences in any system that experiences change. In some cases, there exist analytical

solutions which describe the change in the system; unfortunately, this is rarely the case for all but

the simplest toy systems. Other tools for studying dynamics, such as small-signal analysis and

Lyapunov stability theory, are limited based on the size or perturbation that can be considered and

on the properties of the system. For studying power system dynamics that can be highly non-linear

under large disturbances, we typically rely on transient numerical simulations, wherein a dynamic

model is numerically integrated in time in order to solve for the trajectory of the system and make

conclusions about the stability.

1.3 Transient Numerical Simulations

In light of the ambitious goal of transitioning to a 100% decarbonized power system, the

tools and analyses that are used to study and operate the power system must also change to reflect

the realities of the new system. The adoption of new tools is already well underway in some areas.

For example, the variable nature of wind and solar generation compared to fossil fuel generation

has led to new frameworks for quantifying and valuing flexibility within the power system [8].

Furthermore, the uncertain nature of these resources has led to an increased reliance on weather
 6

forecasting for reliable power systems operations as the coupling between weather and generation

is increased [9]. In a similar manner, the changes in the power system described in the previous

sections create additional challenges for running transient numerical simulations. The overarching

goal of this thesis is to help address these challenges.

When running a transient numerical simulation, there are a plethora of modeling and simu-

lation choices depending on the application of interest. To make these choices, the power system

engineer must consider the inherent trade-off between the complexity of the model and the error

relative to the real system. In general, a more detailed model will require more computational

resources to solve, therefore computational burden can approximately represent model complexity.

In Figure 1.5A, this trade-off is represented simplistically as a straight line. At one extreme, we

might represent millions of individual devices with a single equivalent generator model given by

the second order swing equation; such a simplification is common practice for “external” portions

of a system that are far from the area of study. This model, though exceptionally easy to solve,

is inappropriate for any analysis that requires spatial sensitivity (e.g. how the voltage changes at

a specific bus). On the other extreme, we can develop system-level models that include switching

models of individual inverter-based resources (IBRs). Such a model will be difficult and slow to

solve, but we expect it to more faithfully reflect the real system. Other common simulation frame-

works, such as root mean square (RMS) simulation and waveform simulation with averaged models,

exist in the space between these two extremes. None of these modeling formulations exactly reflect

reality, but each can be useful in an appropriate setting.

For any application that makes use of transient numerical simulations, there is some application-

dependent maximum acceptable level of error and computational burden. Violating these con-

straints is depicted in Figures 1.5B-E as a shaded gray area. In other words, we require a tool

that achieves a trade-off of error and computational burden in the unshaded area. When analyzing

relatively small systems (Figure 1.5B), it is possible to simply add additional modeling detail until

the error is reduced without rendering the model intractable. Of course, for a larger system (Figure

1.5C), the relative computational burden is increased. In this case, it is no longer tractable to
 7
Figure 1.5: Representative trade-off between model error and computational burden for different
types and size of power systems.
 8

use the most detailed available models, but there are still tools that meet both of the constraints.

For example, RMS simulations with an algebraic network model and custom solution methods (i.e

Positive Sequence Simulations) have traditionally been used to model very large systems while

successfully capturing the dynamics of interest [10].

The addition of new technologies into the power system has fundamentally altered the trade-

off between model error and computational burden (Figure 1.5D) for two key reasons. First, given

that large power plants are generally being replaced with smaller DERs, the number of devices in

the system, and therefore the number of dynamic states, tends to increase. The second impact is

more nuanced and is related to the timescale of the controls for IBRs. Figure 1.6 shows the range

of timescales for various physical processes that exist within the power system. The timescale of

inverter-based controls tends to bridge the gap between slow electromechanical phenomena and

faster electromagnetic phenomena. The dominant numerical paradigms for simulating large sys-

tems are based on the separation of distinct group of timescales such that the faster timescales can

be neglected without degrading the overall solution. With the introduction of inverters, these sim-

plifying assumptions are no longer necessarily valid. The end result is that the range of timescales

that must be considered in order to have an accurate simulation is increased. Generally speaking,

problems that contain widely varying timescales (i.e. “stiff” problems) are notoriously difficult to

solve and require specialized numerical methods. In sum, these two effects mean that the current

changes in the power system tend to increase the computational burden of transient simulations.

For large power systems with high penetrations of IBRs, it is not always possible to achieve

reliably accurate results with reasonable computational complexity (Figure 1.5D depicts a simplistic

representation of this shortcoming). The impact of this shortcoming can be seen in both the

research community and the actual operation of the power system. Many researchers have noted the

discrepancy between different modeling paradigms when including IBRs due to the representation

(or lack thereof) of fast network dynamics [12, 11, 13, 14] and/or fast inverter-based resource

controls [10, 15]. Operationally, RMS type simulations are the standard for large scale simulations

because they can be run much faster than electromagnetic transient (EMT) simulations. In multiple
 9
Figure 1.6: The range of timescales relevant to power system dynamic simulations. Image repro-
duced from [11]. (©2023 IEEE).

How to model?

Lightning propagation Electromechanical

phenomena
Switching surges
Inverter-based controls
Thermodynamic
Stator transients and phenomena
subsynchronous resonance
Rotor angle dynamics
Governor and load
frequency control
Voltage control
Boiler dynamics

Wave Electromagnetic
phenomena phenomena

10−7 10−5 10−3 10−1 101 103 105

Time (seconds)

recent real grid events, however, the RMS models were unable to capture the true dynamics of the

system [16, 17]. For medium-sized grids with high penetrations of renewables such as those in Maui

[18], Texas [19] and Australia [20], operators have begun to build EMT models in order to capture

system behavior with higher fidelity. For even larger systems (e.g. the Eastern Interconnection),

however, the computational limitations are restrictive, and building and running EMT simulations

is simply not currently feasible.

The goal of this thesis is to improve the inherent trade-off between model error and computa-

tional burden for simulating large power systems through the addition of novel modeling and sim-

ulation techniques (Figure 1.5E). There are many different strategies for addressing this trade-off,

including novel modeling techniques [21, 22, 23], solution methods [24, 25], software implementa-
 10

tions [25, 26], parallelization frameworks [27, 28], and hardware [29, 30]. The reader is directed to

[31] and the references within for a more detailed review of techniques for accelerating EMT sim-

ulations. Ideally, many small improvements in any of these areas can be implemented in parallel,

resulting in significant overall improvements.

In this thesis, I propose a variety of data-driven machine learning (ML) models that are

used to represent power system components within a simulation along with tools for training such

models within a simulation framework. ML models are considered for this application due to their

ability to capture complex relationships and their speed of evaluation (after training is complete).

Importantly, this thesis does not propose replacing transient simulations with ML models, but rather

integrating ML models into a traditional power system simulation framework. This is an important

caveat as it maintains the use of a well-trusted tool within the industry and helps to address some

long-standing shortcomings of data-driven models, such as limited interpretability. A more in-

depth discussion of the methodology for including data-driven models is included in Section 2.4,

and issues related to the adoption of such models are discussed throughout the thesis.

1.4 Research Impact

New techniques that can improve the trade-off between computational burden and model

error will have a profound impact on power systems in both an operational and research setting. In

essence, the improvement would allow for either running simulations of a more accurate model with

the same computational burden or running simulations with similar accuracy and reduced compu-

tational burden. The following subsections describe selected examples of power system processes

that rely on transient simulations and will therefore benefit from this research.

1.4.1 Contingency Analysis

Contingency analysis refers to the process of checking a selection of critical contingencies in

real time to anticipate any serious stability concerns that might arise from such contingencies and

take preventative actions. A common framework for contingency analysis is N-1, wherein the loss of
 11

any single components (typically lines or generators) is considered. Due to the real time necessity

in contingency analysis, the impact of faster simulations is clear—given a finite set of time, if each

contingency can be analyzed faster, more contingencies can be considered and the resulting analysis

is more comprehensive.

1.4.2 Interconnection Studies

Interconnection studies refers to the studies that are conducted before connecting new gen-

eration resources to the grid to ensure there are not adverse impacts on the system. While the

time constraints are not as extreme as in contingency analysis, there is a need to capture enough

detail in a computationally tractable formulation. The median interconnection wait time was five

years for resources built in 2023 [32]—using improved dynamic simulations, which are computa-

tionally tractable and accurate, is one way to reduce the time spent in the interconnection queue

and alleviate this bottleneck for adding new renewable generation.

1.4.3 Validation of Novel Controls

The impact of this research on future power systems research is somewhat harder to quantify.

Nonetheless, transient simulations are used ubiquitously to validate novel contributions in the field

of power systems. Given the rapid and accelerating transition to IBRs, a major research focus

within the power system community is on developing new controls and methods for not only the

current power system, but also future power systems, which will inevitably have much higher

penetrations of renewables. Given the scale and complexity of real power systems, it is beneficial

to demonstrate novel ideas on large and detailed models as they more closely reflect reality. While

there are fast simulation methods that have worked historically, it is imperative to re-examine

the fundamental assumptions behind such methods, as they are not guaranteed to hold for new

inverter-based systems. In this sense, having the flexibility to employ more detailed models that are

still computationally tractable can be a valuable research tool. Additionally, given the proliferation

of “smart grid” devices, a strong research area is focused on the aggregation and control of large
 12

numbers of heterogeneous devices in order to achieve desirable system-level behavior. These works

require sufficiently large systems to validate the approach. In summary, the impact of this research is

to enable simulations, in both research and operational settings, of large power systems in sufficient

detail such that we can accurately capture important dynamics that can impact the stability of the

system.
 Chapter 2

Preliminaries and Literature Review

2.1 Preliminaries

The topic of this thesis lies at the intersection of three distinct areas: scientific computing,

ML, and power system dynamics (Figure 2.1). The first three subsections serve as a broad overview

of these three areas in order to establish the requisite baseline knowledge for the remaining chapters.

Readers that are familiar with any of these individual areas can skip the appropriate section without

missing the contribution of this thesis. The final subsection addresses the intersections between

these areas, including synergies and challenges for combining them.

2.1.1 Power System Dynamics

Power system dynamics is the study of how the power system evolves over time. A wide

range of timescales are relevant to the operation of the power system. For example, on the decadal

timescale, generation and transmission investments must be made to match the changing electric

load. On the hour to day timescale, specific generators must be dispatched to meet the daily fluc-

tuations in load. While these long term analyses address a changing system and are practically

important, they typically assume a quasi-static condition. Power system dynamics generally refers

to the changes in the system that occur on shorter (sub-second to minute) timescales, and is closely

linked to the concept of dynamic stability. Maintaining stability is a fundamental prerequisite for

accomplishing the overarching goal of the electric power system—to reliably and economically de-

liver electricity from generators to loads. For example, unforeseen instabilities in the system can
 14
Figure 2.1: The three major topics which comprise the content of this thesis are scientific com-
puting, ML, and power system dynamics. Scientific machine learning (SciML) is a rapidly growing
research area with applications in many fields of science and engineering. Operationally, scientific
computing is the dominant paradigm for addressing power system dynamics problems, but the
application of ML is common in a research setting.

Scientific Machine
SciML
Computing Learning

Power System
Dynamics

lead to remedial actions that interrupt service, costly damage to equipment, or worse, cascading

blackouts impacting large portions of the system [33]. There are a diverse collection of dynamic

phenomena that are observed in power systems; the authors in [34] systematically categorize differ-

ent types of stability problems. Each of these problems are important to study and understand as

all can have detrimental impacts on the ability to provide reliable and economical electricity. This

categorization reinforces the idea that stability problems can occur on a wide range of temporal

and spatial scales, from fast local oscillations to much slower inter-area oscillations which can occur

across thousands of miles [35].

The dynamic characteristics of the power system are determined by the network topology

and the devices that comprise the system, including generators, loads, and the network itself
 15

(e.g. transformers, lines, etc.). In response to the need to decarbonize the electric power system,

power systems throughout the world are experiencing a rapid and accelerating adoption of new

technologies. Many of these new technologies (e.g. solar PV, wind turbines, electric vehicle chargers,

etc.) are interfaced to the grid through power electronic converters. In contrast, many fossil-fuel

based technologies, which are in turn being phased out, are interfaced to the grid though electric

machines. A growing body of research suggests that this transition fundamentally changes the

dynamic behavior of the system [36, 12, 13, 14, 15], which in turn leads to new types of stability

problems that were not encountered prior to the introduction of these devices [7]. While the

dynamic stability of the system was historically dependent on the physical inertia provided by

large synchronous machines, the future system comprised of largely converter-interfaced devices

will have greatly reduced inertia. Quantities such as frequency, which are physically interpretable

in the context of synchronous machines, are harder to define in the future system. The successful

transition to a 100% renewable energy power system will depend on our ability to understand the

dynamic behavior of the system in order to anticipate and mitigate stability problems.

2.1.2 Scientific Computing

Figure 2.2 displays a general categorization of methods to study any physical system of

interest. At the highest level, one can experiment on the system itself to deduce its characteristics,

or experiment on a model of the system in order to understand the original system. Furthermore,

models can be split into physical models and mathematical models of the original system. Finally,

some mathematical models can be solved analytically, while others require numerical simulation to

practically find solutions. The field of scientific computing refers to the application of advanced

computing capabilities to understand and solve complex physical problems, typically by formulating

mathematical models to represent the real world and applying numerical analysis to solve such

models [37].

One type of analysis common to scientific computing across many application areas is the

formulation and solution of models based on differential equations (DEs). While DEs are by no
 16
Figure 2.2: A categorization of ways to study a physical system. scientific computing or numerical
simulation (highlighted) involves the simulation of a mathematical model of the system to study
and understand the true system.1

System

Experiment Experiment
with actual with a model of
system actual system

Physical Mathematical
model model

Analytical Numerical
solution simulation

means the only type of problem encountered in scientific computing, they serve as an easy-to-

understand example and are relevant to the types of problems addressed in this thesis. The first

key component is the mathematical model of the physical system, in this case a set of DEs for the

dynamic states of the system that describes how the states evolve in time:

dx
= f (x, t) (2.1)
dt

The second key component is the initial condition or starting state of the model:

x(0) = x0 (2.2)

The third and final component is a numerical method for iteratively solving for the time domain
1
©2010 Eduard Babulak and Ming Wang. Adapted from [38]; Originally published under CC BY-NC 4.0 license.
 17

trajectories of the system states in time:

ODE SOLVE : [x0 ... xt ] (2.3)

Scientific computing is a mature field, and there are therefore a wide array of solution methods for

solving DEs. The choice is not one-size-fits-all; rather, the best solution method depends on the

characteristics of the problem and application requirements (e.g. time to solve, required tolerance,

etc.).

All three of these components represent simplifications of reality and can therefore be a source

of error in our understanding of the real system. For example, it may be difficult to develop an

accurate model of the system due to a lack of understanding of the fundamental processes. In

some cases, the states of the mathematical model might not have physical meaning or be easily

measured, therefore it can be difficult to find valid initial conditions. Finally, all numerical solution

techniques are approximate, therefore care must be taken to control the error while solving to

obtain a sufficiently accurate solution.

2.1.3 Machine Learning

ML can be defined as the development of statistical algorithms that can learn from data and

generalize to unseen data without being explicitly programmed.2 ML is a rapidly evolving research

area that has made fundamental advances in many disparate fields. A comprehensive overview of

ML is well beyond the scope of this thesis; the interested reader is directed to [39] for a foundational

overview of the field. Instead, as in the prior section, we introduce the main components of ML via

a small illustrative example.

At its core, ML is learning from data, therefore a key component is a training dataset.

Consider a dataset consisting of both inputs or “features”, x, and outputs or “targets”, y:

Dtrain = {x, y} (2.4)

2
In this thesis, “machine learning” and “data-driven” are used as adjectives interchangeably, typically to distin-
guish from other techniques that do not rely on learning from data.
 18

This form of dataset with input and output pairs (i.e. labeled data) is the basis of supervised

learning, the most prevalent of three major categories of ML. The aim is to learn a model that

takes the inputs and predicts the output of interest:

ŷ = N (x) (2.5)

Importantly, the goal is not to simply learn a model N which performs well for the training dataset,

but rather a model that can generalize and perform well for unseen data. The process of learning

the model is an optimization problem:

min L(y, ŷ) (2.6)

The aim is to minimize a loss function L which encodes the quality of the prediction compared to the

ground truth data. Given the non-linear nature of many practical ML problems, the optimization

algorithm is often an iterative procedure.

While the type and structure of the model are often informed by the underlying processes that

are being modeled, the power of ML comes from approximating functions by learning from data

without explicitly specifying a model based on known first principles. The universal approximation

theorem for neural networks states that a neural network of sufficient size can approximate any

function [40]. Practically, ML models have outperformed other methods and humans on benchmark

tasks in image recognition and natural language processing [41]. The popularity and success of

such models is largely driven by advances that make the training of large models possible, such as

automatic differentiation (AD), libraries of popular models, and specialized hardware for training

models.

2.1.4 The Intersections

This section addresses a foundational motivating question for the work in this thesis: Why

combine the fields of power system dynamics, scientific computing, and ML? The power system is

considered critical infrastructure that is heavily relied upon for many parts of modern life. As such,

maintaining reliability of the system is of upmost importance, and direct experimentation on the
 19

system is almost never acceptable. This means we must rely on models to understand and operate

the true system effectively. Both scientific computing and ML are frameworks for modeling some

underlying system. The state of practice for operating large power systems is to apply a scientific

computing approach to understand the dynamic behavior of the system. While the application of

ML to predict power system dynamics is an active research area, such methods are not yet common

operationally.

Scientific machine learning (SciML), or the combination of scientific computing and ML, is a

thriving new research area with applications in many fields of science and engineering. The premise

of SciML is that the two areas can be viewed as complementary methods with a similar aim to

model and understand the real world. In this sense, SciML is a method for taking advantage of

the inherent strengths and mitigating the weaknesses of each approach. While this thesis focuses

on SciML for modeling and simulation, the scope of SciML is much broader and can include

tasks such as data-processing, experiment design, and decision automation [42]. The power system

is a high-impact application area that requires interpretability, robustness, and accuracy while

simultaneously dealing with the complexity and scale of the system. In a general sense, these

challenges make the power system an ideal candidate for SciML research [42].

2.2 State of Practice for Power System Transient Simulations

This section provides an overview of two widely used frameworks for power system simulation:

positive sequence simulation and waveform simulation (also commonly referred to as EMT simu-

lation). In reality, there are many potentially useful simulations that do not fit neatly into these

two frameworks. For a comprehensive taxonomy of power system simulation frameworks, including

the theory behind various simplifications, the reader is referred to [11]. Any practical solution for

simulating power systems relies on building system-level simulations by combining a collection of

constituent models according to the topology of the system. This strategy allows for model re-use

and for scaling simulations to systems of arbitrary complexity. Beyond this common approach,

there are many assumptions and design choices that go into building and solving a simulation.
 20

2.2.1 Waveform Simulation

Waveform simulation, also commonly referred to as “point-on-wave simulation” or “abc sim-

ulation”, closely models the real system with minimal simplifications. As the name suggests, the

states of the model are the sinusoidal, time-varying network quantities. For a three phase system,

each phase is modeled independently, therefore it is straightforward to include the effects of un-

balanced dynamics (e.g. single-phase faults) in the model. Waveform simulations are generally

chosen when a high fidelity model is needed. Studies that consider the impact of unbalance in the

system (e.g. single phase to ground faults), studies of very fast transient phenomena (e.g. lightning

strikes), studies that must consider the impact of switching events (e.g. protection studies), and

studies for selecting equipment specifications based on the range of possible operating conditions

(e.g. insulation coordination studies) typically rely on waveform simulation. Depending on the

specific application, common choices for modeling the network are the Bergeron model [43] or more

accurate frequency dependent models [44]. Similarly, at the device level, power electronic convert-

ers can be represented using detailed switching models, or simplified models which represent the

dynamics of the device averaged over the switching cycles. Waveform simulations are typically used

for detailed studies of specific equipment. However, as practitioners increasingly turn to waveform

simulations to accurately capture the dynamics of inverter-based systems, waveform simulations

can also be used to study system dynamics with generic device models.

The solution method for waveform simulation is informed by the time-varying nature of the

network quantities. The dominant integration technique is a fixed time-step method based on the

trapezoidal rule for integration and was originally proposed by Dommel [45]. The level of detail

included in the various constituent models should determine a reasonable simulation time step

and therefore the computational burden of the simulation. Given the lack of simplifications and

inclusion of fast dynamics, waveform simulations become very expensive, if even tractable, for large

systems ([18] reports approximately 600 s wall clock time for 1 s of simulation time for the medium-

sized Maui power system). Initialization of waveform simulations is another challenge for larger
 21

systems. In general, the system model begins in a de-energized state for waveform simulations and

must be carefully guided to the desired initial operating condition through an initialization process,

effectively increasing the needed simulation time span.

The power system community largely relies on commercial software (e.g. PSCAD, EMTP,

Simulink, PLECS, etc.) for waveform simulation. This can present a significant barrier to under-

standing the full simulation implementation as certain assumptions are obscured due to intellectual

property concerns. As waveform simulations have traditionally been employed for smaller systems,

many of the software options are built around a graphical user interface. To accommodate larger

systems, many tools now offer Python interfaces to automate simulation workflows.

2.2.2 Positive Sequence Simulation

Positive sequence simulation, also commonly referred to as “Balanced RMS” simulation or

simply “Transient Stability” simulations, refers collectively to multiple key assumptions. In Positive

Sequence Simulations, the sinusoidal network quantities are represented with quasi-static phasor

quantities, which represent the true signal by its magnitude and phase angle. This representation

means that the fast dynamics of the network itself are neglected as it is represented algebraically

within a system of differential algebraic equations (DAEs). In addition, Positive Sequence sim-

ulations assume a balanced system, allowing for representation of the full system with only the

positive sequence. Practically, these simplifications dramatically reduce the computational burden

compared to waveform simulations and make it possible to simulate very large power systems (e.g.

the entire Eastern Interconnection of North America). For this reason, positive sequence simu-

lations are often built using generic models as opposed to models that directly represent specific

equipment. Over the years, an extensive library of generic models for synchronous machines and

their controls has been developed and validated across multiple software offerings. In recent years

generic models for inverters have similarly been proposed and iteratively improved [46].

The most popular solution methods for positive sequence simulation are partitioned (explicit)

methods; the differential equation portion of the model and the algebraic network model are solved
 22

separately and iteratively. While this approach is computationally simple and can lead to very

fast simulations, choosing an appropriate time step and interfacing models with the network are

key considerations to maintain numerical stability. Unlike waveform simulation, the initialization

of positive sequence simulations follows directly from the power flow solution. In general, it is

more straightforward to achieve a specific operating point because the model is time-invariant and,

strictly speaking, should start in a steady state condition where all derivatives are set to zero.

Similar to waveform simulation, the most popular options for positive sequence simulation (e.g.

PSSE, PowerWorld, PSLF, etc.) are commercial and not open source. Positive sequence software

is commonly built around a tabular interface which is more easily applicable to large systems.

2.3 [Link]

The state-of-practice for power system simulation described in the prior section is charac-

terized by a reliance on two dominant simulation paradigms—waveform simulation and positive

sequence simulation. These paradigms have become immensely popular based on their historical

utility to model and predict power system dynamics with varying level of details. Unfortunately,

given the fundamental changes to the power system brought on by the increase in IBRs, these tools

are increasingly insufficient for studying the dynamics of large systems. Power system operators

and researchers have begun to note the inability of positive sequence simulations to capture rele-

vant power system dynamics, both operationally as in the case of large grid events that were not

predicted by the models in use [16, 17], and in a variety of research papers [12, 11, 13, 14, 10, 15].

When positive sequence simulations are insufficient, the common response is to turn to waveform

simulations in order to capture more detail. As previously discussed, waveform simulations can

quickly become computationally intractable for large systems. In addition, given the myriad differ-

ences between the two paradigms, simply noting that positive sequence and waveform simulations

of the same system have different responses is insufficient for a detailed understanding of how each

assumption and design choice impacts the final result. Despite these shortcomings, power sys-

tem practitioners continue to rely heavily on these frameworks given the convenience provided by
 23

trusted and mature software offerings, even as the dynamics of the power system in practice have

changed rapidly [16, 17, 47].

[Link] ([Link]) [26] is an open source power system simulation

toolbox, written in the Julia Programming Language, which is designed foremost to provide flex-

ibility in building and running simulations of large power systems. Julia is a high level scripting

language that is easy to use, with performance comparable to lower level languages like C, espe-

cially in the area of solving DEs [25]. While positive sequence and waveform simulation tools are

highly specialized and combine many assumptions about the modeling fidelity and solution meth-

ods, [Link] is designed as a research tool to fill the gap between these paradigms by exploring

novel simulation frameworks for the changing power system. This is especially important given that

we cannot be sure a priori what level of detail will be necessary to capture the important dynamics

in an all-inverter system. Multiple dispatch, or the dynamic dispatch of functions based on the

type of inputs, is a key characteristic of the Julia programming language. It allows for writing

libraries of easily extensible component models. Figure 2.3 shows the meta-model of a generic IBR

composed of various sub-components implemented in [Link]. Writing a new model is as easy

as overloading a few function definitions for the sub-component of interest, while all other model

code is reused. Finally, both Julia and [Link] are completely open-source; all assumptions are

available in the source code, which is an essential quality for understanding a research tool.

One key design criteria that makes [Link] unique compared to other software offerings is

a complete separation of the modeling and solution layers. This feature can be seen clearly by

considering the software dependencies for running a simulation shown in Figure 2.4. While the

data handling and model building are handled in [Link] and [Link] respectively, the

models are formulated according to the common interface established by [Link], which

allows direct access to the extensive library of solvers included in [Link] and

other well-established tools such as the Sundials suite of solvers [48]. This separation of modeling

and solution layers allows the user to take advantage of improvements in both modeling and so-

lution algorithms, thereby “stacking” speedups and improving the fundamental tradeoff between
 24
Figure 2.3: The [Link] inverter meta-model. New inverter models can be defined by modifying
a single sub-component.3

Electric Power Grid

Vri Iri

Filter
AC Voltage and Vriflt
Current Dynamics

Vricv cv
Iri
δolc
Reference Frame Conversion Pref
cv cv
Vdq , Idq

Converter Vdc Primary Energy Source

PWM Dynamics Energy-side Dynamics

ref
Vdq Vriflt , Iri
flt cv
, Iri

Inner Loop Control

Frequency Estimator
Voltage Control, Current Control
PLL Dynamics
and Virtual Impedance
ref ref
Vref δolc
ωolc Volc Iolc δolc
Pref
Outer Loop Control Qref
Active and Reactive Power
δpll ωpll
Control Dynamics Vriflt , Iri
flt

speed and accuracy (Figure 1.5). For example, any novel solution methods or improvements to ex-

isting methods in [Link] can be readily applied to power system models with

minimal effort. While many commercial power system software offerings rely on fixed-step solvers,

[Link] enables the use of variable-step methods as well, another potential path for acceleration.

The separation of modeling and solution methods implemented in [Link] is a key design criteria

that informs the work in this thesis.

 25
Figure 2.4: The software dependencies related to [Link]. The software design maintains a clear
separation between modeling packages and solution methods.3

[Link] [Link]
Data Structures Models

System Model
F (x, y), G(x, y)

[Link]
[Link]
Initialization

Jacobian Initial Condition

J(x, y) x0 , y 0
[Link]
Common Interface

Solver

2.4 Methodology Overview

As mentioned in previous sections, the state-of-practice for studying power system dynamics

is to follow the scientific computing approach and build and run numerical simulations. For the

purposes of this thesis, a transient simulation consists of three distinct steps (Figure 2.5A). The

input to such a transient simulation is a dynamic model of the system based on the topology and

individual components. The input includes the operating point of the system (e.g. power and

voltage set-points for generation buses and active and reactive power of loads). The first step is to

solve the nonlinear power flow equations for the given set-points. The power flow solution includes

the network quantities (voltage and currents) for every node in the system. Next, the dynamic

devices at each node are initialized in order to simultaneously satisfy the power flow solution and

the steady state criteria. Starting from the full set of initial conditions, numerical integration is
3
Image from [26]. License: CC BY 4.0.
 26
Figure 2.5: The typical framework for studying power system dynamics is to employ scientific
computing (blue). One alternative is to completely replace the transient simulation using ML
models (orange). In contrast, this thesis focuses on a hybrid methodology in which models are
trained and incorporated into the simulation framework as a means of retaining the interpretability
of a trusted tool while improving the tradeoff between speed and accuracy.

A) Transient Simulation

Physics Power Initialize Numerical Time-domain

Model Flow Devices Integration Trajectories

“Pure” Scientific Computing Approach (state-of-practice)

B)
Physics Time-domain
ML Model
Model Trajectories

OR

Physics Stability
ML Model
Model Classification

“Pure” Machine Learning Approach

Transient Simulation
C)
Physics Power Initialize Numerical Time-domain
Model Flow Devices Integration Trajectories

Model
Training
Transient Simulation
ML Model Power Initialize Numerical Time-domain
Physics Model Flow Devices Integration Trajectories

Scientific Machine Learning Approach (this thesis)

 27

used to solve for trajectories of the system states.

While the pure scientific computing approach is sufficient for small systems, this same ap-

proach cannot easily scale to very large systems. For stiff systems, the computational burden of

numerical integration for a system of n states can scale up to O(n3 ) [49]. Given that the numerical

integration step is computationally expensive, one approach is to simply replace the entire transient

simulation process with a data-driven model (Figure 2.5B). This approach can be useful for very

time-constrained applications (e.g. as a screening tool when evaluating contingencies in real time),

as the speed-up compared to the scientific computing approach can be very large. Many papers

have proposed ML models that predict the stability classification of a specific operating condition

and contingency directly (i.e. Transient Stability Assessment) [50, 51, 52, 53]. For example, in [51]

the proposed model for binary stability classification is based on gated recurrent units (GRUs) for

improved predicted performance and combined with a decision tree for interpretability. In [51], a

stacked auto-encoder model is used for dimension reduction, and the model consists of a voting

ensemble classifier based on a support vector machine model and a random forest model. In [52],

the authors additionally use a generative adversarial network to supplement the training data that

comes from numerical simulations. In addition, the model in [52] predicts a transient stability

index instead of a simple binary classification which offers slightly more detail (e.g. stable but close

to instability). A direct comparison between the methods above is not presented in the literature

largely due to the lack of established benchmarks in this area, which should be considered a critical

area of future work if data-driven models are to become prevalent operationally in power systems.

Despite these examples that employ data-driven models with some level of interpretability, by na-

ture of replacing the transient simulation these approaches inherently lack the the level of granular

interpretability needed for applications where the goal of running simulations is to understand the

cause of a particular instability and suggest actions to mitigate it.

The approach taken throughout this thesis can be seen as a combination of the ML and

scientific computing approaches (Figure 2.5C). ML models are trained and deployed within the

scientific computing workflow along with other physics-based models. By representing a portion of
 28

the power system with a data-driven model with fewer states or numerically advantageous qualities,

the overall speed of the simulation can be improved. In this work, the data for training the models

come exclusively from transient simulations, but it is possible that data-streams from the real power

system could be used to train or improve the data-driven models in the future. This thesis extends

many of the core design principles of [Link] to the task of incorporating data-driven models into

simulations. In chapters 3 and 4, this methodology is implemented for various types of ML models

by directly augmenting the capabilities provided by [Link]. In fact, a key contribution of this

work is formulating ML models so that they can be specified independent of the solution algorithms,

a key design feature of [Link]. Note that this approach is not intended to replace a “pure” ML

or scientific computing approach, but rather provide an additional option for applications where

both model accuracy and computational speed are constrained.

There are several key advantages of using the proposed hybrid approach compared to relying

exclusively on data-driven models. First, this approach allows the power system engineer to retain

interpretability in the areas where it is needed. Physics-based models are typically built from

first principles, therefore the states of the model are imbued with physical meaning (e.g. network

voltage and current, phase-locked loop frequency, etc.). In contrast, ML models are often “black-

box” models where the results are not easily interpretable. In many cases, the more complex and

higher-performing ML models (e.g. neural networks) are especially opaque. In many power systems

applications, transient simulations are performed not simply to obtain time domain trajectories,

but rather as an exercise to understand the system and often to take corrective action. In this

case, it is possible to retain the same interpretable models for an area of interest and represent the

surrounding areas with data-driven models. Second, by using data-driven models to represent a

part of the system instead of the entire system, the task of learning the dynamics becomes modular

and can in theory readily scale to very large systems in the same way that system-level simulations

are built from a library of simple components. In other words, this approach takes advantage of

a fixed network topology to break down and simplify the learning process. In contrast, trying

to learn a single data-driven model for the dynamics of a large system may suffer from the well-
 29

known challenges associated with increasing input data size (i.e. the “curse of dimensionality”).

Finally, this hybrid approach may offer an easier path to industry adoption as it is not a wholesale

replacement of a trusted workflow, but rather an improvement to the existing approach that does

not sacrifice the decades of experience, models, and other domain knowledge in this area.

2.5 Numerical Test Systems

In an ideal case, we could demonstrate the methods of this thesis on models of real power

systems which reflect the devices in use today or planning cases of real systems that reflect the

ongoing transition to renewables. Unfortunately, the power system is considered critical infras-

tructure and therefore access to such models is restricted due to security concerns. Instead, power

system researchers rely on test systems which are intended to reflect certain features of real systems

but do not correspond to physical infrastructure. Dynamic test systems are typically defined for a

single snapshot in time but must include dynamic models for the devices in the grid which adds

an additional level of detail and complexity compared to other test systems used for quasi-static

analyses.

There are a variety of existing test systems for transient analysis within the positive sequence

simulation framework; a collection of such cases are collected and hosted online in a variety of soft-

ware formats [54]. Many of these test systems are derived by augmenting existing IEEE test systems

which were proposed for quasi-static analysis with generic dynamic models and typical parameter-

ization [55]. In addition, recent works have proposed methodologies for automatically deriving

larger test systems. In [56], synthetic power systems are derived based on statistical sampling of

real confidential systems within the continental United States. For representing inverter-based re-

sources, these systems use generic standard library models which have been included in commercial

software packages in recent years [57]. In [58], the authors propose a similar but distinct method

for generating synthetic grids with limited input data, with the goal of easily developing many

synthetic topologies for ML applications. Inverters are represented using the technology agnostic

normal form proposed in [59]. In both cases, it is unclear if inverter-based resources are repre-
 30

sented in enough detail (i.e. including fast control loops) to sufficiently capture potential stability

issues, as is discussed in previous sections. For more detailed waveform simulations, there is lack

of standardization with respect to device models and large test systems. Operators are actively

developing and using waveform models in areas with high penetrations of renewables [18, 19, 20]

however these models are not made publicly available. Furthermore, for the purposes of this the-

sis, waveform models are assumed to be computationally intractable for the largest power systems

because they rely on due fixed-timestep methods to capture the fundamental grid frequency.

There is a lack of suitable existing test systems for our work which is meant to bridge the

gap between positive sequence and waveform simulations in terms of accuracy and computational

tractability for inverter-based systems. Developing test systems from scratch takes significant time

and engineering effort. This is especially true when implementing and validating new inverter device

models or debugging inverter-related stability issues. Therefore, there is a tradeoff between the time

and effort required to develop test systems and the relative benefit of demonstrating methods on

larger systems. In this work, we aim to demonstrate the methods on models that are large enough

to demonstrate the computational benefit of reducing the system size on the wall-clock simulation

time, without requiring unreasonable effort in developing the system itself and/or computational

resources for generating training datasets and training the models. In Chapter 3, we develop a

custom 36-bus system with over 400 dynamic states which includes detailed models of synchronous

machines, grid-following inverters, and grid-forming inverters. In Chapter 4, we use the popular

2-area test system system and modify the system to include inverters [60]. In addition, for a

computational comparison of the proposed methodology, we increase the dimension of the system

by splitting device models into multiple equivalent devices with scaled base powers. This serves to

increase the computational burden of solving the system without changing the time-domain result.

For benchmarking our methods in Chapter 5, we use the 9-bus system and derive larger systems

by connecting copies of the original system with limited modifications. These modifications are not

ideal as they do not reflect real power systems and they introduce unrealistic symmetries compared

to real systems. Despite this, these strategies provide a straightforward way to demonstrate scaling.
 31

Developing a more realistic large test case, possibly by augmenting the methods of [56] [58] with

improved inverter modeling, is an important area of future work and a prerequisite to developing

a full benchmark problem in this area (see Section 6.2.3).

2.6 Literature Review and Contributions

For the purpose of this thesis, we always assume there is a power system model that is

sufficiently accurate for a particular application. This system is considered the “ground truth” and

the goal of the learning problem is therefore to replicate as closely as possible the behavior of this

system. We assume a generic time-invariant power system model:

dx
= F (x, y, ψ) (2.7a)
dt
dy
= G(x, y, ψ) (2.7b)
dt

where x and F are the states and equations of the devices, y and G are the states and equations

of the network, and ψ is the parameters of the model. With this notation, the three steps of a

transient simulation are:

(1) Solve the network power flow equations 0 = G(x, y 0 , ψ) for the initial network states.

(2) Solve the device equations 0 = F (x0 , y 0 , ψ) for the initial device states.

(3) Advance the solution in time through numerical integration to a series of distinct time

points: {x0 , ..., xn }, {y 0 , ..., y n }.

As discussed in Chapter 1, the major shortcoming of the accurate ground truth model is the

computational tractability of running simulations. Therefore, we aim to define an additional system

model which captures the same behavior as the ground truth system with reduced computational

burden:

dx̂
= F̂ (x̂, ŷ, ψ̂) (2.8a)
dt
dŷ
= Ĝ(x̂, ŷ, ψ̂) (2.8b)
dt
 32

There are multiple ways for the model in (2.8) to be faster to solve than the model in (2.7).

One straightforward option is to reduce the size of the state space (i.e. dim(x̂) + dim(ŷ) <<

dim(x) + dim(y)) by aggregating models spatially. This processes is commonly referred to as

“reduced-order modeling’ [61]’, “surrogate modeling” [22], “model aggregation” [23], or “generic

modeling” [62, 63, 64] in the power systems literature. This general strategy—aggregating models

spatially for state reduction—is the primary use case for data-driven models presented in this

thesis. Another mechanism for speed up is by modifying or approximating the fast dynamics of

the system. For example, a positive sequence simulation can be seen as a simplification of an EMT

model wherein the fast dynamic states of the network have been converted to algebraic states.

While the overall dimension of the state space may not decrease, the algebraic network allows for

fast partitioned explicit solution methods. This process and other temporal reductions are based

on the separation of timescales and singular perturbation theory to make model simplifications. A

detailed overview of these assumptions can be found in [11]. Learning approximations of the fast

network dynamics is another potential application of data-driven models, however this mechanism

is not explored in detail in this thesis.

In this thesis, we propose data-driven models as a means to accelerate simulations while

retaining high accuracy. However, there are other approaches. In some cases, it is possible to take

advantage of the underlying structure of the true system to derive an exact reduced order model.

For example, in [61], the authors develop an exact reduced order model of n parallel connected

single-phase inverter models with the same structure. Of course, real models of distribution and

transmission grids are rarely so simple. In [65], the authors consider a more realistic model with

impedances between the inverters, but a single type of inverter control is still assumed. The

general strategy of this work and others is to aggregate many models of the same type into a

single aggregate model. This technique is the basis of widely used “generic” models such as the

composite load model [66] which consists of multiple devices (e.g. inverters, machines, loads) at a

single device and is meant to represent an entire distribution system comprised of multiple feeders.
 33
R
Figure 2.6: A useful categorization of ML models. N represents the data-driven model and
represents a numerical integration algorithm. A time-stepping model is used in conjunction with a
numerical integrator to advance the model state. This approach is similar to physics-based power
systems models that are widely used. In contrast, a direct-solution model predicts directly the
model state based on the previous state. In this thesis, Chapter 3 deals with time-stepping models
while Chapter 4 deals with direct solution models.

Time-stepping model

Direct-solution model

Note that all of these aggregation techniques require knowledge of the structure and parameters of

the underlying system, which is not always available in practice due to a lack of reliable information

or a lack of model transparency from manufacturers. In general, this type of aggregation technique

works well for simple cases, but might not be well-suited to a realistic power system with a wide

diversity of devices and controls. In such cases, data-driven models offer a potential alternative.

It is difficult to directly compare physics-based and data-driven approaches for model aggregation

without considering the application. The goal of this thesis is not to make a definitive claim in

favor of ML surrogate models over physics-based surrogate models, but rather to develop the basic

techniques, which are lacking in the literature to date, necessary to incorporate and evaluate data-

driven models. In the subsequent subsections, we focus exclusively on learning approaches, and

literature related to each of the individual chapters are discussed in more detail.
 34

The specific challenges related to integrating data-driven models in conjunction with existing

physics-based models (Figure 2.5C) depend on the quantity that is predicted by the data-driven

models. The authors in [67] propose a useful categorization of surrogate models into time-stepping

models and direct-solution models (Figure 2.6). A time-stepping data-driven model predicts the

derivatives of the states of the model, and is used in conjunction with a numerical integration

(i.e. time-stepping) algorithm in order to advance the states of the model to the next time point.

In contrast, a direct-solution model bypasses the need for the numerical integration and predicts

directly the states at a future time point based on the initial states. In a power systems context,

the inclusion of time-stepping data-driven models is more straightforward, as the prevalent physics-

based models are formulated as DEs and similarly solved by numerical integration. The inclusion

of direct-solution models is less straightforward and is the primary topic of Chapter 4.

2.6.1 Chapter 3: Relevant Literature

The surrogate proposed in Chapter 3 is based on a combination of two novel implicit ML

models: neural ordinary differential equation (ODE) [68] and deep equilibrium (DEQ) models

[69]. In general, an implicit deep learning layer is specified as satisfying an implicit condition (e.g.

0 = N N (x, y)) as opposed to an explicit relationship between input and output (e.g. y = N N (x)).

A neural ODE, one example of an implicit layer, can be seen as an extension of recurrent neural

networks to continuous time such that the derivatives of the states are simply the output of a neural

network:
dx
= NNΘ (x) (2.9)
dt
At first glance, a neural ODE seems like a perfect fit for power systems modeling; as an ODE it fits

seamlessly with other existing models during the numerical integration step, and the expressiveness

and flexibility of neural networks promises to capture the complex non-linear dynamics that are

possible in the power system. Indeed, since their introduction neural ODEs have been implemented

in various scientific computing toolboxes [70, 71]. In [72], a series of surrogates of varying complexity

are developed around a basic neural ODE building block. The authors found that adding encoder
 35

and decoder blocks around the neural ODE improves performance by either increasing or reducing

the dimensionality of the hidden space with respect to the inputs and outputs. The case studies

for the surrogate models include replacing a single exciter, a PV plant model, and a region of the

full network. In [73], neural ODEs are used to model the uncertainty of networked micro-grids.

The focus is on applying reachability analysis to a hybrid data-driven and physics-based model

to ensure the accuracy of the data-driven portions. Both of these works employ neural ODEs to

approximate the transient dynamics of a portion of a larger system; however, neither sufficiently

addresses the initialization of the models with respect to the larger system.

To address this limitation, we turn to the DEQ Layer proposed in [69]. The layer is inspired

by considering a weight-tied neural network NN with parameters Θ of infinite depth such that the

model converges to a fixed point value z ∗ :

z ∗ = NNΘ (z ∗ , x) (2.10)

The output of a DEQ layer by definition satisfies a fixed point equation. In this work, a DEQ layer

is used to encode the constraints during the initialization of the neural ODE, namely that the model

begin in steady state and match the power flow solution. After the DEQ-based initialization, the

time domain trajectories are solved using the neural ode model in conjunction with an ode solver.

In addition to addressing the initialization to meet the network requirements, the work in

Chapter 3 is also unique in considering the application of simulation acceleration. While other

works characterize the accuracy of the data-driven models with respect to the ground truth, we

characterize both accuracy and simulation time. We demonstrate simulations that are five-times

faster (for a relatively small power system) with accuracy comparable to a physics-based surrogate

model which relies on perfect knowledge of the underlying systems model structure and parameters.

2.6.2 Chapter 4: Relevant Literature

Chapter 4 deals with the integration of direct-solution type surrogate models into a power

system simulation. Recall that a direct-solution model predicts the state values at a future time
 36

without requiring a numerical integration algorithm. There have been numerous works which

incorporate these surrogates into a power system dynamic simulation. In [74], the first work in this

area, the authors use a black-box data-driven model to represent a distribution network while relying

only on terminal measurements. The authors use normalized deviations of voltage and current as

inputs and outputs respectively in order to learn a model that is general across operating points.

The results are promising; however, by using the deviation from steady state voltage as an input

the approach assumes that the dynamic behavior will not depend on the operating point, which

is unlikely to hold in some cases (e.g. an inverter operating with or without headroom). The

model uses past values of the voltage input along with a fixed-step solver with a time step equal

to the prediction interval of the data-driven model. Many works follow a similar strategy to [74]

with various advancements. [75] presents an almost identical methodology with a slightly different

model architecture applied to replacing a portion of the transmission system. In [76], the surrogate

model is used to represent the aggregated response of battery energy storage systems (BESS). The

surrogate uses two separate neural networks for real and reactive power responses, and a genetic

algorithm is used to determine the optimal topology for each. In [77], the surrogate uses the local

frequency measurement as an input in addition to voltage and current. The focus of this work

is on quantifying the uncertainty associated with active distribution networks, however this work

does not fully re-integrate the surrogate into a simulation, but rather relies on a “closed-loop”

configuration for testing the model.

While the works above indicate some promise in using direct solution surrogates within larger

simulations, recent advances have made the possibility more appealing. For example, physics-

informed neural networks (PINNs) are a type of data-driven model which uses a physics-based

model of the system in the loss formulation for training the model [78]. PINNs were first applied to

power systems in [79] as a proof of concept for predicting state values at future time points without

relying on any numerical integration. [80] provides a much more detailed comparison of PINNs

as an alternative to numerical simulations in a power systems context. Compared to “vanilla”

neural networks or other ML models, PINNs benefit from improved training efficiency, accuracy,
 37

and interpretability [80]. While PINNs are not the focus of the work in Chapter 4, they are an

example of a direct solution surrogate, therefore this work offers a path to integrating PINNs into

a simulation as opposed to only using PINNs as a standalone simulation replacement.

The primary contribution of this chapter is to include direct solution data-driven models into

a large simulation by formulating the system as a system of delay differential equations (DDEs).

In an ODE, the dynamics are based on the states at a single point in time:

dx
= f (x(t))) (2.11a)
dt
x0 = x(0) (2.11b)

In DDEs, the dynamic equations can also depend on past values of time. For simplicity, consider

a system of DDEs with a single constant delay of length τ :

dx
= f (x(t), x(t − τ )) (2.12a)
dt
x(t) = h(t), t < τ (2.12b)

In contrast to ODEs, which require an initial condition, a DDE requires an initial history function

(h) which specifies the values of the states for all times less than the largest delay [81]. The work

in Chapter 4 incorporates data-driven models through the formalism of DDEs. The benefit of this

approach is that it retains the separation of model and solution methods—an overarching goal

of this thesis. This makes the resulting model less restrictive; in a power systems context, this

distinction allows for the use of adaptive time-stepping solution methods which may be beneficial

for simulation speed and for solving to a tolerance instead of setting a fixed time-step a-priori.

In addition, adaptive solvers can accommodate the incorporation of multiple surrogates that have

different delay lengths in the same simulation.

2.6.3 Chapter 5: Literature Review

Unlike Chapters 3 and 4 which propose data-driven models and address how these mod-

els can be added to power system simulations, the impetus for the methods of Chapter 5 is to
 38

train such models directly within a power system simulation tool by applying the concepts of dif-

ferentiable programming (DP) to the typical power system simulation workflow. By making the

simulation end-to-end differentiable, the parameters of the data-driven models can be iteratively

updated according to the gradient and the chosen optimization algorithm. There are multiple

areas of literature relevant to this work. First, we briefly review the literature for differentiating

the constituent numerical problems of a power system simulation. Next, we review applications

of trajectory sensitivity analysis (TSA) within power systems, which is a closely related technique

for computing sensitivities of dynamic simulations. Finally, we review the literature related to the

software implementation of the methodology, namely implementations of AD in Julia.

For certain simple analytical functions, mathematically defining the derivative is easy (the

derivative of x2 is 2x). Computing the the gradient of the numerical solution of DEs is more nuanced

and challenging; there is significant nuance in defining the gradient and multiple options for the

computation. One way to categorize the existing approaches is by the order of the differentiation

and discretization steps: Continuous approaches calculate gradients by formulating a new system of

equations (differentiation) before solving the new system numerically (discretization). In contrast,

discrete approaches operate on the original solution (discretization) and propagate derivatives of the

unit operations of the solution (differentiation). In general these two approaches are mathematically

distinct and will not give the same exact numerical result (in other words, the discretization and

differentiation steps do not commute [82, 83]). The choice of sensitivity method can depend on the

solver; in [84], the authors show a simple case where the discrete method is non-convergent in the

case where the iterative solver has adaptive step size control that depends on the parameters that are

differentiated. Other works have demonstrated numerical oscillations in the gradient calculations

for discrete approaches to partial differential equations [85, 86]. In [87], the authors generalize that

continuous methods tend to be more efficient, while discrete methods tend to be more stable [87].

Approaches for computing the gradient of differential equation solves can also be categorized

as forward or reverse depending on the direction of propagation of the derivatives. In the DEs

community, reverse methods are commonly referred to as adjoint methods (e.g. [88, 89]). The key
 39

consideration for choosing a forward or reverse method is the relative input and output dimensions

of the function being differentiated. Generally, for both continuous and discrete methods, when

the output dimension is larger, forward methods tend to be computationally favorable, whereas

when the input dimension is larger reverse methods are preferable [84]. In addition to performance,

memory requirements can be another factor in choosing a method. Unlike forward methods, reverse

methods require access to intermediate values during the gradient calculation. To avoid saving the

entire solution, checkpointing techniques can offer a balance between storage requirements and

re-computation [90].

One example of a forward method is TSA, which has been widely used in power systems

applications since it’s introduction and efficient implementation in [91]. Again consider a simple

ODE with initial conditions:

dx
= f (x(t))) (2.13a)
dt
x0 = x(0) (2.13b)

We can define the solution to (2.13) as the trajectory of x:

x(t) = ϕ(x0 , t) (2.14)

The trajectory sensitivities are the time-varying partial derivatives that are derived by taking the

Taylor series expansion of (2.14) and neglecting all but the linear term:

∂x(t)
∆x(t) ≈ ∆x0 (2.15)
∂x0

These sensitivities quantify the variation of the trajectory in response to small changes in the initial

conditions. Sensitivities with respect to parameters are also available through the addition of trivial

parameter equations to (2.13)[91].

Since its introduction, TSA has been used extensively in the power systems community;

we include a non-exhaustive list of examples to demonstrate the range of applications. TSA has

been used for parameter estimation (e.g. for wind-turbine generators [92] and thermostatically
 40

controllable loads [93]) in addition to a Bayesian framework for estimation from noisy measurements

[94]. TSA has been used to determine which parameters are identifiable based on measurement

[95] and therefore reduce the number of parameters to be identified for composite load models

[96]. TSA has been used to predict the critical value of parameters such as the critical fault

clearing time [97, 98]. TSA has been used to bound the uncertainty of a nominal trajectory by

efficiently generating sets of trajectories based on perturbed parameter sets [99, 100]. TSA has been

used for optimizing control strategies [101, 102, 103, 104, 105], or for optimizing the placement of

components within the system to promote stability [106, 98]. Finally, TSA has been used as a

means of predicting instability and therefore allowing preventative control actions [107, 108, 109],

such as generation rescheduling [110].

The work of Chapter 5 can be viewed as an extension of TSA in a few key areas: 1) We include

an objective function with the state trajectories as inputs, 2) we provide flexibility in determining

how the gradients are computed (e.g. forward/reverse, continuous/discrete methods), and 3) we

explicitly consider the initialization of the models within the function that is differentiated. The

extension to include reverse methods for differentiation is especially key to ML workflows as TSA

is a forward method that will scale poorly with the number of parameters [88].

An important part of Chapter 5 is the software design and implementation for computing the

gradients of simulations with respect to parameters. One general approach is to start with a ML

library which is designed to be compatible with AD and implement the domain specific functionality

within this framework. There are many high-quality and open ML frameworks to choose from

(PyTorch, Tensorflow, and JAX are three popular options). Other libraries are built on these

frameworks to provide scientific computing functionality. For example, the NeuroMANCER library

is designed to combine scientific computing and ML functionality for various system identification

and optimization problems [70]. Another approach is to make generic code which is not explicitly

written in a ML framework compatible with automatic differentiation. The first approach is likely

more straightforward to achieve a working example. In [111], the authors implement a “neuro-

physical” load model for representing distribution systems using the NeuroMANCER library. The
 41

authors are limited to a single operating point, presumably because the functionality of a full

power system simulation tool is not available. In this work, we pursue the second option of making

existing generic code compatible with AD. In this way, we can retain all the functionality that

already exists within [Link], including the libraries of validated models. In addition, the potential

exists for including differentiable simulations within other capabilities implemented in Sienna such

as economic dispatch and unit commitment problems. In this way, the ability to differentiate

simulations becomes just an additional feature of these existing tools.

Multiple tools are employed to differentiate a full simulation. The core sensitivity methods

for the numerical sub-problems (solving nonlinear equations, solving DEs) are implemented in the

[Link] package [112]. The modularity of this approach mirrors that of the forward

pass used for [Link]. A general purpose AD library is used to to differentiate the higher level

functionality involved in running a simulation, such as the coordinated initialization of the devices,

parameter handling, and the user-defined objective function. A system for manually overriding the

derivative computation for certain functions allows for interfacing the general purpose AD with

[Link]. We considered two packages for this purpose: [Link] [113] and En-

[Link] [114]. Zygote is the more mature framework, but it does not have first-class support for

mutation, which is used extensively both in the Julia base library and within [Link] for handling

parameters and initial conditions. Therefore, we choose [Link], which has built-in support

for mutation and the potential for better performance due to differentiating compiler-optimized

code [114], despite being less mature. In addition, the [Link] application programming in-

terface (API) provides more nuanced handling of inputs that can be considered constant and do

not need to be differentiated, which is common in our application (we generally only care about

differentiating with respect to parameters). While the implementation in Chapter 5 is geared to-

ward ML/optimization applications due to a scalar objective function, using [Link] allows for

expanded flexibility to incorporate applications such as TSA as part of a more general purpose

differentiable simulation. In general, Chapter 5 serves to bridge the gap between ML and scientific

computing for power systems through DP.

 42

2.7 Summary of Publications

The following manuscripts comprise the foundation for this thesis (Chapter 3–5):

• M. Bossart, J.D. Lara, C. Roberts, R. Henriquez-Auba, D.S. Callaway, and B.-M. Hodge,

Acceleration of Power System Dynamic Simulations using a Deep Equilibrium Layer and

Neural ODE Surrogate, under review.

• M. Bossart, S. Chatzivasileiadis, and B.-M., Integrating Data-Driven Dynamic Power

System Models with Delay Differential Equations, under review.

• M. Bossart and B.-M., Differentiable Programming for Transient Power System Simula-

tions, under review.

The following publications are not included directly in this thesis:

• J.D. Lara, R. Henriquez-Auba, M. Bossart, D.S. Callaway, C. Barrows, PowerSimula-

[Link]—An Open Source Modeling Package for Modern Power Systems with

Inverter-Based Resources, preprint arXiv:2308.02921, 2023.

• M. Marković, M. Bossart, B.-M. Hodge, Machine learning for modern power distribution

systems: Progress and perspectives, in Journal of Renewable and Sustainable Energy, 2023.

• R.W. Kenyon, A. Sajadi, M. Bossart, A. Hoke, B.-M. Hodge, Interactive Power to Fre-

quency Dynamics Between Grid-Forming Inverters and Synchronous Generators in Power

Electronics-Dominated Power Systems, in IEEE Systems Journal, 2023.

• C. Roberts, J.D. Lara, R. Henriquez-Auba, M. Bossart, R. Anantharaman, C. Rackauckas,

B.-M. Hodge, D.S. Callaway, Continuous-time echo state networks for predicting power

system dynamics, in Electric Power Systems Research, 2022.

• A. Sajadi, R.W. Kenyon, M. Bossart, B.-M. Hodge, Dynamic interaction of grid-forming

and grid-following inverters with synchronous generators in hybrid power plants, in 2021

IEEE Kansas Power and Energy Conference Proceedings, 2021.

 43

• A. Ivemeyer, M. Bossart, R.W. Kenyon, A. Sajadi, B.-M. Hodge, D.K. Molzahn, Assessing

the Accuracy of Balanced Power System Models in the Presence of Voltage Unbalance, in

2021 IEEE Power and Energy Conference at Illinois Proceedings, 2021.

• R.W. Kenyon, M. Bossart, M. Marković, K. Doubleday, R. Matsuda-Dunn, S. Mitova,

S.A. Julien, E.T. Hale, and B.-M. Hodge, Stability and control of power systems with high

penetrations of inverter-based resources: An accessible review of current knowledge and

open questions, in Solar Energy, 2020.

• M. Bossart, R.W. Kenyon, D. Maksimović, B.-M. Hodge, The effect of power electronic

loads on western interconnection stability, in 2020 IEEE Power & Energy Society General

Meeting Proceedings, 2020.

 Chapter 3

Acceleration of Power System Dynamic Simulations using a Deep Equilibrium

Layer and Neural ODE Surrogate

This chapter presents the contents of the following manuscript:

Matthew Bossart, Jose Daniel Lara, Ciaran Roberts, Rodrigo Henriquez-Auba, Duncan Call-

away, and Bri-Mathias Hodge, “Acceleration of Power System Dynamic Simulations using a Deep

Equilibrium Layer and Neural ODE Surrogate,” under review.

3.1 Introduction

The transition from fossil fuel-based generation to renewable energy sources is accelerating.

Unlike traditional generation based on rotating synchronous machines, wind, solar PV, and storage

are IBRs that interface with the grid through power electronics. Synchronous machines and IBRs

have fundamentally different dynamic behavior, and therefore the adoption of IBRs has profound

impacts on the dynamic behavior of power systems.

Given the fundamental shift in the dynamics of the system brought on by IBRs, it is not

guaranteed that historically useful frameworks for the study of system stability will provide ade-

quate results for future systems. Recent work has shown that including line dynamics can change

system-wide stability classification [12, 11, 13, 14] and that it is difficult to accurately capture fast

inverter controls, such as the inner current loop, in the positive sequence domain [15]. Events in

Australia and Texas demonstrate the limitations of positive sequence tools to capture observed

system oscillations [16, 17].

 45

These recent challenges in IBR simulations highlight the need for more detailed modeling

(i.e. more dynamic states) which in many cases includes the representation of electromagnetic

transients. Furthermore, the combination of EMT dynamics and the fast control loops of IBRs

imposes a smaller time step requirement for solving systems with IBRs [115]. Considering that the

time to solve systems of DEs is largely determined by (1) the number of states in the system and

(2) the time step required for solving the system [116], simulating systems with large penetrations

of IBRs will be more computationally expensive. For a more in-depth discussion on the underlying

assumptions of various time domain simulation approaches, refer to [11]. The need for increased

detail in simulation models is at odds with the need to carry out studies that require simulating

multiple scenarios, such as contingency analysis or system impact studies [117].

A useful technique to accelerate simulations is to create simpler versions of complex models

called “reduced-order models” or “surrogates”. These surrogates can capture intricate system dy-

namics with fewer states, significantly reducing the computational demands during simulation. In

this manuscript, we use the term “physics-based surrogate model”1 2 to describe a surrogate with

a specified structure that is designed to replicate specific observed dynamics. Prior works have

developed “physics-based surrogate models” by representing many devices with a single aggregate

model of similar structure [61, 65, 66]. Developing these surrogate models requires knowledge of

the constituent models’ structure and parameters, which is not always available. For instance, sim-

ulations that involve manufacturer-provided black-box models will not have the internal structure

revealed to the model user. In addition, power systems with inverter control and DERs can demon-

strate complex dynamic behaviors that a structured model may not fully encapsulate in a general

sense. In this case, the physics-based surrogates must rely on distinct parameterizations in order

to capture the dynamic behavior accurately for different simulations [64, 118]. These drawbacks

can be mitigated by employing other types of surrogate models that are capable of emulating the
1
While we use “physics-based surrogate model” in this paper to distinguish from the data-driven model we
propose, this is not a term that is uniformly established or applied in the literature. Other terms such as “reduced
order models” or “aggregate models” are also used to describe this type of model.
2
“Physics-based surrogate model” should not be confused with “physics-informed neural networks”. The former
is a surrogate model based on the known structure of the true model. The latter is a neural network model that is
regularized during training using the DEs [79].
 46

underlying system’s behavior without relying on structural assumptions.

data-driven
model

Figure 3.1: For large systems, modeling all devices in detail (above) becomes computationally
intractable. This paper proposes a data-driven surrogate model that is designed to be integrated
into simulations containing physics-based models (below).

The “data-driven surrogate model” presented in this paper is capable of creating a reduced-

order representation with no underlying structural assumptions (Fig. 3.1). Unlike physics-based

surrogate models that require extensive design and validation, data-driven surrogates can be trained

on simulation data to capture the dynamics of a sub-section of a system. These surrogates are

adequate for modeling complex, non-linear behavior, and the number of states can be adjusted

to balance computational cost and accuracy, based on the complexity of the underlying system.

The framework also captures faster time-scale phenomena, such as network dynamics and inner

control loops of IBRs, when necessary. Both physics-based and data-driven surrogate models have

limitations and may not be appropriate depending on the level of detail needed. For example, as a

black-box model which aims to capture terminal behavior, the model presented in this paper does

not represent voltages and currents of internal buses and cannot be used to study faults within the

surrogate region. Data-driven surrogates have proven to be effective in accelerating simulations in

other fields, such as chemical reactions and heating systems [119].

In previous applications of SciML to power system dynamic simulation, the focus was on

completely replacing the transient simulation process. This involved either classifying stability (Fig.
 47
A)
System
Model Machine Learning Classify
Device Algorithm Stability
Setpoints

B)
System
Model Machine Learning Time Classify
Domain Stability
Device Algorithm Solution
Setpoints

C) Transient Simulation
System Physics + data-driven models
Model Time Classify
Power Initialize Numerical Domain Stability
Device Solution
Flow Devices Integration
Setpoints

Figure 3.2: Types of applications of ML to power system dynamics. A: Direct classification of

stability [50, 51, 53]. B: Predicting time domain trajectories [79]. C: Incorporating ML models
within the existing transient simulation workflow [72, 73, 77, 120]. This work fits into Fig. wherein
data-driven models are integrated into dynamic simulations.

3.2A) or predicting the time domain solution directly (Fig. 3.2B). We take a different approach by

integrating data-driven models into an existing transient simulation workflow (Fig. 3.2C). Prior

works that follow a similar approach have developed data-driven models based on variations of

recurrent neural networks to predict the current state based on a fixed sequence of prior states

[77, 120]. These surrogates are trained based on prior values recorded at a fixed simulation time

step; integrating them into a larger system therefore eliminates the possibility of using adaptive

time-stepping. A more natural choice for seamless integration is a neural ODE, which can be

seen as an extension of recurrent neural networks to continuous time [68]. The authors in [73]

use neural ODEs to model the uncertainty of DERs in systems of networked microgrids; in [72],

neural ODEs are used for parameter estimation and black-box modeling. However, neither of these

works addresses essential considerations related to integrating data-driven models (i.e., End-to-
 48
Section II. A Section II. B-D Section III. A-C

Solve Initialize Numerical Time

Physics
Power Dynamic Integration Domain
Models
Flow Devices Data

Model Training Evaluation

Data-driven
Surrogate Solve Initialize Numerical Time
Power Dynamic Integration Domain
Physics Flow Devices Data
Models

Figure 3.3: An overview of the proposed methodology. The different colors correspond to the focus
of each subsection of the methodology.

End implementation), such as initializing the models to start in steady-state and respecting the

boundary conditions imposed by physics-based devices and the power flow solution.

The novel surrogate structure proposed in this paper is based on DEQ layers [69] and neural

ODEs [68], which are both examples of implicit layers [121]. Implicit deep learning is a novel field in

which layers are formulated as satisfying implicit conditions (e.g. NN(x, y) = 0), rather than as an

explicit relationship between input and output (e.g. y = NN(x)). The key insight in this work is that

these novel layers can be used to encode the steady-state and network boundary conditions. These

layers resolve the initialization challenges of data-driven surrogate models allowing for seamless

integration into a system model.

The complete End-to-End (E2E) integration of ML surrogates is a seldom-discussed step

in other works and is one of the main focuses of this work. The contributions target the entire

data-driven surrogate pipeline, including development, training, initialization, and integration into

existing transient simulation workflows as shown in Fig. 3.2C. Further, we aim to build tools that

enable system operators to evaluate dynamic phenomena in grids with significant participation of

IBRs that employ heterogeneous control techniques and structures.

 49

The E2E methods developed in this manuscript are shown at a high level in Fig. 3.3 and

consist of the following steps: (1) A simulation is built and solved using physics-based models and

numerical integration to generate time domain trajectories. (2) The time domain simulation data

is used as the “ground truth” data to train a data-driven surrogate model. (3) The data-driven

model is used along with physics-based models to build and solve a new simulation. (4) The two

datasets are compared on the basis of accuracy and simulation speed.

The main contributions of this paper are:

• A data-driven surrogate structure, consisting of both DEQ and neural ODE layers (Section

3.2), which can be integrated into time domain simulation procedures.

• An E2E training, validation, and integration procedure for data-driven surrogate models

in time domain simulations including EMT dynamics.

• Numerical results for a medium-size case study demonstrating the implementation of the

proposed approach in an existing simulation platform (Section 3.4).

3.2 Surrogate Structure

The block diagram of the surrogate model is shown in Fig. 3.4. The proposed surrogate design

is informed by the requirements of the procedure followed by standard time domain simulation tools

described in Section 3.2.1. This work focuses exclusively on single-port surrogates (i.e. surrogates

that are connected to the rest of the system through a single point of connection). Extensions to

multi-port surrogates are left to future work.

3.2.1 Power Systems Time Domain Simulation Procedure

A time domain simulation requires the specification of two components: (1) a system model,

including the DAEs that describe the system’s physics and controls, and (2) a time-stepping (or

integration) algorithm [11]. In this work, we assume the system model consists of a set of time-
 50

invariant DEs:

dxsys
= fsys (xsys , y sys , ssys , η sys ) (3.1a)
dt
dy sys
= gsys (xsys , y sys , ψ sys ) (3.1b)
dt

where xsys , fsys (·), ssys , and η sys represent, respectively, the states, equations, references, and

parameters of the devices (e.g., inverters, machines, loads, surrogate models, etc.). The circuit

dynamics of the network are represented as the sub-system y sys and gsys (·) with network parameters

ψ sys .

Given the system model (3.1a)–(3.1b), we define the procedure to execute a time domain

simulation as the following three steps:

(1) Solve the network model power flow gsys (·), resulting in an initial condition for the network

states y 0sys .

(2) Based on the power-flow solution (y 0sys ), solve the steady-state equations for each device

to determine x0sys and obtain the system-level steady-state.

(3) Apply a disturbance and numerically integrate, starting from the initial conditions and

using a stepping algorithm.

The way in which each device interfaces with the network, and therefore the rest of the

system, is also an important feature of the simulation definition. We assume a current balance

framework wherein the input to each device is the voltage at the point of common coupling and

the output is an injected current [122]. With this formulation, a dynamic surrogate model can be

written as:

ẋ = f (x, s, v, η) (3.2a)

i = h(x, v, η) (3.2b)

where x are the dynamic states, v is the voltage input to the model, i is the model output current,

η are parameters of the model, and s are references which are determined during initialization.
 51
Implicit Reference
Transient Simulation: i′0 i′′ DEQ Layers Frame
0 x0
1. Solve Powerflow → v 0 , i0 i0 Nout x̂0 , ŝ
NNθ1 0 = NNθ2 (x0 , s, v 0 ) Explicit
2. Initialize Devices → x0 , s i0 = h(x0 ) Normalization
3. Numerical Integration Nin s Layer
v0
v ′0 v ′′
0
M(θ)

dx
Rest of v(t) → v(t) v′ v ′′ x i′′ i′
ODE [x1 , x2 ]⊺
dt
System
Nin NNθ2 N−1
out M−1 (θ) i(t)
← i(t) SOLVE

Neural ODE

Figure 3.4: The proposed surrogate model. The blocks are color-coded according to their function—
a key novelty lies in the combination of DEQ and neural ODE implicit layers with shared parame-
ters.

Function f represents the model’s DEs, and h is the mapping from the dynamic states to the

output current.

Definition (Initial Network Interface Requirements). Given v 0 , i0 as the initial voltage and current

from the power-flow solution, we define that a surrogate model satisfies the network interface

requirements if the following conditions are met:

0 = f (x0 , s, v 0 , η) (3.3a)

i0 = h(x0 , v 0 , η) (3.3b)

This definition implies that for a surrogate model to satisfy the initial network interface require-

ments, the differential terms must be initialized to a steady state while also matching the initial

voltage and current values of the power flow solution.

This requirement is critical because if the surrogate is not in steady-state, it is not possible to

find a full system condition that is in steady-state, which is typically the desired condition before

a perturbation is applied. Furthermore, if the surrogate does not match the power flow quantities,

the devices outside of the surrogate model have to compensate to achieve system balance leading

to a different initial state than intended.

 52

3.2.2 Network Interface

To reintegrate a surrogate into a simulation, the input and output network quantities (v and i)

need to be in the system reference frame. Internally, these quantities are first transformed to a local

reference frame (v ′ and i′ ) and then normalized (v ′′ and i′′ ) in order to enhance the performance

of the data-driven models. Without a local reference frame transformation, the solution of the

surrogate model would change depending on the arbitrary definition of the system reference voltage.

We define a local reference frame rotation to align one component of the input voltage to zero.
 
vi
v ′ = M (θ)v, i′ = M (θ)i, θ = arctan (3.4)
vr
 
 sin(θ) −cos(θ)
where M (θ) =  , v = [vr , vi ]⊺ , i = [ir , ii ]⊺ , and θ is the time-varying angle of the
cos(θ) sin(θ)
input voltage in the system reference frame. This technique is analogous to the reference frame

transformations which are common in many inverter and machine models [123, 11].

ML workflows include some data preprocessing steps so that inputs to neural networks are

normalized (e.g. bounded to [−1, 1]). Min-max normalization is included to normalize both the

inputs and outputs based on the range of values seen in the training data:
2(v ′ − v ′min )
v ′′ = Nin (v ′ ) = −1 +
v ′max − v ′min
(3.5)
′′ ′ 2(i′ − i′min )
i = Nout (i ) = −1 + ′
imax − i′min
where v ′min and v ′max are vectors of the minimum and maximum values respectively of each element

of v ′ across the entire training dataset. The inverse of the normalization is included at the output

of the model. Care should be taken to ensure the ranges of values in the validation and test datasets

are similar to those expected to be found in both the training dataset and in practical usage.

3.2.3 Initialization

The initialization portion of the proposed surrogate is a key differentiating feature with

respect to prior work. As mentioned in Section 3.2.1, the network interface requirements (3.3a)-

(3.3b) are explicitly considered during initialization. In [72], the initial conditions for the neural
 53

ODE are calculated by passing the power flow quantities to an explicit feed-forward neural network.

However, there is no guarantee that the resulting initial conditions satisfy the network interface

requirements. In contrast, the initialization of the proposed surrogate takes place through an

implicit layer—a DEQ layer to be precise—such that the forward pass of the model ensures the

relevant conditions are met. This feature implies that even the untrained surrogate matches the

network interface requirements (see Fig. 3.5). Unless the initial conditions of a model can be

calculated precisely in closed form (which is sometimes possible with physics-based models), the

initialization process is fundamentally implicit. Therefore it is reasonable to model using an implicit

ML layer. This approach embeds the network interface requirements into the structure of the

surrogate, and satisfying the requirements is not reliant on the training process.

The initialization portion of the surrogate is defined as:

(i0 , v 0 , p) → [x⊺0 , s⊺ ]⊺ (3.6a)

p = [p⊺1 , p⊺2 ]⊺ (3.6b)

s.t.

0 = NNp2 (x0 , s, v ′′0 ) (3.6c)

i′′0 = [x01 , x02 ]⊺ (3.6d)

where the system of equations are solved numerically with the initial condition:

[x̂⊺0 , ŝ⊺ ]⊺ = NNp1 (i′′0 , v ′′0 ) (3.7)

The initialization consists of an explicit prediction from the first neural network (3.7) followed

by a “correction” by the implicit layer (3.6c)-(3.6d). The equations that are solved in the implicit

layer guarantee that the surrogate satisfies the network interface requirements. The explicit pre-

diction layer is given by NNp1 —a densely connected feed-forward neural network with parameters

(weights and biases) p1 and the initial operating point as the input:

[i′′⊺ ′′⊺ ⊺ 1
0 , v 0 ] =z0 (3.8a)
 54

z1k+1 = σ(Wk+1
1
z1k + b1k+1 ) ∀k = 0, ..., K − 1 (3.8b)

1
[x̂0 , ŝ] = σ(WK+1 z1K + b1K+1 ) (3.8c)

The number of hidden layers K and the width of the hidden layers dim(bk ) ∀k = 1, ..., K − 1 are

chosen as hyperparameters of the training problem. The output dimension dim ([x̂0 , ŝ]) depends

on the number of states, which is also a hyperparameter. In summary, the initialization mapping

of the surrogate is given by an implicit-explicit layer similar to the layer proposed in [124] as an

isolated deep learning model for non-power systems applications.

The explicit layer is important in the overall structure of the surrogate as it allows the

surrogate to generalize over different operating points (contingent on good training). The explicit

layer is trained to give a good initial guess so that [x̂0 , ŝ] ≈ [x0 , s]. This implies that a different

operating point i0 , v 0 will yield different initial conditions x0 for the neural ODE and therefore a

different dynamic response. In comparison, prior works have dealt with multiple operating points

by defining the data-driven states relative to the steady state [77], implying that the dynamic

response of the model will be the same regardless of the operating point. This assumption is not

always true. For instance, an inverter operating without any headroom will behave much differently

from an inverter with sufficient headroom when faced with an increase in load.

3.2.4 Differential Equations

The dynamic portion of the proposed surrogate (f ) is given by a neural ODE [68] wherein

the derivatives of the hidden states are defined as the output of a neural network:

f (x, s, v, p) : = NNp2 (x, s, v) (3.9)

The inputs to the neural ODE include the dynamic states (x), the references which are determined

during initialization (s), and the time-varying exogenous input (v) which comes from the interface

with the rest of the system. The neural ODE neural network (NNp2 ) is also a full-connected

architecture with parameters p2 :

[x⊺ , s⊺ ,v ⊺ ]⊺ = z20 (3.10a)

 55
neural ODE states ground truth output neural ODE states ground truth output
predicted initial conditions surrogate predicted output predicted initial conditions surrogate predicted output
Neural ODE States

Neural ODE States

2 1

0 0

−2 −1

0 2 4 6 8 10 0 2 4 6 8 10

−0.85
−0.86
Real current (p.u.)

Real current (p.u.)

−0.88
−0.9 −0.9
−0.92

−0.95 −0.94

0 2 4 6 8 10 0 2 4 6 8 10

0.02 0.02
Imag. current (p.u.)

Imag. current (p.u.)

0
0

−0.02
−0.02
−0.04
−0.04
0 2 4 6 8 10 0 2 4 6 8 10
Time (s) Time (s)

Figure 3.5: Surrogate before and after training. The left column shows the behavior of the surrogate
before training, while the right column shows the same quantities after training. The top row shows
the trajectories for the hidden states (x) along with the predicted initial conditions (x̂0 ) shown as
green circles. The bottom two rows show the model outputs (real and imaginary current).

2
z2k+1 = σ(Wk+1 z2k + b2k+1 ) ∀k = 0, ..., K − 1 (3.10b)

2
ẋ = WK+1 z2K (3.10c)

The key difference is that the output layer is not bounded by an activation function and does not

include a bias. From a numerical perspective, the neural ODE describes a non-linear relationship

which defines the derivatives of the states and therefore the neural ODE architecture is compatible

with any numerical integration algorithm.

The final portion of the model is the mapping from the dynamic states (x) to the outputs

which interface with the rest of the system. In this case, the simplest possible option is chosen—the
 56

first two dynamic states are selected to represent the output and the inverse of the normalization

and reference frame transformations are applied to obtain the output current in the system reference

frame and units:

−1
h(x, p) := M −1 (θ)Nout ([x1 , x2 ]⊺ ) = i (3.11)

The complete surrogate model (Fig. 3.4) is given by (3.4) – (3.11).

One way to visualize the novel properties of the proposed surrogate is to consider the behavior

of an untrained surrogate model. While this is not common practice in a traditional ML model, in

this case the untrained behavior helps demonstrate the parts of the model performance which rely

on training and the parts which are inherent due to the structure. A sample result is shown for

the proposed surrogate both before and after training (Fig. 3.5). The left side shows the surrogate

before training, and the right side shows the same quantities after training. On the top row are

the dynamic states of the neural ODE and on the bottom two rows are the outputs of interest

(real and imaginary current). The predicted initial states (x̂0 ) are plotted as circles at time zero.

As expected, before any training, there is a significant difference between the predicted and actual

initial states x̂0 and x0 , and the dynamic response does not match the ground truth. However,

the hidden states are in steady state (flat line before the perturbation is applied) and the initial

outputs match exactly the ground truth data, precisely because both these constraints are encoded

into the DEQ layer. After the surrogate is trained, the prediction of the initial conditions is much

closer, as is the dynamic response.

3.3 Training Methodology

While the parameters of a physics-based model can often be derived from the physical system

itself, the parameters of the proposed data-driven models must be learned from data during the

training process. In practice, this means that the surrogate structure described must be fully

differentiable so that the gradient of each parameter with respect to the loss can be calculated

during training. The details of the training process are discussed in the following section.
 57

3.3.1 Data Generation

The dataset used for training and evaluation of the model, D, consists of n pairs of input

and output time-series trajectories where each trajectory contains m time points:

D = {{(v (i,j) , i(i,j) )}m n

i=1 }j=1 (3.12)

For the purpose of this paper, the input is the bus voltage at the point of common coupling

(v = [vr , vi ]⊺ ∈ R2 ) and the output is the injected current from the surrogate (i = [ir , ii ]⊺ ∈ R2 ).

The full dataset of n trajectories is randomly split into train, validation, and test datasets.

In this work all training data is generated from numerical simulations of the full system model

including EMT dynamics as employed in [13] (see Section 3.4.1 for implementation details). The

dataset consists of the network quantities (voltages, currents) at the point of common coupling

between the surrogate portion and the remainder of the system. Generating data by simulating

perturbations of interest in the full system is an easy way to ensure that the relevant dynamics

that exist in the test dataset are captured during training, even though the train and test datasets

consist of distinct trajectories that are generated independently. While the black-box structure of

the surrogate lends itself to applications in which an accurate physical model is not known and

only data is available (e.g. a real distribution feeder with measurement data at the sub-station

transformer), the extension to measurement data is beyond the scope of this work.

3.3.2 Loss Function

The purpose of the loss function is to quantify the performance of the surrogate model

during training so that the parameters can be iteratively updated (trained), thereby improving the

performance. The loss function is given by:

L(p1 , p2 ) = (1 − α)Linit. + αLdyn.

Linit. = RMSE ([x⊺0 , s⊺ ], [x̂0 ⊺ , ŝ⊺ ]) (3.13)

⊺ ⊺
 
Ldyn. = RMSE [i⊺(1) ... i⊺(m) ], [î(1) ... î(m) ]
 58

where RMSE is the well known root mean square error:

s
Pn 2
i=1 (xi − x̂i )
RMSE(x, x̂) =
n

p1 and p2 are the parameters of the explicit initialization prediction layer and neural ODE,

respectively. The initialization loss, Linit. , quantifies how well the explicit initial condition predic-

tion layer matches the actual initial conditions. The dynamic loss, Ldyn. , quantifies how well the

predicted output trajectories match the ground truth data. Both Linit. and Ldyn. depend on both

sets of parameters, therefore all parameters are trained simultaneously according to the total loss.

α is a tradeoff factor between the initialization and dynamic loss and is tuned as a hyperparameter

of the training.

3.3.3 Training Algorithm

The key challenge in training the proposed surrogate is to maintain the convergence of both

the DEQ and neural ODE layers throughout training. To the best of the authors’ knowledge, this

work is unique in combining a DEQ layer and neural ODE with shared parameters, hence there

are unique training challenges. In [125], the authors note the importance of equation scaling to

promote stability while training neural ODEs. For DEQs, it is desirable for the neural network to

be constrained in order to promote stability [69]. In this case, however, the output layer of the

DEQ neural network must be unconstrained because it represents the derivative of the model states

which should be allowed to take on all real values (3.10c).

In order to address the convergence of the DEQ layer, the training algorithm (Algorithm

1) includes a pre-processing step to ensure that the starting random initialization of the neural

networks does not lead to divergent behavior when the training begins. To achieve this, the neural

networks are randomly initialized and the convergence of the DEQ layer is checked for each entry

in the training set. This ensures that for the first training iteration the surrogate will converge and

we can calculate a meaningful value for both Linit. and Ldyn. . For the actual training portion, the

loss is calculated stochastically for a single trajectory for each training iteration and the parameters
 59

are updated according to the first-order Adam algorithm [126]. In order to prevent over fitting,

the loss on a validation set is checked at regular intervals during training. The training ends when

a maximum number of iterations is met. The hyper-parameters of the training process (Table

3.2) are tuned and selected based on the best performance on the validation dataset, and the final

performance results are evaluated on an unseen test dataset.

Algorithm 1 Surrogate Training Procedure

Require: Dtrain = (u(i,j) , y(i,j) ), kmax
while DEQ does not converge for all of Dtrain do
Randomly initialize p1 , p2
end while
while k < kmax do
Choose a perturbation at random to calculate loss.
update p1 , p2 according to Adam.
end while

[Link] [Link]
(power system simulator) (machine learning library)

1. Data Generation 2. Training

Ground Train / Surrogate
Truth Simulate
Validation Model
Power
System Simulate Test Update
Model
params

Ground Surrogate
Truth Packaging
Compare Surrogate
Response Power
Simulate
from System
Simulation Model
4. Evaluation 3. Integration

Figure 3.6: Overview of the computational setup for generating results.

 60

î

data-driven
surrogate
gfm
î gfl

physics
rest of system surrogate full model
Figure 3.7: The case study system. The right side shows the full-order model (unshaded) and the
two surrogate models (shaded) which can each be connected to the remainder of the system. The
performance criteria is based on how well the predicted current (î = [îr , îi ]⊺ ) matches the ground-
truth result (i = [ir , ii ]⊺ ).

3.4 Results

3.4.1 Computational Set Up

All of the code for generating data, training the surrogates, re-integrating the surrogates into

a simulation, and evaluating the performance is open source and written in the Julia programming

language. Julia is chosen for this work for its unique blend of high performance and ease of use.

The steps involved in generating the results that follow are summarized in Fig. 3.6. [Link]—a

package for simulating the transient response of power systems with large penetrations of IBRs—

is used to generate the train, validation, and test datasets [26]. The surrogate model shown in

Fig. 3.4 is built using the packages of the Julia SciML ecosystem [127] along with [Link] (a

ML library) so that the entire process of initializing and numerically integrating to solve for time

domain trajectories is fully differentiable. The model is trained via back-propagation through the
 61
Table 3.1: Summary of dynamic models.

Component Model Dynamic States

Transmission lines π-model 6
Machines GENROE 4
Excitation System Type II AVR 4
Turbine Governor Type II 1
Inverter A Grid-following [128] 16
Inverter B Grid-forming - Droop [129] 15

numerical solution of the DEs of the surrogate model. After training, the composable nature of

Julia along with the design of the surrogate make it possible to re-integrate trained data-driven

models into [Link] by defining a new device model. Finally, the resulting system, which includes

both physics-based models and the data-driven surrogate is simulated to generate a dataset of

predictions for comparison to the ground truth data. For the full details of the packages used and

the experiments presented in this section, the reader is directed to the paper repository3 .

3.4.2 Case Study

In this section, the training and evaluation of the proposed surrogate are demonstrated on

the synthetic 36-bus power system shown in Fig. 3.7. The system is modeled in detail such that

the inverter inner control loops and line dynamics are not neglected. The dynamic models are

summarized in Table 3.1 and the full detailed descriptions can be found in the [Link] documen-

tation [26]. The 36-bus system as constructed contains 719 dynamic states in total, making it a

medium-sized system with non-trivial computational burden to solve.

[Link] Data Generation

The system is simulated across a range of operating conditions and perturbations to generate

the ground truth trajectories that make up the train, validation, and test datasets. For each

trajectory, the operating point of the base system is changed by randomly scaling the set points

between a fixed range for each generator and load that is part of the external system (i.e. the
3
[Link]
 62
1.02

Real voltage (p.u.)

1

0.98

0.96
0 2 4 6 8 10
0.1
Imag. voltage (p.u.)
0

−0.1

−0.2

−0.3
0 2 4 6 8 10
Real current (p.u.)

−0.8

−1

0 2 4 6 8 10
Imag current (p.u.)

−0.5

0 2 4 6 8 10

Time (s)

Figure 3.8: The trajectories of the train dataset. The top row shows the model inputs (voltage) in
gray, and the bottom row shows the model outputs (current) in black. The red traces correspond
to the largest single load step (85% step) in the train dataset.

left side of Fig. 3.7). The generator voltage set points are randomly sampled from a uniform

distribution [0.96, 1.04], the generator active power references are randomly sampled from [0.0, 1.0],

and the load active and reactive powers are scaled by a factor randomly sampled from [0.5, 1.5]

as in [72]. The perturbations for the base study consist of step load changes wherein the active

and reactive power of a randomly selected load are both scaled by a factor sampled from [0.0, 2.0].
 63

Each simulation is run for a total of 10 seconds with the perturbation applied at 1 second. The

recorded states are the real and imaginary voltage and real and imaginary current at the point of

common coupling between the portion of the system which will be modeled with the surrogate and

the remainder of the system. For each trajectory, 100 data points are saved at a linear spacing of

0.1 seconds (m = 100). This sampling interval is distinct from the solver time step used to generate

the datasets and evaluate the surrogates. If necessary for a particular application, the sampling

interval can be optimized as an additional hyper-parameter. The trajectories of the train dataset

are shown in Fig. 3.8.

[Link] Physics-Based Surrogate

The performance of the proposed surrogate is compared against benchmark physics-based

surrogate models. The physics-based surrogates are based on the common practice of aggregating

multiple similar devices by representing the collective behavior with a single model. This general

methodology is similar to the aggregations proposed in [61, 65] and [66], however the exact devices

vary. In our case study, the portion of the system to be replaced with a surrogate contains two types

of inverters and one type of load, therefore the most complex physics-based surrogate is comprised

of two aggregate inverter models and one aggregate load model located at a single bus (the dynamic

equations of the aggregate models and the individual device models are the same), and is shown

in gray in Fig. 3.7. The initial parameters of this model are determined by a base-power weighted

average of the same parameter for the individual devices. The base-power of the aggregate devices

is the sum of the individual base powers. While this set of parameters is a reasonable starting point,

there is no guarantee that this particular parameterization best captures the behavior of the true

model across the test dataset. Therefore, in order to make a more fair comparison, the physics-based

surrogate model is also trained using the Adam algorithm on the same training dataset. If all of

the physics-based parameters are included in the training, the underlying model is easily perturbed

to an unstable parameterization, therefore only a key subset of the parameters are learned. These

include the proportion of load modeled as constant impedance/constant power/constant current,

 64
Table 3.2: Summary of tunable hyperparameters considered for data-driven and physics-based
surrogates.

Surrogate Lower Upper Chosen

Parameter Description
type limit limit value
ninit number of hidden layers in NN1 1 3 2
winit neurons in hidden layers of NN1 3 23 17
ndyn number of hidden layers in NN2 1 3 1
wdyn neurons in hidden layers of NN2 6 39 16
data-driven
xdim hidden states 2 15 14
log10 (η) log of ADAM learning rate -3.0 -1.0 -1.52
α scaling of init. and dyn. loss 0.1 0.9 0.76
physics-based log10 (η) log of ADAM learning rate -7.0 -1.0 -7.0

the proportional and integral gains for active and reactive power control in the grid-following

inverter, and the active and reactive power droop coefficients of the grid-forming inverter. In a

realistic power system modeling scenario, the full knowledge of the structure and parameters of

the true model might not be known. Therefore, for completeness we also include surrogate models

based on aggregations of sub-sets of the devices (e.g. Load + GFL inverter, GFM inverter only,

etc.).

[Link] Hyperparameter Tuning

Hyperparameter tuning is conducted using a basic random search algorithm in order to

explore a range of options as summarized in Table 3.2. Most of the relevant hyperparameters

are related to the size of the neural networks within the data-driven surrogate. The structure of

the physics-based surrogate does not have any hyperparameters, however the initial step norm for

the Adam algorithm is also included as a tunable hyperparameter. For the data-driven surrogate

a total of 200 random hyperparameter combinations are trained and the top performing set of

hyperparameters on the validation set are selected. Given that the physics-based model has only

one hyperparameter, only 33 options are evaluated.

 65

Load + GFL + GFM

Load + GFL + GFM (trained)
Load + GFM
Load + GFL
GFL
GFM

1e-1
MAE (p.u. current)

1e-2

1e-3

1e-4

real current imag. current

Figure 3.9: MAE for the various permutations of physics-based surrogates compared to the full
order model. The error is plotted on a log scale, and each point represents the mean error for a
single perturbation trajectory (90 data points total).

[Link] Accuracy and Simulation Time Results

After the hyperparameters are chosen, the surrogates are finally evaluated on the test dataset,

which consists of 90 trajectories generated by the same process used for the train and validation

datasets. Fig. 3.9 shows the mean absolute error (MAE) of real and imaginary current for the

various physics-based surrogates for each trajectory in the test dataset. The training process results

in negligible improvements for the physics-based surrogate, indicating that the starting parameters

capture the behavior well and that the structure of the model limits the ability to improve the

performance across the entire dataset through training. As expected, the more complex surrogate

model which includes all of the devices performs best. By only including a single inverter, however,

the average error in the real current increases by almost an order of magnitude. In the results that

follow, we include the model with all three devices along with the model of only the GFM inverter
 66

as points of comparison for the data-driven surrogate.

All of the subsequent results are contingent on the encoding of the initial network inter-

face requirements in the surrogate design, which allows for running system-level simulations with

the data-driven surrogate included. Fig. 3.10 compares a single sample trajectory from the test

dataset. For this trajectory, the oscillations of the ground truth data are captured with both types

of surrogates, however, the imaginary current error is significantly lower for the data-driven sur-

rogate and the GFM only surrogate is significantly under damped. This single trajectory is not

representative of the performance overall; for more meaningful results, the surrogates are compared

on the basis of simulation speed (Fig. 3.11) and error relative to the full order detailed model (Fig.

3.12) across the entire test dataset. The level of accuracy or error that is “acceptable” is specific

to the downstream analysis enabled by the simulation. In this work, we demonstrate the accuracy

of the surrogate compared to the physics-based alternative without considering a specific workflow;

expanding the metrics and quantifying the quality of the model based on various use cases is left

to future work.

Fig. 3.11 shows the distribution of simulation times for the test dataset. The results are

straightforward and expected—the simulation time is dramatically reduced for each of the surrogate

models due to the significant reduction in total number of states. The exact speedup is dependent

on the properties of the system model, including, importantly, the stiffness and the type of solver

used. In this case study, the speed up is approximately a factor of four for the data-driven and

physics-based surrogates. One of the benefits of the proposed approach, to learn only part of a

larger system, is the scalability. For example, learning and integrating many relatively small sub-

systems (e.g. individual distribution systems) into a single larger simulation offers the potential to

scale to very large systems in a modular manner.

Fig. 3.12 shows the mean absolute error (MAE) of real and imaginary current for each

trajectory in the test dataset. The overall accuracy of the data-driven and physics-based surrogate

is similar; the data-driven surrogate has slightly lower average error for real current and slightly

higher average error for imaginary current. In contrast to the physics-based surrogate, which
 67

Data-driven Surrogate
Load + GFL + GFM
GFM
Ground Truth

Real current (p.u.) −1

−1.05

0 5 10
Imag. current (p.u.)

−0.38
−0.4
−0.42
−0.44
0 5 10
Time (s)

Figure 3.10: A single sample trajectory from the test dataset. The surrogate trajectories are the
results of numerical simulations that are run with the corresponding surrogates re-integrated as
a device component. The MAEs corresponding to this trajectory are highlighted in black in Fig.
3.12.

starts with almost optimal parameterization, the data-driven surrogate requires training to give a

reasonable output, but is able to learn the behavior without imposing a model structure a priori.

Table 3.3 compares the full order model with the two surrogates on the basis of number of

Table 3.3: Comparison of model complexity in terms of number of states and number of parameters.

Model States Parameters *Simulation time (s)

(speedup)

Full order model 432 ˜600 5.86

Physics-based surrogate 31 54 1.28 (4.58x)
Data-driven surrogate 14 1322 1.23 (4.76x)
*Average simulation time across the test dataset (90 trajectories).
 68

Data-driven Surrogate
Load + GFL + GFM
GFM
Ground Truth

Simulation time (s)

10

Figure 3.11: Simulation time comparison for the system of Fig. 3.7 with the full order model and
each of the surrogate models. The times are plotted on a log scale, and the systems with surrogate
models are approximately four times faster on average.

Data-driven Surrogate
Load + GFL + GFM
GFM
MAE (p.u. current)

1e-1

1e-2

1e-3

1e-4
real current imag. current

Figure 3.12: MAE for the data-driven and physics-based surrogates compared to the full order
model. The error is plotted on a log scale, and each point represents the mean error for a single
perturbation trajectory (90 data points total). The black dots correspond to the trajectory shown
in Fig. 3.10.

parameters, number of states, and simulation time for the test dataset. While the physics-based

surrogate dramatically reduces the number of states and parameters, the data-driven surrogate has
 69

Data-driven Surrogate
Load + GFL + GFM
GFM

MAE (p.u. current)

1e-1

1e-2

real current imag. current

Figure 3.13: MAE for the out-of-sample test dataset. The perturbation for this dataset is tripping
of half of the line connecting the surrogate with the remainder of the system.

more parameters than the original model. The expressiveness of ML models in general is related

to the large number of trainable parameters [130]; in fact, neural networks of sufficient size are a

class of universal function approximators [40]. In the context of modeling dynamic systems, the

flexibility and expressiveness of the data-driven model comes largely from the parameters and the

representation of non-linear relationships between states rather than a large number of dynamic

states. This distinction is key for reducing simulation time because the number of states largely

determines the simulation time, whereas the number of parameters has a minimal impact.

[Link] Out-Of-Sample Results

In the results presented thus far, the train, validation, and test datasets are derived from

the same processes for randomly changing the system operating point and applying a random load

step. In order to test the generalization of the surrogate models, the models are also evaluated

on a second out-of-sample test dataset. In this dataset, the operating point selection remains

unchanged and the system is perturbed by tripping one half of the connecting line between the

surrogate and the rest of the system. The error results are shown in Fig. 3.13. In this case, the
 70

data-driven surrogate outperforms the physics-based surrogate in terms of imaginary current error.

This result suggests that data-driven surrogate models with a large number of parameters can

potentially generalize across a larger number of use cases compared to physics-based surrogates,

which often require frequent re-parameterization [64, 118]. More work is needed to fully evaluate

how the surrogate generalizes to unseen perturbations.

It is important to note that the performance of physics-based and data-driven surrogates will

be highly dependent on the underlying model, therefore one method will not always outperform

the other. For example, for modeling many devices that are electrically close and very similar

(e.g. individual wind turbines of a wind plant), physics-based aggregation is likely to provide

sufficiently accurate results. For systems with a diversity of devices with varying control structures

and parameters, such as are expected in the transition towards a carbon-free power system, the

flexibility of data-driven approaches can be beneficial.

3.5 Conclusion

In this paper we define a set of initial network interface requirements and propose a data-

driven surrogate structure which satisfies the requirements by design—ensuring that the trained

models can be re-integrated into system-level simulations. The surrogate methodology is evaluated

on a medium-sized power system with 719 states, and achieves similar accuracy (mean average

current error of 8.3 × 10−3 ) and speedup (4 times) as a physics-based surrogate. However, unlike

the physics-based alternative, the proposed surrogate does not rely on detailed knowledge of the

underlying device models, which is unlikely to be available for some applications. One avenue

of future work is to apply the proposed methodology to distribution system modeling with high

levels of DERs where incomplete network knowledge is particularly relevant. For other potential

applications, such as representing multiple boundary buses surrounding an area of interest, an

extension to multi-port surrogates would be necessary.

One of the downsides of data-driven methods as compared to physics-based approaches is the

need for training, which can be a computationally expensive process. In this work, ground truth
 71

training data is generated through simulations of power system perturbations. An alternative

approach is to define an exogenous perturbation and simulate the response of only the surrogate

portion of the system. Such an approach can potentially improve the training efficiency; however,

the exact form of the perturbing signal, such that the surrogate is accurate when exposed to real

power system perturbations, is an open research question.

Another key area of future work is elucidating the exact requirements for how general a

data-driven surrogate must be in order to replace physics-based models in certain applications.

For example, in this work the data-driven surrogate is shown to perform well across different

operating conditions and types of perturbations. However, the surrogate cannot capture changes

due to internal changes to the surrogate structure or the parameters of the ground truth model.

Establishing benchmark problems will be key to determining the aspects of realistic performance

that are needed to make data-driven approaches a practical alternative in the future.

Finally, while this paper compares a purely data-driven and purely physics-based surrogate,

future work should also consider hybrid models which can take advantage of the benefits of each

method. The potential benefits of hybrid models include reduced data requirements for training,

improved generalization, and improved interpretability [130].

 Chapter 4

Integrating Data-Driven Dynamic Power System Models with Delay

Differential Equations

This chapter presents the contents of the following manuscript:

Matthew Bossart, Spyros Chatzivasileiadis, and Bri-Mathias Hodge, “Integrating Data-Driven

Dynamic Power System Models with Delay Differential Equations,” under review.

4.1 Introduction

Power systems throughout the world are rapidly introducing IBRs, motivated by the need to

decarbonize the system to meet climate targets. These additions fundamentally change the dynamic

behavior of the power system and simultaneously present challenges for how we model and analyze

the dynamic stability of the system. Given the immense scale and complexity of interconnected

power systems, one of the primary challenges is the computational burden of dynamic analyses

for large system models. In this work, we focus on transient simulations, a ubiquitous tool in

both research and operational settings for studying the non-linear behavior of large systems in

response to disturbances. The computational burden of transient simulations is exacerbated by an

increase in the total number of devices as numerous small DERs, such as battery storage, rooftop

PV, and electric vehicles, displace large, centralized power plants. Furthermore, as more IBRs

with fast control loops are added to the system, the simplifying assumptions based on separation

of timescales, which have historically allowed for very fast simulation techniques, are no longer

necessarily valid [14, 15]. In response, many techniques have been proposed to help accelerate
 73

transient simulations including aggregate reduced-order models [65, 66], parallel computing [27, 28],

dedicated hardware [29, 30], and novel numerical solvers [24, 25].

Data-driven models have been proposed to address many of the current challenges in power

system dynamic simulation, both those related to simulation speed and otherwise. For example,

data-driven models have been proposed as a means of developing both accurate and computationally

efficient model aggregations [22]. Once trained, the inference of ML models1 is typically much faster

than solving large systems of DEs, offering a pathway to simulation acceleration. Data-driven

models such as neural networks can capture complex non-linear behavior and therefore are good

candidates for representing an aggregated collection of devices. Additionally, the large number of

parameters makes it possible to represent behaviors over many operating points, avoiding the need

for excessive reparameterization [64, 118]. ML models are promising for applications where detailed

models are not known due to lack of system visibility or are otherwise unavailable, as the models

can be trained and updated directly from data. ML models can also offer potential advantages

over traditional models such as quantifying uncertainty [77] or protecting the privacy of models by

obscuring the details via a black-box model [131]. Given the wide scope of potential use-cases and

advantages of ML models, it is important to consider how these models can be incorporated into

simulations without restricting the methods for solving such systems.

4.1.1 Scope and Literature Review

In this work, we propose a methodology for integrating data-driven models into transient

simulations; the extent of the scope is detailed in Fig. 4.1 and Fig. 4.2 and described below.

Some prior works have used ML models to completely replace the transient simulation processes

[51, 50] (top row of Fig. 4.1). In contrast, we pursue a modular approach wherein data-driven

models are combined with existing physics-based models and used within the existing simulation

framework (bottom row of Fig. 4.1). The modular approach has two major advantages. First,
1
Machine learning models are a promising subset of data-driven models. In this work, we use these terms
interchangeably.
 74

Time Replace
System
Data-Driven Model Domain transient
Model
Solution simulation

Transient Simulation
System
Include
Model Solve Initialize Time
Numerical data-
Power Dynamic Domain driven
Data-Driven Integration
Flow Devices Solution models
Model

Focus of this work

Figure 4.1: Two different modes of using data-driven models for dynamic analysis. One option is to
replace the entire transient simulation with a data-driven model and predict the solution directly
(top row). In contrast, this work focuses on integrating data-driven models into the transient
simulation pipeline (bottom row).

integrating data-driven models within an existing trusted tool (i.e., numerical simulation), rather

than a complete replacement, offers a more realistic path to adoption. Power system engineers

can continue to leverage the interpretability of existing models while using using novel data-driven

models according to the specific application and engineering judgement. Second, the modular

nature of the approach allows for the potential to scale to large systems by partitioning the learning

task among many independent models. The “full replacement” approach is limited to small systems

due to the increasing data requirements and complexity of the learning task as the system size

increases [132].

In [67], the authors propose a categorization of dynamic data-driven models according to

the quantity that is predicted by the model. The data-driven model either predicts the derivative

of states and is used in conjunction with a numerical solver to generate time-domain trajectories

(referred to henceforth as time-stepping models) or directly predicts the states of the system at some

future time (referred to henceforth as direct-solution models). A visual summary of this distinction

is shown in Fig. 4.2. A typical power system time-domain simulation consists of both a system
 75

“Time-stepping
models”

“Direct-solution
models”

Focus of this work

Figure 4.2: Two different types of data-driven models. Time-stepping models are used in con-
junction with a numerical solver to advance the solution in time (top row). Our work focuses on
direct-solution models that directly predict the solution at a future time step (bottom row).

model of DEs and a time-stepping algorithm. This means including time-stepping data-driven

models is relatively straightforward as the model output (the derivative of the states) conforms

with the existing framework. In contrast, the integration of direct-solution models with existing

time-stepping algorithms is not straightforward.

There are prior works that integrate data-driven direct-solution models within a larger system

(i.e. match the focus of this work as presented in Fig. 4.1 and Fig. 4.2), however, the method of

integrating such models limits the overall flexibility and usefulness of the resulting simulation. Many

of these works are proposed for modeling active distribution systems [74, 76, 77], although some treat

the problem more generically as representing any sub-portion of a network [75]. While the specific

features of the models vary for these works, each proposes fixed time step solution methods to

match the prediction interval of the data-driven models. This coupling of the modeling and solution

methods significantly constrains the solution techniques that are available for solving the resulting

system. If a data-driven model is used as a surrogate, a computationally cheap approximation of

a detailed model, this strategy eliminates the ability to use adaptive time-stepping methods that

could otherwise be a key avenue toward acceleration [26]. In [131], the authors propose a distinct
 76

method for integrating data-driven models, specifically PINNs, into conventional numerical solvers

by approximating part of the numerical work of the solver using a trained model via “hybrid”

solvers. This approach does offer the potential for using adaptive time-stepping solvers, however

the results are limited to fixed time steps.

4.1.2 Contribution

The key contribution of this work is a methodology for incorporating direct-solution data-

driven models (Fig. 4.2) that is compatible with variable-timestep numerical solvers. Although

historically underutilized in power systems, the use of modern variable-timestep solvers offers a po-

tential path to acceleration. Variable-timestep solvers can efficiently handle the multiple timescales

and stiff dynamics introduced by fast-acting inverter controls [11, 26], adapting to rapid changes

with small time steps and taking larger steps during other periods, thereby improving the compu-

tational efficiency. While prior approaches for incorporating direct-solution models either constrain

the solver to a fixed time step ([74, 76, 77, 75]) or require the use of novel hybrid solvers ([131]),

our approach encodes the model as a delay relationship that is valid at any time and therefore does

not explicitly constrain the solution method or time step. This feature allows for using adaptive

time-stepping methods and cohesively including models with distinct prediction intervals in the

same simulation. The proposed approach is demonstrated on a collection of case studies derived

from the well known two-area test system [60] as a means of highlighting the advantages compared

to prior approaches.

We note that the focus of this work is on the specific methodology for integrating data-driven

models into simulations via delay relationships. We expect that such a method can unlock a series of

opportunities for the use of data-driven methods. This can lead to a new field of research, increased

interest by industry, and faster adoption of data-driven models in practice. We demonstrate the

methodology using relatively simple black-box models and small datasets. Future research on more

sophisticated data-driven models (including incorporating physics-based knowledge) and training

processes for large, diverse datasets are necessary to make data-driven models a viable modeling
 77

tool in practice.

4.2 Methodology

4.2.1 Power System Time-Domain Simulation

A time-domain simulation requires the specification of two components: (1) a system model,

including the DAEs that describe the system’s physics and controls, and (2) a time-stepping (or

integration) algorithm. In this work, we assume the system model consists of a set of time-invariant

DEs:
dx
M = f (x(t))
dt (4.1)
(0)
x = x(0)

where x are the states of the system (both dynamic and algebraic combined), f (·) represents the

equations of the system model, and M is a diagonal mass matrix, which is used to encode the time

constants of the dynamic states (M [i, i] = 0 if the ith state is algebraic). Note that the states x in

(4.1) include both the states of the devices (e.g. inverters, machines, loads) and the states of the

network itself.

Given an initial condition x(0) , a time-stepping algorithm is used to generate time-domain

trajectories of M points for each of the states {x(0) , ..., x(M ) }. In general, the time points of the

trajectories, {t(0) , ..., t(M ) }, will depend on both the system dynamics and the specific solver and

can have variable step-sizes for modern implicit solvers.

4.2.2 System Structure

The goal of this work is to replace a portion of the power system model in (4.1) with a

direct-solution surrogate model. To do so, we must make certain assumptions about the structure

of the original system. We assume that the system can be split into two sub-sections a and b:

dxa
Ma = fa (xa , xba )
dt (4.2)
dxb
Mb = fb (xb , xab )
dt
 78

Figure 4.3: The assumed system structure for incorporating a surrogate model.

where xa and xb are the subset of states that only appear in the equations for sections a and b

respectively, xab are the subset of states in section a that also appear in the equations for section

b, and xba are the subset of states in section b that also appear in the equations for section a. The

structure of the Jacobian for such a system is shown in Fig. 4.3 where the shaded areas indicate

areas with potentially non-zero entries and the dark shading indicates the portion of the Jacobian

representing the interface between the two sections. We note that this assumption is similar to

Assumption 1 in [131], however we do not make any assumptions about the division of algebraic

and dynamic states. Fig. 4.3 shows the diagonal of the mass matrix, which can be non-zero

anywhere, indicating that any of the system states can be considered dynamic or algebraic. As in

[131], we note that it is possible to separate the system into many sections by applying the division
 79
i12 → i23 →
... ...
v1 v2 v3

Interfacing
Section a Section b
States

Figure 4.4: An example of splitting the system at a single bus of common coupling.

iteratively as long as the assumed structure continues to hold. In practice, many of the benefits

of including data-driven models will be realized by including many distinct models within a single

system.

One intuitive way to divide a power system that is consistent with this assumption is at the

point of common coupling between two otherwise isolated areas. This scenario is shown in Fig. 4.4,

where the voltage at the connecting bus and the current injected into that bus from the adjacent

section are the interfacing states. Given the sparse interconnected structure of the power system, it

is not uncommon to be able to separate the system in this way. For example, individual distribution

systems are typically connected to the transmission system through a single point of coupling at a

substation. While we demonstrate this type of division in the numerical results, the methodology

is not limited to divisions along a single physical point of connection, as long as (4.2) is satisfied.

4.2.3 Integrating Surrogate Models with Delay Equations

The focus of this work is a method for incorporating data-driven surrogate models into the

system model shown in (4.1). The general strategy is to add delay constraint equations to the

original system thereby converting the original system to a system of DDEs. As a reminder, this

strategy enables the use of variable-timestep solvers by formulating the model independent of the

solution method. We follow the notation developed in Section 4.2.2 assuming that we will replace
 80
Table 4.1: Summary of differences between proposed direct-solution surrogate models.

Surrogate Training Method N xs

Physics-informed Gradient Optimization +

1 xa
Neural Network Physics Regularization
Terminal Data Gradient Optimization
>1 xsb
Surrogate

section a with a surrogate model. The surrogate model is given by:

xs (t) = S(xba (t), xba (t − τ ), ... , xba (t − N τ ),

xs (t − τ ), ... , xs (t − N τ )) (4.3)

xab (t) = s(xs (t))

where xs are the states of the surrogate model, S is the data-driven surrogate model, τ is a constant

time delay, N is the number of delayed inputs to the model, and s is the mapping from the surrogate

state space to the states required to interface with the remainder of the system. It is possible for

the model to predict the value of the states at time t based solely on the value at time t−τ (N = 1).

Depending on the details of the model, additional inputs at prior multiples of τ may improve the

model performance (N > 1). The surrogate model is fundamentally different than the DEs that

it replaces in that it is an algebraic relationship that advances the solution in time, as opposed to

numerical integration.

Many of the surrogates that have been proposed in the literature fit the general structure of

(4.3). For example, if s = I (the identity transformation), then the surrogate model predicts only

the interfacing states (the injected current in the example of Fig. 4.4). In this work, we refer to

this as a terminal data surrogate model because the inputs and outputs are only the interfacing

quantities at the terminus of the surrogate area. If, however, xs = xa , the surrogate model predicts

the entire state-space of the portion of the system to be replaced. PINNs are an example of this

type of surrogate wherein the underlying dynamic equations are used during the training process

to penalize deviations from the known physics [78, 79]. In this case, s is known from the underlying

dynamic model. These differences are summarized in Table 4.1.

 81

In this work, we use a terminal data surrogate in the numerical case study. This choice is

not meant to suggest that a black-box terminal data model is preferable to PINNs, but rather that

the model is a better choice to thoroughly compare the numerical properties of a power system

model with and without the surrogate model included. One of the potential reasons to use a direct-

solution surrogate is to reduce the size of the state space, which is not possible with PINNs. PINNs

are advantageous due to their improved training performance and interpretability compared to a

fully black-box model, however, the training and modeling performance, while important, are not

the focus of this work.

Regardless of the exact form of the surrogate, we can include the surrogate from (4.3) in the

physics-based formulation as an algebraic constraint.

dxb
Mb = fb (xb (t), s(xs (t))
dt
0 = xs (t) − S(xba (t), xba (t − τ ), ... , xba (t − N τ ), (4.4)

xs (t − τ ), ... , xs (t − N τ ))

Then, the entire system can be written as a set of DDEs:

dx̂
M̂ = fˆ(x̂(t), x̂(t − τ ), ... , x̂(t − N τ ))
dt (4.5)
x̂ = h(t) ∀ t < Nτ

where M̂ , fˆ(·), and x̂ are the mass matrix, equations, and states respectively of the system with

the surrogate model. The system in (4.5) is meant to approximate the ground truth system of

(4.1). Because (4.5) is a DDE (i.e. the derivatives can depend on delayed values of the states), an

initial history function, h, is required in place of an initial condition [81].

4.2.4 Solving Delay Differential Equations

In this section we discuss various considerations for solving the system of DDEs in (4.5).
 82

[Link] Method of Steps

The most popular method for solving DDEs, and the method used in this work, is the method

of steps. The basic idea of the method is to break down a system of delay equations as a sequence

of initial value problems (IVPs) where the size of the delays determine the size of the interval of

the component problems [81]. For simplicity of notation, consider the case of (4.5) with a single

constant delay (N = 1). An IVP for the kth interval of the solution is given by:

dx̂k
M̂ = fˆ(x̂k (t), x̂k−1 (t − τ )), t ∈ [tk−1 , tk ] (4.6a)
dt
x̂k (tk−1 ) = x̂k−1 (tk−1 ) (4.6b)

where tk := kτ , x̂k is the solution for the kth interval, and x̂k−1 is the previously computed solution

for the prior interval. Equation (4.6b) indicates that the solution for adjoining IVPs is continuous

(i.e. the end point of the prior solution segment equals the starting point for the current segment).

The full solution is defined by the solutions to the constituent IVPs:

x(t) := x̂k (t) ⇐⇒ t ∈ [tk , tk+1 ] (4.7)

For the first interval, the delayed values for t < 0 are given by the history function, and subsequently

the delayed values come from prior solution intervals. A key feature of the method of steps is

that results from the theory of DEs (e.g. guarantees for existence and uniqueness of solutions)

in addition to numerical solution techniques are transferable and applicable to DDEs. This is a

critical advantage of this work; while we require additional complexity we retain access to many

of the advanced features available from the vast body of literature and software for DEs, such as

modern adaptive time-stepping solvers.

One of the numerical challenges that is unique to solving DDEs is the propagation of discon-

tinuities [133]. It is possible that the overall solution can have discontinuities even if the history

function and the individual IVP solutions are smooth, arising from the transition between the two.

In practice, the initial conditions in (4.1) are found by solving the power flow equations and initial-

izing each of the individual device models such that the entire system is in a steady state condition.
 83

This means that the history function of (4.5) will typically be a constant function:

h(t) = x̂(0) ∀t < N τ (4.8)

Provided that the system is not perturbed from the steady state operating condition before time

τ N , the transition between the history function and the solution intervals will be smooth and the

problem of propagation of discontinuities can be avoided. While we consider surrogates with fixed

delays, we note that an extension to surrogate models with time-dependent or state-dependent

delays is possible given that there are solution methods for these types of problems [134].

[Link] Comparison with solving DEs

In this work we emphasize computational acceleration as a key motivating feature for employ-

ing surrogate models. Of course, there are additional benefits of employing data-driven surrogate

models that are not solely based on acceleration. For example, data-driven models can be bene-

ficial when physics-based models are unavailable or unknown, or as a means to updating models

based on measured data. However, we believe that given the current strain to conduct transient

simulations for large power systems with sufficient accuracy, data-driven models will not be useful

if they significantly increase the computational burden. Therefore, in this section we describe the

trade-offs that impact the computational burden of solving the original ground truth system in

(4.1) versus the system with a surrogate model of (4.5).

By converting from a system of ODEs or DAEs to DDEs, we inherently add some computa-

tional overhead related to storing and accessing past values of the system state. In addition, the

solution of DDEs requires a solution time step at repeating intervals of the delayed values as part

of the method of steps. This means that if everything else is equal, the delayed system will likely

require more time steps. For neural network-based models, there is additional overhead related

to the forward evaluation of the model. This overhead can be significant as the model will be

evaluated many times during the solution process.

In order to overcome the computational disadvantages described above, there are two po-
 84

tential mechanisms for speedup by including direct-solution surrogates. First, is reduction of the

the overall system order. In (4.1), dim(x) = dim(xa ) + dim(xb ), whereas in (4.5), dim(x̂) =

dim(xb ) + dim(xs ). Therefore this approach is more likely to be computationally beneficial if

the dimension of the surrogate is much less than the portion of the system that it replaces

(dim(xs ) << dim(xa )). This condition can potentially be achieved by using multiple independent

surrogate models in a single system (e.g. representing many distribution grids with surrogates).

The second mechanism for acceleration is the conversion of dynamic states represented by DEs to

an algebraic relationship defined by the surrogate model. This conversion can be especially bene-

ficial to the overall simulation time if the dynamic states that are eliminated are constraining the

time step of the numerical solver. This is the mechanism for acceleration that is demonstrated in

[131]. In Section 4.3.4 we characterize the trade-offs described above to determine the conditions

for achieving simulation acceleration.

4.2.5 Computational Setup

The computational set up for generating numerical results is summarized in Fig. 4.5 and

described in this section. The overarching goal is to compare the solution characteristics of the

original system with a modified system that includes the direct-solution surrogate (shown in orange

and gray respectively).

[Link] Surrogate Model

We assume a terminal data surrogate model where the interfacing states are the current and

voltage at the bus separating the two sub-systems (see Fig. 4.4). Before the terminal quantities

are passed to the neural network model, the input quantities are normalized for a mean of zero and

standard deviation of one as is common in many ML workflows. The core part of the surrogate

model can take on a variety of forms. We choose long short-term memory (LSTM) networks

[135] as the model to process the delayed inputs and predict the output. Unlike feed-forward

neural networks (FFNNs), which feature unidirectional flow of information from inputs to outputs,
 85
[Link] [Link]
Differential
Run
Equation
Simulation
System
(Solve)
Model

Surrogate
Training
[Link] Update [Link]
params.
Terminal Delay
Run
Data Differential
Surrogate + Equation
Simulation
(Solve)
Model Model

Figure 4.5: High-level block diagram of the computational setup to generate the numerical results.
The colors of the blocks indicate the primary software package used in that portion.

LSTMs include feedback loops, which makes them well suited to deal with sequential time series

prediction tasks [136]. Based on initial experience, LSTMs outperform FFNNs and tend to remain

numerically stable throughout the training process.

[Link] Training

The datasets for training the surrogate models are generated from running a series of M

dynamic simulations of the full ground truth system of DEs:

m −1 M
D = {{t(n) , v (n) , i(n) }N
n=0 }m=1 (4.9)

where v = [v r , v i ], i = [ir , ii ], M is the number of simulations in the dataset, and Nm is the

number of timesteps per simulation (n is the index of the data point within a simulation). The

total number of training data points is given by |D| = M

P
m=1 Nm . For the numerical results, we

train the surrogates on a single perturbation (M = 1).

The training process consists of iteratively updating the parameters of the model to minimize

a scalar loss function. In the case study, we use the popular Adam algorithm [126]. For the terminal

data surrogate, the input is the voltage at the point of coupling and the output is the injected
 86

current. We choose the loss to be the RMS error of the simulated current in the system reference

frame with respect to the ground truth data:

v
|D|
u
u 1 X r
u 
L= t r 2 i i
(ij − îj ) + (ij − îj ) 2 (4.10)
2|D|
j=1

One option for generating the predictions, î, is to apply the surrogate model in isolation to the

values of the sampled dataset. In this case the model is only trained at a discrete set of points and

is likely to perform poorly when used in conjunction with a variable-step solver, as many of the

evaluation time steps will not correspond to the sampled times of the training dataset. Instead, we

emphasize that î comes from the solution trajectory of the DDE system. Intuitively, this will benefit

the performance as the surrogate is trained and deployed in the same environment (i.e. within a

numerical simulation). For this reason, training the model requires back-propagation through the

numerical solution of the model. We use a discrete reverse mode sensitivity method due to its

numerical stability and scaling with respect to the number of parameters [112, 84].

[Link] Software

The entire software stack for generating the numerical results is written in the Julia pro-

gramming language. Julia is chosen for this work for its unique blend of high performance and

ease of use. [Link] is chosen for the power systems modeling [26]. As part of this work, [Link]

was extended to support delays. In addition, [Link] has experimental support for computing the

gradient of a transient simulation with respect to the system parameters. This feature allows the

training to take place within the simulation environment by iteratively updating the parameters

that define the surrogate model while keeping the remaining system parameters fixed. This is a sig-

nificant improvement over more complex workflows that require isolated training and re-integration

into the simulation environment ([77, 131, 74, 76, 75, 22, 72]).

The numerical solution of the models is handled by the packages of the DifferentialEqua-

[Link] ecosystem [25], specifically [Link] and [Link]. These packages

include many state of the art equation solvers, including the ability to automatically compose many
 87

solvers with the method of steps to extend to DDEs [133]. The surrogate model uses [Link] [137] for

the neural network implementations, however the model is packaged as a component corresponding

to the [Link] API so that it can be seamlessly added to the system model in [Link].

One of the key features of the computational setup is the composability of the software

packages shown in Fig. 4.5. Each of the main functions of the methodology (e.g. power systems

component modeling, neural network implementations, optimization algorithms, solving DEs, etc.)

is handled by a dedicated package that can be composed for this application. This is critical to

ensure that this work can be extended and improved as ML models become more commonplace in

power systems applications and as new solvers are developed.

4.3 Numerical Results

4.3.1 Test System

The system used (with various modifications) throughout this section is the well known two-

area test system [60] shown in Fig. 4.62 . This system is commonly used to study inter-area

oscillations between two distinct areas connected by relatively long transmission lines. The two

areas interact via the oscillations, therefore this is a compelling case to demonstrate a direct-
2
The code and parameters for all numerical results are available online: [Link]
bossart/PowerSystemDelaySurrogates

Area 1 Area 2

1 3

2 4

Figure 4.6: Single line diagram of the 2-area system used in the case study. The first two results
demonstrate replacing portions of the system with a surrogate model.
 88

î(t)

: direct-solution
surrogate

Figure 4.7: A DDE system with a single direct-solution surrogate replacing a full area of the original
system. The labeled current is compared in Fig. 4.8.

solution surrogate that interacts with the DEs of the unmodified portion of the system and can

capture the same behavior as the ground truth system.

4.3.2 Single Surrogate Model

In this section we demonstrate replacing area 1 of the two-area system with a direct-solution

surrogate model. We train and evaluate the dynamic performance of the model in response to a

25% positive step in the active power setpoint of generator 3 (the starting active power output

is 700 MW). There are two important caveats to emphasize. First, the faulted generator exists

outside of the surrogate model. It is not possible to represent a specific fault within the surrogate

as the individual buses have been abstracted and replaced with the neural network based model.

While this is a limitation, this type of aggregation is commonplace in power systems [65, 66] and

can be useful depending on the application. Second, the computational benefit of using data-driven

surrogates is realized by first training the model on a dataset and then deploying the model for

unseen simulations. However, we do not train and evaluate the model on separate datasets as the

focus of this work is the re-integration of the models, not the model training. Robust training

over realistic sized datasets is an important avenue of future work, and can likely be improved by

incorporating physics-based knowledge into the surrogate models (e.g. via PINNs).

The time-domain trajectories of the interfacing current between the surrogate and the rest
 89
Ground truth
Voltage Source Surrogate
Constant Impedance Surrogate
Terminal Data Surrogate (pre-training)
Terminal Data Surrogate (trained)

Real current (p.u.) 4.5

3.5

3
0 2 4 6 8 10
−1
Imag. current (p.u.)

−1.1

−1.2

−1.3

−1.4
0 2 4 6 8 10
Time (s)

Figure 4.8: Comparison of current trajectories for the ground truth systems and various surrogates,
including the proposed terminal data surrogate, an ideal voltage source, and a constant impedance.

of the system are shown in Fig. 4.8 for the original and surrogate systems in black and orange

respectively. As a point of comparison, traces are also included from simulations where one area

is instead replaced with an ideal voltage source (blue) and constant impedance (green). An ideal

voltage source or infinite bus is commonly used to represent surrounding areas that are unable

to be modeled in detail. Fig. 4.8 shows that the direct-solution surrogate can be incorporated

and trained within a simulation and capture the interactions between two areas of the original

model. In contrast, the voltage source and constant impedance models fail to accurately capture

the inter-area oscillations.

 90

Figure 4.9: Work precision diagram for the DDE system shown in Fig. 4.7. The x-axis shows the
final error relative to a high accuracy reference solution of the same system and the y-axis shows
the simulation time for various solvers and tolerances. Note that both axes are in logarithmic scale.

While a result similar to Fig. 4.8 would be possible with a number of prior data-driven

models, the key contribution of this work is the flexibility of solution methods for solving the DDE

model. Fig. 4.9 shows a work-precision diagram comparing multiple solvers at various tolerances

for the same problem. This result demonstrates the automatic access to a plethora of numerical

solvers via the existing ODE solver ecosystem together with the method of steps; the result of Fig.

4.9 which includes multiple variable-timestep methods is not possible with prior approaches. With

access to variable-timestep solvers, we can importantly trade off the accuracy of the solution and

the computational burden depending on the needs of a specific application. In addition, we can

automatically take advantage of advances in numerical solvers that arise from other disciplines, as

these large-scale numerical problems are by no means unique to power systems.

4.3.3 Multiple Surrogate Models

In this section we demonstrate an extension to including multiple surrogate models within a

single simulation. In particular, we incorporate two surrogate models that are trained at different
 91
î1 (t)

î2 (t) ↑

Figure 4.10: A DDE system with two direct-solution surrogates replacing individual generators of
the original system. The labeled current is compared in Fig. 4.11.

prediction intervals. This demonstrates one of the key features of the proposed methodology; the

prediction interval of the data-driven model is not directly coupled to the time step of the numer-

ical solver. This feature allows for the interoperability of component models, which is especially

important in building large scale power system simulations where models might come from distinct

entities. On the modeling side, the prediction interval is related to the time-scale of the dynamics

of the underlying component(s), therefore it is plausible to include differing prediction intervals

when replacing different component(s).

We modify the two-area system by replacing generator 2 with a droop-based grid-forming

inverter and develop surrogates for individual devices (Fig. 4.10). For the synchronous generator

surrogate we choose τ = 500 ms, whereas for the grid-forming inverter surrogate we choose τ = 200

ms due to the faster dynamics of the inverter control. N = 4 for both models (τ and N are defined

in (4.3)). The interfacing currents for both of the surrogates are shown in Fig. 4.11. For the real

current, the two surrogates successfully capture the differing behavior of the synchronous generator

(more oscillatory) and the grid-forming inverter (less oscillatory). As in the prior section, the

terminal data surrogates can be trained to significantly outperform the voltage source and constant

impedance surrogates.
 92
Ground truth
Voltage Source Surrogate
Constant Impedance Surrogate
Terminal Data Surrogate (pre-training)
Terminal Data Surrogate (trained)

7.2
Real current (p.u.) 7
6.8
6.6
6.4
6.2
6
0 2 4 6 8 10
Imag. current (p.u.)

−1

−1.5

−2

−2.5

0 2 4 6 8 10
Time (s)

Figure 4.11: Comparison of current trajectories for the ground truth systems and various surrogates,
including the proposed terminal data surrogate, an ideal voltage source, and a constant impedance.

4.3.4 Characterizing Computational Trade-offs

In this section we numerically characterize the trade-offs discussed in Section [Link] to ex-

amine the potential for simulation acceleration. Our goal is to isolate the numerical differences

for solving the DE and DDE systems without considering the error related to imperfect training

of the data-driven models and its impact on the solution. Therefore, for the following results we

employ an ideal surrogate model where the predicted output current comes directly from a prior

solution of the DE system. The implementation is identical to the terminal data-surrogate (i.e. the

neural network is still evaluated based on the delayed state values), with the addition of an ideal
 93

Simulation time (s) 10

1
Line trip , DDE solution
Line trip , DE Solution
Power step , DDE solution
Power step , DE Solution

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Reduction in # of states (%)
Figure 4.12: Comparison of simulation time between DE and DDE systems for different reductions
in the number of states. The DE system has 356 states, and the number of states replaced by the
surrogate is varied by modeling certain devices as many parallel connected devices with scaled base
powers.

correction term to match the output from the DE system. This methodology allows us to contin-

uously modify the defining characteristics of the surrogate without a) training many independent

models and b) obscuring the underlying trade-offs due to the differing levels of error for the trained

models. All of the following results use Rodas5P [138], a fifth order Rosenbrock method, as this

solver is the highest performing in our work (Fig. 4.9) and prior power system benchmarks [26].

Each simulation is run for 10 s with the same perturbation as prior results applied at 5 s.

Fig. 4.12 compares the simulation time of a ground truth DE system with 356 states with

DDE systems with varying reductions in the number of states. The network structure of the DE

system is identical to the original two-area system (Fig. 4.6), however each generator is replaced

by multiple equivalent parallel connected generators with scaled base powers [65] to increase the

number of states and make the system more computationally expensive. The type of perturbation
 94
Line trip , DDE Solution
Line trip , DE Solution
Power step , DDE Solution
10 Power step , DE Solution

Simulation time (s)

0.1

0.001 0.01 0.1 1

Delay length, τ (s)
Figure 4.13: Comparison of simulation time between DE and DDE systems for different lengths of
the delay for the terminal data surrogate.

can significantly impact the dynamics that are excited and therefore the solution time. For each

of the results in this section, a step in the active power reference is significantly faster to solve

than a branch trip as the latter excites the fast dynamics of the network and requires many more

time steps immediately following the perturbation. For both perturbations, reducing the number

of states reduces the overall simulation time, however for the more computationally challenging line

trip, the crossover point occurs at around 10% state reduction whereas the crossover for the power

step occurs at around 25% state reduction.

Fig. 4.13 compares the simulation time for various prediction intervals of the terminal data

surrogate (i.e. different delay sizes). Recall that for the method of steps the size of the delay

impacts the number of solution intervals. For the longest delay size of 5 s, the overhead of the

method of steps is negligible and the DDE system is faster to solve for both faults.3 Once again,
3
A 5 s delay size is unrealistic to capture much faster changes in the system, and is only possible to characterize
using the ideal surrogate.
 95
10
Line trip , DDE Solution
Line trip , DE Solution
Power step , DDE Solution
Power step , DE Solution

Simulation time (s)

0.1
10 100 1000 10k 100k
Number of NN parameters
Figure 4.14: Comparison of simulation time between DE and DDE systems for different number of
parameters in the terminal data surrogate. The neural network is evaluated as part of the surrogate
model to generate these results, but the output is corrected according to the ideal model.

the crossover point when the DE and DDE systems take equal time to solve is more advantageous

for the more difficult line trip perturbation. The limiting behavior for both perturbations as the

delay size decreases converges to the same approximate simulation time as the numerical time step

becomes primarily constrained by the delay intervals rather than the step-size needed to solve the

original system.

Given the importance of the delay on the computational efficiency, the question remains on

how to choose the delay size for a particular application. Intuitively, the delay should have some

relation to the timescale of the dynamics that are desired to be captured. However, there is some

evidence that the feasible prediction horizon for PINNs can far exceed the time steps of numerical

solvers for solving the same system [139]; this is the exact condition that is necessary for our

methodology to be computationally advantageous.

Fig. 4.14 compares the simulation time for various neural network sizes. For neural networks
 96

with less than 10,000 parameters, the impact of the forward evaluation time on the total simulation

time is minimal. For much larger networks, however, the impact becomes severe. For large neural

networks, the memory requirements for storing the models can also be a limiting consideration.

The impact of converting dynamic states to algebraic relationships is not isolated in these

results. The authors in [131] demonstrate this impact, but do not consider a system with dy-

namic lines, which are often the fastest states, and are increasingly important to model with the

growth of IBRs. An interesting research direction is to develop surrogates for the network itself

that can capture dynamic phenomena that are not possible with a traditional algebraic network

representation.

4.4 Conclusions and Future Work

Data-driven models have been shown to produce high-quality reduced-order models that can

help address the computational burden of simulating inverter-based power systems. This paper

presents a modular and scalable methodology for combining the promise of data-driven models

with existing mature state-of-the-art numerical solvers. Prior approaches for using direct-solution

data-driven models in power system simulations were restricted to fixed-timestep solution methods.

In contrast, our methodology enables the use of direct-solution surrogate models within transient

power system simulations through the introduction of delay relationships to the system of DEs,

thereby converting the system of DEs to DDEs. Our method takes advantage of the wide range

of existing adaptive timestep solution methods via the method of steps. In addition, multiple

surrogates with different delays can seamlessly be introduced in the same system resulting in a

DDE system with multiple fixed delays. This paper unlocks new opportunities to use data-driven

models, in turn contributing to the extensive validation that is needed to develop trust and adoption

in industry.

There are areas of future work that are needed to make the proposed methodology opera-

tionally feasible, but are not fully addressed in this work. Data-driven models must be shown to
 97

sufficiently4 capture the behavior of the original system. This includes training models on large and

diverse datasets without getting stuck in local minima, which is a common challenge for training

dynamic models. Incorporating physics-based knowledge (e.g. PINNs) is one avenue to poten-

tially address both the modeling and training challenges, however there are many opportunities for

progress in this area.

4
The definition of sufficient is application dependent and should be quantified through application-specific bench-
mark problems
 Chapter 5

Differentiable Programming for Transient Power System Simulations

This chapter presents the contents of the following manuscript:

Matthew Bossart and Bri-Mathias Hodge, “Differentiable Programming for Transient Power

System Simulations,” under review.

5.1 Introduction

The overarching goal of the modern electric power system is to reliably deliver electrical

energy from generation resources to loads at low cost. One key aspect of reliably operating the power

system is to continually maintain dynamic stability, or the ability of the system to return to a stable

operation point following a perturbation to the system (e.g. fault, device trip, etc.). The power

system is critical infrastructure that is vital to everyday life, therefore power system researchers and

operators have minimal opportunities to run experiments on the real system. Instead, we must rely

on mathematical models to assess the stability of the system. The status quo across many important

power system analyses such as contingency analysis and interconnection studies is to conduct

numerical simulations to assess the response of a specific model with a specific parameterization.

Such analysis is necessary to understand the transient large-signal behavior of the system in response

to credible contingencies and to plan preventative actions accordingly.

The process of simulating a model of the power system can be understood generically as the

forward execution of a numerical model. In contrast, DP refers to the methods for computing the

gradient of a numerical model’s output [84]. Gradient information has been instrumental in solving
 99

inverse problems in a variety of scientific disciplines that, like power systems, rely on numerically

solving DEs. For example, gradients have been used for aerodynamic shape optimization [140] in

computational fluid dynamics, improving numerical weather prediction accuracy [141], optimizing

mechanical beam design [142], and improved ice sheet modeling [143] to name a few. This work

focuses on applying DP to power systems through the computation of gradients for dynamic sim-

ulations, however, the ultimate goal is to improve or enable key power system applications such as

parameter tuning, optimal control, and uncertainty quantification.

The idea of calculating and applying gradients in power systems is not new. The seminal

works in this area developed TSA—the sensitivity of trajectories to small changes in either pa-

rameters or initial conditions—for hybrid power system models [91]. TSA is an example of DP in

which the outputs of interest are the state trajectories and sensitivities are calculated by propagat-

ing derivatives forward while simultaneously solving the original problem. Trajectory sensitivities

have enabled numerous algorithms for solving inverse problems including parameter estimation and

quantifying uncertainty [144]. The reader is referred to [145] and [88] for a comprehensive review

of the power system applications of TSA. The main drawback of the TSA approach is the poor

scaling with the number of parameters and/or size of the system. The authors in [88] note the

scaling limitations of the TSA approach and develop a discrete adjoint method that is compatible

with events.

One of the unique features of a power system simulation, in contrast to other domains which

also utilize systems of DEs, is the process for initializing the model. In most power system ap-

plications, the user does not directly set the initial conditions of the model. Instead, the model

is initialized according to meaningful inputs (e.g. setpoints) and a predefined routine in order to

obtain a feasible and steady state operating condition [122, 26]. From a numerical perspective, the

initialization can consist of solving systems of nonlinear equations before solving the DEs. In this

work, we propose a method for differentiating a power system simulation with respect to parameters

that explicitly includes the initialization process. While prior works have developed techniques for

differentiating the solution of nonlinear equations [146] and DEs [84], including extensions to DAEs
 100

[89] and equations with discrete events [91, 88], these are all sub-problems of a complete power

system dynamic simulation. This paper builds on these prior works by applying the concepts

of differential programming holistically to the process of power system simulation, ensuring that

the gradient is calculated based on simulations which retain the desired steady-state behavior and

match the power flow solution. We demonstrate the importance of including the initialization in the

gradient calculation by evaluating the performance on a parameter tuning optimization problem

using a gradient-based algorithm. When the tuned parameter impacts the initial operating point,

including the initialization becomes critical to converging to the optimal parameterization.

In addition to more traditional power systems optimization problems, recent work has sug-

gested DP as a bridge between traditional scientific computing and ML training workflows [130].

In this vein, by differentiating through the entire simulation, our work enables the training of

data-driven models directly within a mixed simulation environment containing both traditional

physics-based models and data-driven models. Such mixed approaches including PINNs [78, 79],

universal differential equations [87], and others promise to take advantage of the expressiveness

and performance of ML techniques while retaining the interpretability and generalization of first-

principles models. While it is possible to train data-driven models in a bespoke environment and

then implement them within a simulation package for validation and testing [22, 72, 80], this siloed

process presents a significant practical impediment to the adoption of data-driven models. Instead,

we successfully train ML models using gradient-based optimization directly within an existing power

system simulation environment.

While this work emphasizes end-to-end differentiation of the simulation process, the most

performance-critical portion of the simulation is typically the solution of the DEs. Just as there are

many algorithms for numerically solving DEs, there are a variety of approaches for differentiating

programs that solve DEs [84], many of which have not yet been applied to power system problems.

Each approach can vary based on accuracy, stability, efficiency, memory requirements, and diffi-

culty of implementation. For power system simulations, the size of the system, type of simulation

(e.g. quasi-static phasor, EMT), and other modeling assumptions can vary widely. As such, the
 101

choice of method for differentiating a power system simulation is non-trivial and highly system

and application dependent; there is no single best method for all cases. This paper addresses this

challenge by introducing a modular and flexible architecture for calculating gradients, making it

possible to test multiple methods for a specific system with minimal code changes. The flexibility

is enabled by a composable software ecosystem with separate modeling and solution layers. This

capability is key for choosing the optimal method that is suited to the downstream application of

the gradients.

The remainder of the paper is organized as follows. In Section 5.2, we develop the definition

and notation for the forward transient numerical simulation. Section 5.3 defines mathematically

the gradient quantity we are interested in computing. Section 5.4 describes the software ecosystem

for computing the gradients. Section 5.5 presents numerical simulation results for a number of po-

tential power system applications that demonstrate the effectiveness and versatility of the proposed

method. Finally, the paper is concluded with a summary and discussion of future work in Section

5.6.

5.2 Background and Notation

In this section we establish a definition and the necessary mathematical notation for transient

power system simulations. Note that the definition presented applies to power system models

which employ phasor quantities to represent three-phase voltage and current and can therefore

be initialized to steady state conditions. Waveform simulations which represent the time-varying

sinusoidal network quantities as simulation states directly are not included. The notation of this

section is shown pictorially in Fig. 5.1.

We assume a generic time-invariant power system model:

dx
= F (x, y, ψ d , sd )
dt (5.1)
dy
= G(x, y, ψ n , sn )
dt
where x, F , ψ d , and sd are, respectively, the states, equations, parameters, and setpoints of
 102
Device Transient Simulation
Parameters

Network Numerical
Setpoints Power Device Objective
Flow Initialization Integration
Function
Network
Parameters

Figure 5.1: The steps of a power system transient simulation.

the devices. Similarly, y, G, ψ n , and sn are the states, equations, parameters, and setpoints of the

network. The initial states of the system are given by x0 and y 0 .

In this work, we consider a transient simulation to consist of three distinct steps:

(1) Solve the network power flow equations 0 = G(y 0 , ψ n , sn ) for the initial network states

(y 0 ).

(2) Solve the device equations 0 = F (x0 , y 0 , ψ d , sd ) for both the initial device states (x0 ) and

the device setpoints (sd ).

(3) Given the input parameters (ψ d , sn , ψ n ) along with the values found in the previous steps

(x0 , y 0 , sd ), the parameters and initial conditions of Eq. (5.1) are fully specified. The

solution is advanced in time using a numerical integration algorithm to a series of discrete

time points: xt := {x0 , ..., xn }, y t := {y 0 , ..., y n }. The model is explicitly required to start

in steady state (Step 2), therefore an event E is applied to the system during the numerical

integration in order to incur some non-stationary dynamic response.

For simplicity of notation we define a single function fsim which encapsulates the three steps

described above:

xt , y t = fsim (E, ψ n , sn , ψ d ) (5.2)

Furthermore, we define a scalar-valued objective function which takes as inputs the time-
 103

domain trajectories of the numerical simulation, the parameters, and additional data D:

L = fobj (xt , y t , ψ n , sn , ψ d , D) (5.3)

In a classical optimization framework, D consists of measurements or simulated data from a high-

fidelity model and the objective function quantifies the quality of the match between the trajectories

(e.g. root-mean-square error). For our purposes, however, fobj can take any differentiable form.

Finally, we define a top-level function as the composition of the numerical simulation function and

scalar-valued objective function:

g = fobj ◦ fsim (5.4)

5.3 Differentiating Power System Simulations

The overarching goal of this work is to develop a method to calculate the gradient of function

g—which consists of both the numerical solution of a simulation and a user-defined scalar objective

function—with respect to parameters p:

 
∂g(p)
 ∂p1 
 . 
 
∇g(p) =  .. 
(5.5)
 
 
∂g(p)
∂pn

p = [ψ Tn , sTn , ψ Td ]T ∈ Rn

Note that the parameters of Eq. (5.5) include the network setpoints (power and voltage setpoints of

generation buses, active and reactive power setpoints of load buses), the network parameters (e.g.

the impedance of network lines), and the device parameters (e.g. the inertia of a generator), but

not the device setpoints (e.g. the voltage reference of an automatic voltage regulator). In order to

maintain the feasibility of initializing the dynamic model to the power flow solution and starting in

steady state, we only consider parameters that can be known and set a priori by the model user,

not those which are determined during the initialization. The primary contribution of this paper

is a flexible framework and software implementation for the non-trivial computation of Eq. (5.5).
 104

There are multiple existing approaches that can be applied to compute the gradient of Eq.

(5.5). The finite differences approach numerically approximates the gradient via multiple forward

evaluations of the function at different values [147, 148]. This approach is straightforward to imple-

ment given the ability to execute the forward function. However, the accuracy of the approximation

is limited by the size of the step and decreasing the step size leads to floating point cancellation

errors [149]. Furthermore, finite differences scales poorly with the parameter size, as each par-

tial derivative must be computed separately. Another option is symbolic differentiation wherein

functions are represented and differentiated algebraically. The main downside of symbolic differ-

entiation is expression swell; the algebraic derivative expressions can quickly become exceedingly

large as they are composed via the chain rule for complex programs [150]. Instead, we focus on

AD which algorithmically propagates derivative information based on known derivatives for simple

functions and the chain rule. AD is well known for its role in revolutionizing the field of deep

learning [150, 151] by allowing for efficient gradient computation of models with large numbers of

parameters. Like ML applications, we focus on reverse AD—propagating derivatives from output

to input—due to its computational advantage for problems where the input dimension is much

larger than the output dimension (in our case the input dimension is the size of p and the output

dimension is one due to the scalar objective function). AD is “automatic” in the sense that a func-

tion is differentiated by simply propagating derivative information through the component steps

without any manual intervention by the user.

While it is possible to apply AD directly to the solution of numerical equations, it is not

the only option [146, 84]. In the following sections, we discuss the methods for differentiating the

numerical sub-problems of a simulation: solving nonlinear equations and solving DEs.

5.3.1 Nonlinear Equations

Solving the power flow equations and initializing the dynamic device models are both exam-

ples of solving nonlinear equations which are necessary for transient simulation. For demonstration,

consider the generic solution x of a nonlinear system of equations f with parameters p along with
 105

a scalar objective function:

0 = f (x, p)
(5.6)
L = g(x, p)

By the chain rule we can write the derivative of g with respect to the parameters:

dg ∂g ∂g ∂x
= + (5.7)
dp ∂p ∂x ∂p

Differentiating the f equation with respect to p and substituting gives:

∂f ∂f ∂x
0= + (5.8)
∂p ∂x ∂p
dg ∂g ∂g ∂f −1 ∂f
= − (5.9)
dp ∂p | ∂x{z∂x } ∂p
:=λT

The trick to efficiently calculating the gradient is to define the adjoint as shown in Eq. (5.9). Then

solving Eq.(5.7) is accomplished by solving the adjoint equation for λ, which is simply a linear

system of equations of the size of the input dimension:

∂f T ∂g T
λ= (5.10)
∂x ∂x

Note that with this adjoint method, the gradient of g with respect to all parameters is solved

in a single step, and is independent of the steps of the numerical algorithm used to solve the

original problem. For complicated problems requiring many iterations, this adjoint method for

solving nonlinear equations is more efficient than direct AD through the steps of the computation,

underscoring the value of alternative methods for differentiating solutions of numerical problems

within an AD framework.

5.3.2 Differential Equations

The options for differentiating the solution of DEs are more varied than the options for

nonlinear equations. We note that in the power systems context this choice of sensitivity method is
 106

especially consequential given that the differential equation will dominate the overall computation

time due to the scaling of numerical solvers for large and/or stiff systems [116, 152] . In this section

we survey the existing approaches with a focus on the features most relevant to power system models

(the reader is referred to [84] for a more in-depth review). The composable and modular design

and implementation of our method allows application specific testing of the methods described in

the remainder of this section.

[Link] Discrete vs. Continuous Approaches

One way to categorize the existing approaches is by the order of the differentiation and

discretization steps: Continuous approaches calculate gradients by formulating a new system of

equations (differentiation) before solving the new system numerically (discretization). In contrast,

discrete approaches operate on the original solution (discretization) and propagate derivatives of the

unit operations of the solution (differentiation). In general these two approaches are mathematically

distinct and will not give the same exact numerical result (in other words, the discretization and

differentiation steps do not commute [82, 83]). The choice of sensitivity method can depend on the

solver; in [84], the authors show a simple case where the discrete method is non-convergent in the

case where the iterative solver has adaptive step size control that depends on the parameters that are

differentiated. Other works have demonstrated numerical oscillations in the gradient calculations

for discrete approaches to partial differential equations [85, 86]. In [87], the authors generalize that

continuous methods tend to be more efficient, while discrete methods tend to be more stable [87].

[Link] Forward vs. Reverse Approaches

Similar to AD, approaches for computing the gradient of differential equation solves can also

be categorized as forward or reverse depending on the direction of propagation of the derivatives.

In fact, directly applying forward/reverse AD to the numerical solve of the original system results

in a corresponding forward/reverse discrete method. In the DEs community, reverse methods are

commonly referred to as adjoint methods (e.g. [88, 89]).

 107

The key consideration for choosing a forward or reverse method is the relative input and

output dimensions of the function being differentiated. Generally, for both continuous and dis-

crete methods, when the output dimension is larger, forward methods tend to be computationally

favorable, whereas when the input dimension is larger reverse methods are preferable [84]. The

well known TSA [91] is an example of a forward sensitivity method, therefore the method will scale

poorly with the number of parameters. This makes it incompatible with large power system models

or models that have large numbers of parameters (e.g. many modern ML models).

In addition to performance, memory requirements are another factor in choosing a method.

Unlike forward methods, reverse methods require access to intermediate values during the gradient

calculation. To avoid saving the entire solution, checkpointing techniques can offer a balance

between storage requirements and re-computation [90].

[Link] Power System Considerations

One unique feature of power system simulations is the inherent inclusion of discrete events

to represent perturbations to the steady state condition. Therefore, prior works on sensitivity

methods for power system dynamics have focused largely on properly handling discrete events in

the forward simulation for both the forward [91] and discrete adjoint [88] cases. In addition to the

trade-offs discussed in the prior sections, the numerical stiffness of the model plays a key role in

how each method scales to large systems (the same is true for the forward simulation) [116, 84].

Stiffness can vary significantly depending on the inclusion of fast electromagnetic dynamics in the

model [11].

In summary, for the numerical sub-problems which make up a power system simulation, there

are complex trade-offs between performance, stability, memory usage, and ease of implementation

for existing methods to compute gradients, any of which could be critical in the context of running

power system dynamic simulations (see [84] Table 1 for a comprehensive comparison). For this

reason, we argue that the best method should depend on the application context. In [88], the

authors focus exclusively on discrete methods for use in optimization problems that require high
 108

accuracy. However, not all gradients are used in an optimization context or require the same level

of accuracy.

5.4 Software Implementation

The gradient calculations described in the previous sections are implemented as an extension

to [Link], an open source power system simulation toolbox [26]. The reader is referred to [11]

for background on the dq0 modeling framework employed in [Link] and to [26] for additional im-

plementation details. Unlike many commercially available tools which offer bespoke combinations

of modeling assumptions and specialized solution methods, [Link] retains full separation of the

modeling and solution layers. This design choice is key for enabling flexibility in changing model-

ing choices and trying novel solution methods with minimal code changes, and is born out of the

realization that there is not a single best solution method for all power system models. Figure 5.2

depicts this separation through a simplified schematic of the software dependencies for executing

a simulation. [Link] handles the power system modeling, including the efficient implementation

of repeated component models within a system-level simulation. The numerical solution of the

component nonlinear and differential equation problems are handled by the generic packages Non-

[Link] and [Link], respectively, which are part of the Julia SciML

ecosystem [25]. The SciML ecosystem in Julia is chosen over others such as PETsC [153, 154] and

Sundials [48] for its combination of performance, scope of available algorithms (both solvers and

sensitivity methods), support of events, and composability. The objective function is implemented

in Julia code from the Base library. We choose to extend [Link] as opposed to using a general

purpose differential programming library (e.g. NeuroMANCER [70]) in order to retain access to

the library of power system specific models and processes for initialization.

The software dependencies for the gradient calculation are shown in Fig. 5.3 (analogous to

Fig. 5.2 for the forward simulation). The high level AD interface for differentiating through a

simulation is provided by Enzyme and the corresponding Julia bindings ([Link]) [114]. En-

zyme offers first-class support for mutation which is used extensively in the base Julia library and
 109
Device [Link]
Parameters

Network
Setpoints Objective
[Link] Function
[Link]
[Link]
Network
Parameters

Figure 5.2: The software dependencies involved in the forward numerical solution of a power system
simulation. The modeling layer consists of [Link] and related packages while the numerical sub-
problems are solved by [Link] and [Link].

within [Link] for efficient parameter and initial condition handling [26]. In addition, unlike many

language-specific AD tools, Enzyme operates on compiler-optimized code, leading to impressive

performance benefits compared to other tools [114]. Finally, Enzyme supports branching code

which is used extensively to implement power system models. While the parameter and initial

condition handling of the forward simulation is handled natively by Enzyme without intervention,

propagating derivatives through the solution of the numerical problems is accomplished via the

Enzyme custom rules interface. The methods for calculating sensitivities of nonlinear equation

solutions and differential equation solutions are included from the [Link] package

and include methods from each of the categorizations discussed in Section 5.3.2.

One of the key practical implementation challenges addressed in this work is the parameter

passing between numerical sub-problems. On one hand, given the complexity of power system sim-

ulations and the need to re-use models, it becomes desirable to organize system-level parameters

hierarchically based on the device. On the other hand, the generic sensitivity methods implemented

in [Link] rely on a flat vector of parameters as an input to the problem. The Com-

[Link] package serves to bridge this gap by simultaneously satisfying the requirements

of the [Link] API while also allowing for a hierarchical structure. The overall
 110
Device [Link]
Parameters

Network
Setpoints
[Link] [Link] [Link]
[Link]
Network
Parameters

Figure 5.3: The software dependencies involved in the gradient calculation of a power system
simulation. [Link] is the top level AD package, and the modularity of the approach comes
from the generic options for differentiating numerical programs from [Link]. The
[Link] API for computing gradients is designed to minimize the amount of modeling code that
must be differentiated by [Link].

parameters are organized first by device and then by type (parameter or setpoint). This allows

for efficient passing of the relevant parameters to each sub-problem (e.g. the power flow solution

or the initialization of a single device model) during the forward pass and correct propagation of

derivatives when calculating the overall gradient.

The code sample of Fig. 5.4 elucidates some of the key features of the API for calculating

gradients. The first line shows the user-defined loss function based on the solution of the simulation

and external data. The second line shows the get sensitivity functions function which returns

the top-level functions for computing the forward simulation and gradient. Note that the the

simulation object (part of the [Link] API) is an input in addition to the parameter(s) of interest

and the system states that are needed for the calculation of the loss. By building on the existing

[Link] API, the forward and gradient functions minimize the amount of code that must be

differentiated with Enzyme; the structures of the simulation are built once and remain largely

unchanged, only the parameters and initial conditions are modified. This is helpful not just for

performance of the gradient calculation but also for making [Link] compatible with [Link]—

a major effort of this work. With the contributions of this paper, gradient calculations for power
 111

Figure 5.4: Code sample demonstrating the computation of the gradient of a simulation.

system simulations in [Link] are made composable and flexible, making it possible to compare

methods for a given system with minimal code changes (e.g. changing the sensealg keyword

argument).

5.5 Results

In this section we showcase a selection of numerical results for various power system models

and applications which rely on the end-to-end differentiable simulations proposed in this work.

We start in Section 5.5.1 with a pedagogical example of the canonical single-machine infinite-bus

system. In Section 5.5.2, we demonstrate how the gradient of a transient simulation can enhance

our understanding of the dynamic behavior of a system. In Sections 5.5.3 and 5.5.4 we apply the

gradient calculation to an optimization context, for tuning physical parameters and training data-

driven models respectively. Finally, in Section 5.5.5 we demonstrate the flexibility of the proposed

architecture by benchmarking the gradient calculation for various sensitivity methods and power

system models. We note that all of the numerical results deal with device level parameters, while

the extension to network parameters (i.e. parameters that change the power flow solution) are left

as future work.
 112

5.5.1 Single Machine Infinite Bus

We start with a simple example of an infinite bus connected to a classical two-state generator

model (without excitation control or governor):

dδ
= Ωb (ω − 1.0) (5.11)
dt
dω 1
= (τ m − τ e − D(ω − 1.0)) (5.12)
dt 2H 0
       
r ′
r −xd   0  r
i  −1  a v 
 =M     − M   (5.13)
ii x′d ra e′q0 vi

where δ and ω are the states, vr and vi are the terminal real and imaginary voltages, and ir and

ii are the real and imaginary injected currents. The initialization of the model involves solving for

δ0 , ω0 , τ0m , e′q0 such that the device is in steady state and the injected current matches the current

from the power flow solution.

Figure 5.5 compares the gradients of the numerical solution of the one machine infinite bus

model with respect to the shaft inertia (H) and the machine d-axis transient reactance (x′d ). The

value of x′d will clearly impact the initial conditions, whereas H, by virtue of appearing in the

denominator of Eq. (5.12), will have no impact. In this case, the objective function is simply

δ(t), making this case analogous to the well known TSA [91]. On the left pane of Fig. 5.5 is

the trajectory of δ in response to a change in the power setpoint at t = 0.5 seconds. In the

middle is the gradient of the simulation with respect to H. The inertia does not impact the initial

conditions, therefore the gradient before the system is perturbed is zero; once the perturbation is

applied the gradient changes. In contrast, the right side shows the gradient for x′d . The non-zero

but constant value before the perturbation is indicative of the parameter’s impact on the initial

condition and is a result of incorporating the initialization into the function that is differentiated.

If the initialization is not included in the gradient computation even though the parameter impacts

the initial conditions, the result is nonsensical; the dashed trajectory is immediately time-varying

even before the perturbation is applied.

 113
Include initialization
Ignore initialization
0.015 1
0.178
0.01 0.8
0.176
0.005
0.6

∂δ/∂xd'
∂δ/∂H
0.174
0
δ

0.4
0.172
−0.005
0.2
0.17 −0.01

0.168 0
−0.015
0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2
Time (s) Time (s) Time (s)

Figure 5.5: Machine angle and gradient with respect to different parameter types for one machine
infinite bus system.

5.5.2 Enhancing Understanding

The gradient of a power system simulation with respect to the parameters can provide in-

formation that is supplemental to the time-domain trajectories. In this section, we demonstrate

how the gradient can enhance the model user’s understanding of the system by encoding which

parameters are most impactful for a specific quantity of interest. For this case study, we consider

a three-bus system with one synchronous generator and one grid-forming inverter.

The left side of Fig. 5.6 (black trace) shows the time-domain trajectory of the voltage

magnitude at the inverter bus in response to a step in the voltage reference; this is the output of

the forward pass of the numerical problem. We formulate the objective function as the maximum

voltage deviation:

fobj = max (|v(t) − v0 |) (5.14)

{t}

The voltage deviation is a key quantity as large voltage deviations can lead to equipment

damage. In the right pane, we show the gradient for this combination of simulation and objective

function with respect to each parameter in the system. The virtual inductance of the inverter

inner control and the reactive power droop gain are the two parameters with the largest magni-

tude gradients; this confirms our basic power system intuition of reactive power being critical for
 114

1.02 Original parameterization Inverter virtual inductance (+)

Inverter virtual inductance (+ 10%) Inverter reactive power droop gain (-)
Voltage magnitude (p.u.)
1.019 0.01
Inverter reactive power droop gain (+ 10%) Other parameters

Magnitude of gradient
1.018
0.001
1.017

1.016 100μ

1.015 10μ
1.014
1μ
1.013

1.012 100n
0 0.02 0.04 0.06 0.08 0.1
Time (s) Parameter (by magnitude of gradient)

Figure 5.6: Time-domain voltage trajectory (left) and gradient for an objective function of the
maximum voltage deviation (right). The magnitudes of the gradients indicate the parameters with
the largest impact on the local voltage.

1.006 Inverse of Governor Droop (+)

Governor Time Constant (-)
0.01
1.005 Other parameters
Magnitude of gradient
Frequency (p.u.)

0.001
1.004
100μ
1.003
Original parameterization 10μ
1.002 Inverse of Governor Droop (+ 10%)
Governor Time Constant (+ 10%) 1μ
1.001
100n
1
0 1 2 3 4 5
Time (s) Parameter (by magnitude of gradient)

Figure 5.7: Time-domain frequency trajectory (left) and gradient for an objective function of the
maximum frequency deviation (right). The magnitudes of the gradients indicate the parameters
with the largest impact on the system frequency.

controlling the local voltage. We run two additional transient simulations where we increase the

highlighted parameters by 10% from the nominal value. For the virtual inductance, which has a

positive gradient, we see that the maximum voltage deviation increases, whereas for the reactive

power droop gain, which has a negative gradient, the effect is opposite. The gradient encodes these
 115
Include initialization Ignore initialization Ground Truth Starting Condition
25 Include initialization Ignore initialization

0.5 0.49
20
0.48
MAE of Xq (%)

0.45
15 0.47
0.4
0.46
0.35

δ1

δ2
10 0.45
0.3 0.44
5 0.25 0.43
0.2 0.42
0
0.15 0.41
10 20 30 40 50 0 5 10 0 5 10
Iteration Time (s) Time (s)

Figure 5.8: Parameter tuning problem where the inclusion of the initialization in the gradient
calculation is critical. The left side shows the loss across iterations. The right side shows the
starting results, along with the best result for each method.

impacts without the need for multiple simulations.

A similar example for the system frequency is included in Fig. 5.7. In this case we consider a

step in the active power reference and formulate the objective as the maximum frequency deviation:

fobj = max (|ω(t) − ω0 |) (5.15)

{t}

The governor droop and a time constant of the governor model have the largest gradients. Again,

the gradient confirms our intuition that the governor is critical to the active power response of the

generator and therefore the frequency.

5.5.3 Parameter Tuning

In this section, we show how our end-to-end method can be critical to obtaining good per-

formance for an optimization problem. We formulate a power system parameter tuning problem as

an optimization problem in which the objective is to minimize a loss function that quantifies the

mismatch between data and predicted trajectories. For ease of demonstration, the ‘ground truth’

data is obtained from the results of a simulation. For a power system operator, the same tuning

could be applied against measured data from the real system. We demonstrate this methodology

on the IEEE 9 bus test system. We tune the q-axis machine reactance by formulating an L2 norm
 116

loss function on the generator angles. We choose this combination of parameters and states because

the q-axis reactance will impact both the initial condition of the generator angle in addition to the

shape of the transient trajectory. This makes it a simple case to demonstrate the value of differ-

entiating through the initialization process. Machine reactances are randomly set for the starting

parameterization between ±40% of the nominal values.

The results of the optimization are shown in Fig. 5.8. On the left side is the mean average

error of the q-axis reactance for the three machines over 50 iterations of gradient descent using

our method for the gradient calculation. The inclusion of the initialization is critical to converging

to the correct trajectory; without it the loss reaches a local minima that does not correspond to

the true parameter values. The right side of the figure shows the trajectories corresponding to the

minimum loss for both methods, in addition to the starting and ground truth trajectories. We note

that while the initialization is not included in the gradient calculation, the model is re-initialized

with the updated parameter value for generating the trajectories in Fig. 5.8, therefore the green

trajectory still starts in steady state.

In this case including the initialization is critical to converging to the correct optimum, how-

ever, there are significant challenges to gradient-based tuning of power systems parameters. For

systems that are near the stability boundary, it is possible to move the system into an unstable

region when updating the parameters. In these cases, a more conservative step size may be nec-

essary. When tuning different types of parameters simultaneously, it may be necessary to scale

the step size given the large range of scales for different power system parameters. Finally, the

problem of parameter identifiability must be addressed depending on the available measurements

[93]. The use case presented here bypasses these challenges by tuning the same single parameter

across multiple devices for a relatively stable system. Strategies for robust tuning of many power

system parameters, either simultaneously or in multiple stages, is an interesting direction for future

research that is enabled by this work.

 117
i(t)
î(t)

v(t) v̂(t)

Figure 5.9: The ground truth 2-area system (left) and the surrogate system with one area replaced
with a neural-ode based surrogate. The predicted and actual trajectories of the interfacing voltage
and current are compared via the loss function and in the following results.

Ground truth Voltage Source Constant Impedance

Neural ODE (untrained) Neural ODE (trained)

Current magnitude (p.u.)

Voltage magnitude (p.u.)

0.994 4.4

0.992 4.2
0.99 4
0.988 3.8
0.986
3.6
0.984
0 2 4 6 8 10 0 2 4 6 8 10
Time (s) Time (s)

Figure 5.10: Comparison of the interfacing voltage and current trajectories for the ground truth
model, data-driven surrogate model, and common power system surrogates.

5.5.4 Training Data-driven Models

In this section we demonstrate the training of a data-driven model within a dynamic power

system simulation. Prior works have begun to explore the benefits of including data-driven models

in dynamic simulations—for example, simulation acceleration [22, 139, 131] and quantifying un-

certainty [77]. Data-driven models require training before use. Prior works have relied on a siloed

approach wherein the models are trained and then subsequently added into a simulation environ-

ment for evaluation. These bespoke solutions require significant time and engineering effort to

implement and represent a barrier to adoption of such models. We address this challenge directly

through the DP paradigm so that the gradients needed for training are available and models can
 118
Voltage Source Constant Impedance Neural ODE

Mean average error

0.1

0.01

0.001

100μ

10μ

v̂ᵣ v̂ᵢ îᵣ îᵢ

Figure 5.11: Mean average error across the test dataset for the data-driven surrogate, a voltage
source surrogate, and a constant impedance surrogate.

be trained and evaluated in the same simulation framework by simply updating the parameters.

We demonstrate on the two-area system [60] by replacing one area with a black-box neural

ode-based surrogate model [22]. The original ground truth and modified surrogate systems are

shown together in Fig. 5.9. We note that the surrogate model in [22] was designed to be com-

patible with existing power system models by having the means to be initialized to a steady state

operating condition that matches the power flow. For this reason, the gradient which incorporates

the initialization is required to train such a surrogate within a simulation.

We train the data-driven model by formulating the objective (loss) function as the RMS

error between the actual and predicted voltages and currents at the interfacing buses (the labeled

quantities in Fig. 5.9). The training dataset consists of five simulations with an active power step

of varying magnitude. The gradients are computed using a reverse and discrete adjoint method and

the parameters are iteratively updated according to the Adam optimization algorithm [126]. Figure

5.10 compares the response of the data-driven surrogate system with the ground truth response

both before and after training. As an additional point of comparison, the response of the system

if the area is replaced with an ideal voltage source (infinite bus) or constant impedance is also
 119

included. Figure 5.11 shows the aggregated errors over a test dataset consisting of 50 random

active power perturbations within the range defined by the train dataset, demonstrating that the

data-driven model can dramatically outperform the voltage source and constant impedance on each

of the interfacing quantities.

While this example shows a single surrogate, a realistic use case of data-driven models for

a large system is likely to involve training and deploying many surrogate models. In this case,

training the models within the simulation workflow can dramatically simplify the workflow; it is

possible to train multiple models simultaneously by simply specifying the collection of parameters

that should be updated according to the computed gradient. Training data-driven models such as

the neural ode-based surrogate over larger and more varied datasets without getting stuck in local

minima is a key direction for future work.

5.5.5 Benchmarking

In this section we demonstrate the flexibility of the framework by benchmarking the forward

and gradient computation across various power system models, fault types, options for differ-

entiating through the differential equation solve, the differential equation solver, and tolerances.

We include results from four systems of various sizes, along with methods from each category of

approaches discussed in Section 5.3. Note, however, that only one forward method is included.

In practice, the forward continuous and discrete methods give the same result, as the continuous

method is implemented by solving the augmented equations simultaneously with the original system

[87]. The discrete method is preferred due to the ease of handling system events. Each simulation

is run for 5 s with a 0.1 p.u. step in the active power setpoint of a generator applied at 1 s.

The forward and gradient execution times for this benchmarking exercise are plotted in Fig.

5.12. Each dot represents a unique combination of the one of the inputs described above. For

simplicity, the data is color-coded according to the option for differentiating the solution of the

ODEs. Although slightly obscured by the variety of options represented in the data, we can see

that the forward method is advantageous for easier problems with fewer parameters, but does
 120
Forward Reverse (Continuous) Reverse (Discrete)

Gradient Execution Time (s)

100

10

0.001 0.01 0.1 1 10

Forward Execution Time (s)

Figure 5.12: Execution time benchmarking. Each point represents a distinct combination of system
model, fault type, sensitivity algorithm for differentiating the differential equation solve, differential
equation solver, and tolerance.

not scale favorably. This result demonstrates the flexibility of the methodology; for a specific

application with specifications for computation time and/or accuracy, the ability to benchmark

across many methods is key to choosing the best solver/sensitivity options. Open benchmarks for

adjoint techniques are common in other fields (e.g. [140]), and are an important mechanism for

identifying the most promising techniques and implementations.

5.6 Conclusion

This paper proposes a method for differentiable power system dynamic simulations. Unlike

prior works which focus on the solution of the DEs, our end-to-end method also includes the steps

to initialize the system model to a feasible operating condition. In addition, the proposed method

is flexible and modular, based on existing modeling and differentiation tools, allowing for efficient

benchmarking of the wide array of techniques for differentiating numerical problems.

As part of an open source tool, there will be many avenues for extensions and future im-

provements. One key ongoing area of future work is the extension to parameters which impact the
 121

power flow solution in addition to the device level initial conditions. Another key improvement is

to simplify the API for computing gradients and take further advantage of existing functionality

in [Link]. Future efforts will also focus on improving the performance through novel adjoint

implementations, including the direct application of Enzyme reverse mode AD to a power system

simulation. Such improvements and novel methods will be necessary to scale to larger and more

stiff power system models. In addition, there are many potential avenues of future work which use

gradients for downstream problems, including the inclusion of data-driven models within dynamic

simulations.
 Chapter 6

Conclusions and Future Work

The high-level motivation for this thesis is to mitigate the negative impacts of anthropogenic

climate change. Decarbonization of the electric power system in conjunction with electrification of

other fossil-fuel based industries is one popular strategy to economically and efficiently decarbonize

society. Luckily, the constituent technologies necessary to make this transition—wind turbines, solar

panels, electric vehicles, among others—exist and are already in many cases the most economical

options [4]. The transition is underway and gaining momentum. Given the fundamental nature of

this transition, many of the tools and processes for studying and operating the power system must

also evolve. This thesis focuses on transient numerical simulations, one of the key tools for studying

the large-signal dynamic behavior of the power system. In this area, we need improved tools that

are both high fidelity and computationally tractable to successfully navigate the energy transition.

Incorporating data-driven or ML models into transient simulations is one promising approach to

improve such tools.

Despite the promise of ML for modeling dynamical systems, there are significant barriers to

their adoption. First, there are incongruencies between ML models that have been proposed in

the literature, often in other fields, and the typical workflow of a power system simulation. Two

of the major contributions of this thesis, the introduction of a neural ODE surrogate that can

be appropriately initialized and the integration of direct solution surrogates via DDEs, address

these incongruencies and make data-driven models coherent with existing practices. Second, the

siloed process for training and deploying data-driven models requires extensive bespoke effort and
 123

is not operationally feasible. We address this limitation via the development of an end-to-end

differentiable simulation tool that enables online training of data-driven models.

6.1 Summary of Contributions

In Chapter 3, I propose a novel data-driven surrogate model consisting of a DEQ layer and

a neural ODE to accelerate time-domain simulations via a reduction in the number of states. In a

typical power system simulation workflow, the power flow solution determines the initial states of

the network and individual devices are initialized to meet the network interface requirements (i.e.

matching the power flow and being in steady state). Prior applications of neural ODEs to power

systems have not considered the initialization requirements. We address this shortcoming with a

DEQ layer which solves the nonlinear equations associated with the neural ODE and provides the

initial conditions which satisfy the network interface requirements. For a medium size system with

both synchronous generators and inverters, we demonstrate an average of 5x simulation acceleration.

Our black-box model achieves accuracy comparable to a physics-based surrogate model without

relying on full knowledge of the underlying system structure and parameters which might not be

available in practice.

In Chapter 4, I propose a methodology for cohesively combining two seemingly distinct types

of dynamic models; data-driven models which predict directly the solution of a system at a future

point in time and differential equation-based models that are used in conjunction with a numerical

solver to advance the solution in time. These types of models have been used together in the

past by severely restricting the solution methods available for the resulting system. In some cases,

the methodologies impose a fixed-step solution method and in others require that the data-driven

models are implemented as part of the solution method. We propose a method to integrate such

models by adding delay equations and thereby converting the original system of DEs to a modified

system of DDEs. In this way, we retain full separation of the novel data-driven models and the

solution techniques, which is a key design principle of this thesis. The composability of the method

of steps allows us to retain access to existing ODE solvers including high-performance variable
 124

time step methods. Finally, we comprehensively characterize the numerical tradeoffs between the

ODE and DDE systems and explore and demonstrate various mechanisms for acceleration, such

as reducing the order of the state space and increasing the prediction horizon of the data-driven

model.

In Chapter 5, I propose a methodology for end-to-end differentiable power system simula-

tions. In the context of data-driven models, this work dramatically simplifies the process of model

training. Prior works have relied on a siloed approach of first training data-driven models in iso-

lation and then re-integrating them into a simulation for testing. In our approach, data-driven

models with random initial parameters are added to the system model akin to traditional models

and the parameters are iteratively updated according to the gradient. Models are automatically

trained in the same environment as they are deployed (i.e. within a simulation) without additional

engineering effort or techniques. This contribution also allows for flexibility during the training

process; multiple models can be trained simultaneously, physics-based parameters can be trained

in conjunction with data-driven models, or models can be re-trained online as new data becomes

available. Overall, this improvement lowers the barrier to using data-driven models in industrial

practice by reducing the specialized effort needed for training. The scope of Chapter 5 also ex-

tends beyond data-driven models. Prior applications have used gradients of dynamic simulations

for parameter tuning, uncertainty quantification, and incorporating dynamic constraints into opti-

mization problems. Given the wide range of potential applications, we employ composable software

packages to maintain flexibility over the myriad numerical techniques for solving and differentiating

through numerical problems. This allows for choosing the appropriate methods for a specific appli-

cation, whether it be enhancing our understanding of a illustrative test system or simultaneously

updating thousands of parameters of a deep neural network surrogate model.

6.2 Limitations and Future Work

This thesis explores using data-driven models in conjunction with traditional first-principles

models to enhance our ability to accurately and efficiently simulate and understand power systems.
 125

While the results show promise, the bar for adopting such models operationally is high. The power

system industry is reluctant to change from existing approaches which have been validated over

many years. This thesis alone does not comprehensively make the case for adopting data-driven

models. In this section we highlight some of the primary limitations of the work. Each item suggests

a promising direction for future work which is detailed in the subsequent sections, with the hope

of inspiring researchers to continue this line of work.

• Training methods for robustly training data-driven models across large and diverse power

system datasets are underexplored in this thesis. In Chapter 3, significant time and hyper-

parameter tuning was required to achieve satisfactory performance. A more fundamental

look at the training methods and how they can be adapted to the unique challenges of

power systems datasets is proposed in Section 6.2.1.

• This thesis focuses on acceleration of simulations while maintaining acceptable simulation

accuracy compared to a ground truth system. While this is one promising application,

there are other potential benefits of using data-driven models that can help strengthen the

case for adoption. Exploring data-driven models for uncertainty quantification, especially

in representing distribution systems, is proposed in Section 6.2.2.

• This thesis proposes techniques to incorporate data-driven models and improve simulations,

but does not explicitly consider downstream applications. In Section 6.2.3, I suggest con-

tingency analysis as a promising downstream application of simulations with data-driven

models. Establishing common benchmark problems has led to rapid improvement of ML

in other areas [155, 156] and should be repeated for this problem in order to gain a more

fair and complete comparison between approaches [157].

• This thesis combines black-box data-driven device models and traditional device models

at the system level. This approach allows for flexibility in deploying data-driven models

according to engineering judgement. Future work should continue to explore a more tightly
 126

coupled approach wherein physics-based constraints or knowledge are directly embedded

into the data-driven models. In Section 6.2.4 I propose the adoption of PINNs as an

automated modeling option for [Link].

6.2.1 Robust Training of Data-driven Models for Power System Dynamic Simu-

lations

One of the general challenges of training ML models via iterative gradient-based algorithms

is avoiding local minima in the loss landscape. Random batching of the training dataset, decreasing

the batch size (stochastic gradient descent), or increasing the stepsize are all common approaches

to alleviate this problem. Unfortunately, there is no ideal value for each of these hyper-parameters,

and therefore they must be tuned for a specific learning problem to achieve optimal performance.

While optimization algorithms such as stochastic gradient descent and Adam are most prevalent,

there are also specialized algorithms such as simulated annealing [158] which are designed with the

avoidance of local minima in mind.

Training models which are a component of dynamic simulations (i.e. the training problems

described in this thesis) presents unique challenges with respect to local minima. Deep learning

models are prone to learn low frequency behavior first [159]; in dynamic simulations, this manifests

as models learning the average behavior of a dynamic trajectory but failing to capture the higher

frequency details. This behavior can be especially problematic for power systems applications

where stability issues often manifest as sustained oscillations. Techniques such as extending the

timespan of the simulation throughout training or breaking the training problem into multiple

conjoined sections [160] have been proposed to deal with the local minima issue for dynamic data-

driven models, but have not been studied with respect to power system problems. Another unique

challenge in power systems is capturing the distinct behavior of the same dynamic model depending

across different perturbations to the system. Training multiple models for each type of fault or

encoding the fault as an input to the model can potentially alleviate this modeling burden and

improve the training performance.

 127

Future work should focus on comparing the techniques above and/or proposing novel methods

for power systems applications. We also note that the difficulty of training is closely linked to the

form of the data-driven model. It is likely that by incorporating more physics-based knowledge

into the data-driven models we can help alleviate the training burden, at the expense of requiring a

more bespoke model and additional knowledge of the true system (see Section 6.2.4). For example,

learning a small correction term within the control of an inverter will be easier from a training

perspective than learning the dynamic response of the entire device via a black-box model.

6.2.2 Capturing Uncertainty in Power System Transient Models

Model aggregation is required in some form to make simulations of very large power systems

computationally tractable. Aggregation is especially common at the interface of the distribution

and transmission systems. A typical approach is to derive a deterministic aggregate model that

approximates a high-fidelity model (e.g. [66]), but there is no quantifiable guarantee as to how

well the aggregate model will capture the behavior of the high-fidelity model in future simulations.

In this context, data-driven models that can numerically quantify the uncertainty associated with

a prediction [161] have the potential to give a more complete and useful representation. This

additional information can provide safeguards when using aggregate models, by capturing instances

where additional training or a more detailed simulation is necessary.

In addition to capturing the modeling uncertainty, data-driven models can also provide a

mechanism for capturing the inherent uncertainty of the system. The energy transition introduces

new forms of uncertainty; wind and solar generation are subject to uncertain forecast errors and

with more controllable devices at the grid edge, the dynamics of distribution systems are increas-

ingly dependent on complex and uncertain human behavior. The assumption that we can safely

ignore these uncertainties while evaluating the stability of power systems at fast time scales is not

guaranteed to hold. Deep learning models such as Bayesian neural networks can capture and quan-

tify both aleatoric and epistemic uncertainty [162] and have begun to be applied to power systems

simulation applications [163]. Given the computational limitations that already exist on transient
 128

simulations and the additional computational overhead of training Bayesian deep-learning models,

an important task in this area will be to understand the conditions when uncertainty quantification

is critical and the conditions when it can be safely ignored.

6.2.3 Benchmarking Data-driven Models for Contingency Analysis

Throughout this thesis, we quantify the performance of data-driven models by comparing the

accuracy relative to a ground truth model and the time required for simulation. Transient simula-

tions do not occur in isolation, but rather are ultimately used to better understand the system or

inform decision making. One promising application for data-driven models is contingency analysis,

the continual screening of feasible contingencies in close to real time as the system operating condi-

tions evolve. This process is key to identifying potential stability issues that could occur in response

to a credible contingency and taking remedial actions to limit the impact. Contingency analysis

is a promising application for the methodologies of this thesis for two key reasons. First, there

is a straightforward benefit for simulation acceleration; in a time constrained application faster

simulations enables more simulations and therefore a more robust evaluation of the contingency

space. Second, contingency analysis involves repeated simulations of a system over a long period of

time. This means that the initial cost of training data-driven models can be amortized over many

simulations making the tradeoff more attractive.

The first step in this line of work is to clearly quantify a set of requirements via a benchmark

problem. A clear benchmark which is developed independently of solutions enables fair comparisons

between data-driven models in addition to solutions outside of the scope of ML [164]. Ideally, this

process would include input from stakeholders with knowledge of current operational practices for

contingency analysis. Consistent benchmarking has been a hallmark of the ML community and

is a significant factor in the rapid success and improvement of model architectures. Developing a

comprehensive benchmark for a contingency analysis problem will be challenging as there are many

hard to quantify criteria (e.g. ease of adoption). The methods for incorporating data-driven models

developed on top of [Link] as part of this thesis make it an ideal environment for benchmark
 129

problems in this space.

6.2.4 Physics Informed Neural Networks for Transient Simulations

Even if researchers can demonstrate robust training methods for data-driven models (Section

6.2.1), acceleration of simulations when using such models (Chapter 3), and additional features such

as uncertainty quantification (Section 6.2.2), the issue of trusting data-driven models to capture

the behavior of the system generally and not only for select examples remains. To address this

concern requires extensive benchmarking, including on challenging test problems that are near the

stability boundary and involve multiple data-driven models interacting.

One potential path to address these concerns is to systematically incorporate physics informed

neural networks into [Link] such that every device can be interchangeably modeled with DEs

or PINNs. This setup would allow for the extensive and robust demonstration of data-driven

models necessary for adoption. PINNs provide additional benefits for this use case compared to

the black-box models which are the focus of this thesis. First, the state space of PINNs matches

the differential equation-based model; this is an important feature for building trust in the model

as the individual states retain their physical interpretability regardless of which components are

modeled with neural networks. Second, the physics-based regularization which defines PINNs can

improve training performance, either through improved generalization [78] or by requiring less data

[165]. The training processes within [Link] can be highly automated with minimal additional

modeling effort as the equations required for regularization are already implemented as part of

the differential equation-based models. The work of this thesis has already addressed many of the

obstacles to enabling this proposal: PINNs are an example of a direct solution surrogate so the

methods of Chapter 4 are directly applicable to building simulations with PINNs and DEs. In

addition, Chapter 5 enables the online training of PINNs within [Link].

 Bibliography

[1] “Paris Agreement to the United Nations Framework Convention on Climate Change,” Dec.
2015. T.I.A.S. No. 16-1104.
[2] V. Masson-Delmotte, P. Zhai, H.-O. Pörtner, D. Roberts, J. Skea, P. Shukla, A. Pirani,
W. Moufouma-Okia, C. Péan, R. Pidcock, S. Connors, J. Matthews, Y. Chen, X. Zhou,
M. Gomis, E. Lonnoy, T. Maycock, M. Tignor, and T. Waterfield, “Global Warming of 1.5
ºC,” special report, The Intergovernmental Panel on Climate Change, 2018.
[3] “Net Zero by 2050 - A Roadmap for the Global Energy Sector,” tech. rep., International
Energy Agency, Oct. 2021.
[4] M. Roser, “Why did renewables become so cheap so fast?,” Our World in Data, 2020.
[5] “Renewable power generation costs in 2022: Executive summary,” tech. rep., International
Renewable Energy Agency, Abu Dhabi, 2023.
[6] Y. Lin, J. H. Eto, B. B. Johnson, J. D. Flicker, R. H. Lasseter, H. N. V. Pico, G.-S. Seo, B. J.
Pierre, and A. Ellis, “Research Roadmap on Grid-Forming Inverters,” Tech. Rep. NREL/TP-
5D00-73476, National Renewable Energy Laboratory, 2020.
[7] N. Hatziargyriou, J. Milanovic, C. Rahmann, V. Ajjarapu, C. Canizares, I. Erlich, D. Hill,
I. Hiskens, I. Kamwa, B. Pal, P. Pourbeik, J. Sanchez-Gasca, A. Stankovic, T. Van Cutsem,
V. Vittal, and C. Vournas, “Definition and Classification of Power System Stability – Re-
visited and Extended,” IEEE Transactions on Power Systems, vol. 36, pp. 3271–3281, July
2021.
[8] J. Cochran, M. Miller, O. Zinaman, M. Milligan, D. Arent, B. Palmintier, M. O’Malley,
S. Mueller, E. Lannoye, A. Tuohy, B. Kujala, M. Sommer, H. Holttinen, J. Kiviluoma, and
S. K. Soonee, “Flexibility in 21st Century Power Systems,” Tech. Rep. NREL/TP-6A20-
61721, May 2014.
[9] “Forecasting increasingly critical to harnessing wind and solar for power systems,” AEMO,
Nov. 2023.
[10] “Beyond Positive Sequence RMS Simulations for High DER Penetration Conditions,” tech.
rep., North American Electric Reliability Corporation, Oct. 2022.
[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,” IEEE
Transactions on Power Systems, pp. 1–16, 2023.
 131

[12] U. Markovic, J. Vorwerk, P. Aristidou, and G. Hug, “Stability Analysis of Converter Control
Modes in Low-Inertia Power Systems,” in 2018 IEEE PES Innovative Smart Grid Technologies
Conference Europe (ISGT-Europe), (Sarajevo, Bosnia and Herzegovina), pp. 1–6, Oct. 2018.

[13] R. Henriquez-Auba, J. D. Lara, C. Roberts, and D. S. Callaway, “Grid Forming Inverter

Small Signal Stability: Examining Role of Line and Voltage Dynamics,” in IECON 2020 The
46th Annual Conference of the IEEE Industrial Electronics Society, (Singapore, Singapore),
pp. 4063–4068, IEEE, Oct. 2020.

[14] J.-S. Brouillon, M. Colombino, D. Groß, and F. Dörfler, “The effect of transmission-line dy-
namics on a globally synchronizing controller for power inverters,” in 2018 European Control
Conference (ECC), pp. 2242–2247, June 2018.

[15] W. Wang, G. M. Huang, P. Kansal, L. E. Anderson, R. J. O’Keefe, D. Ramasubramanian,

P. Mitra, and E. Farantatos, “Instability of PLL-Synchronized Converter-Based Generators
in Low Short-Circuit Systems and the Limitations of Positive Sequence Modeling,” in 2018
North American Power Symposium (NAPS), pp. 1–6, Sept. 2018.

[16] B. Badrzadeh, “Electromagnetic transient simulation models for large-scale system impact
studies in power systems having a high penetration of inverter-based resources,” Energy
Systems Integration Group, June 2019.

[17] “Odessa Disturbance Texas Events: May 9, 2021 and June 26, 2021 Joint NERC and Texas
RE Staff Report,” tech. rep., Sept. 2021.

[18] R. W. Kenyon, B. Wang, A. Hoke, J. Tan, C. Antonio, and B.-M. Hodge, “Validation of
Maui PSCAD Model: Motivation, Methodology, and Lessons Learned,” in 2020 52nd North
American Power Symposium (NAPS), (Tempe, AZ, USA), pp. 1–6, IEEE, Apr. 2021.

[19] E. Rehman, M. Miller, J. Schmall, and S. H. Huang, “Dynamic Stability Assessment of High
Penetration of Renewable Generation in the ERCOT Grid,” Tech. Rep. Version 1.0, Apr.
2019.

[20] R. de Silva, “Review of AEMO’s PSCAD Modelling of the Power System in South Australia,”
Tech. Rep. JA6921 – Final Report, Dec. 2017.

[21] G. Gurrala, D. L. Dinesha, A. Dimitrovski, P. Sreekanth, S. Simunovic, and M. Starke, “Large

Multi-Machine Power System Simulations Using Multi-Stage Adomian Decomposition,” IEEE
Transactions on Power Systems, vol. 32, pp. 3594–3606, Sept. 2017.

[22] M. Bossart, J. D. Lara, C. Roberts, R. Henriquez-Auba, D. Callaway, and B.-M. Hodge,

“Acceleration of Power System Dynamic Simulations using a Deep Equilibrium Layer and
Neural ODE Surrogate,” May 2024. arXiv:2405.06827 [[Link]], preprint.

[23] R. Galarza, J. Chow, W. Price, A. Hargrave, and P. Hirsch, “Aggregation of exciter models
for constructing power system dynamic equivalents,” IEEE Transactions on Power Systems,
vol. 13, pp. 782–788, Aug. 1998.

[24] Y. Liu and K. Sun, “Solving Power System Differential Algebraic Equations Using Differential
Transformation,” IEEE Transactions on Power Systems, vol. 35, pp. 2289–2299, May 2020.
 132

[25] C. Rackauckas and Q. Nie, “[Link] – A Performant and Feature-Rich Ecosys-

tem for Solving Differential Equations in Julia,” Journal of Open Research Software, vol. 5,
May 2017.
[26] J. D. Lara, R. Henriquez-Auba, M. Bossart, D. S. Callaway, and C. Barrows, “PowerSim-
[Link] – An Open Source Modeling Package for Modern Power Systems with
Inverter-Based Resources,” Aug. 2023. arXiv:2308.02921 [[Link]], preprint.
[27] M. A. Carreño, D. F. Rodrı́guez M, and J. M. Gers, “A Simple Application for Executing
Dynamic Simulations in Parallel on Cluster Machines,” in 2022 IEEE Power & Energy Society
Innovative Smart Grid Technologies Conference (ISGT), pp. 1–4, Apr. 2022.
[28] M. Xiong, B. Wang, D. Vaidhynathan, J. Maack, M. J. Reynolds, A. Hoke, K. Sun, and J. Tan,
“ParaEMT: An Open Source, Parallelizable, and HPC-Compatible EMT Simulator for Large-
Scale IBR-Rich Power Grids,” IEEE Transactions on Power Delivery, vol. 39, pp. 911–921,
Apr. 2024.
[29] Y. Song, Y. Chen, S. Huang, Y. Xu, Z. Yu, and W. Xue, “Efficient GPU-Based Electromag-
netic Transient Simulation for Power Systems With Thread-Oriented Transformation and
Automatic Code Generation,” IEEE Access, vol. 6, pp. 25724–25736, 2018.
[30] W. Chen, J. Xu, K. Wang, P. Wu, J. Ding, and D. Fan, “Steady State Power System Simu-
lation Based on GPU-Accelerated Holomorphic Embedding Method,” in 2022 Power System
and Green Energy Conference (PSGEC), pp. 890–894, Aug. 2022.
[31] S. Subedi, M. Rauniyar, S. Ishaq, T. M. Hansen, R. Tonkoski, M. Shirazi, R. Wies, and
P. Cicilio, “Review of Methods to Accelerate Electromagnetic Transient Simulation of Power
Systems,” IEEE Access, vol. 9, pp. 89714–89731, 2021.
[32] J. Rand, N. Manderlink, W. Gorman, R. Wiser, J. Seel, J. Kemp, S. Jeong, and F. Kahrl,
“Queued Up: 2024 Edition, Characteristics of Power Plants Seeking Transmission Intercon-
nection As of the End of 2023,” tech. rep., Lawrence Berkeley National Laboratory, Apr.
2024.
[33] G. Vassell, “Northeast Blackout of 1965,” IEEE Power Engineering Review, vol. 11, pp. 4–,
Jan. 1991.
[34] P. Kundur, J. Paserba, V. Ajjarapu, G. Andersson, A. Bose, C. Canizares, N. Hatziargyriou,
D. Hill, A. Stankovic, C. Taylor, T. Van Cutsem, and V. Vittal, “Definition and classification
of power system stability IEEE/CIGRE joint task force on stability terms and definitions,”
IEEE Transactions on Power Systems, vol. 19, pp. 1387–1401, Aug. 2004.
[35] “Interconnection Oscillation Analysis: Reliability Assessment,” tech. rep., North American
Electric Reliability Corporation, July 2019.
[36] R. W. Kenyon, M. Bossart, M. Marković, K. Doubleday, R. Matsuda-Dunn, S. Mitova, S. A.
Julien, E. T. Hale, and B.-M. Hodge, “Stability and control of power systems with high
penetrations of inverter-based resources: An accessible review of current knowledge and open
questions,” Solar Energy, vol. 210, pp. 149–168, Nov. 2020.
[37] G. H. Golub and J. M. Ortega, Scientific Computing and Differential Equations: An
Introduction to Numerical Methods. Academic Press, 1992.
 133

[38] E. Babulak and M. Wang, “Discrete Event Simulation,” IntechOpen, 2010.

[39] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning. Springer
Series in Statistics, Springer New York, 2009.

[40] K. Hornik, “Approximation capabilities of multilayer feedforward networks,” Neural

Networks, vol. 4, pp. 251–257, Jan. 1991.

[41] D. Kiela, M. Bartolo, Y. Nie, D. Kaushik, A. Geiger, Z. Wu, B. Vidgen, G. Prasad, A. Singh,
P. Ringshia, Z. Ma, T. Thrush, S. Riedel, Z. Waseem, P. Stenetorp, R. Jia, M. Bansal,
C. Potts, and A. Williams, “Dynabench: Rethinking Benchmarking in NLP,” in Proceedings
of the 2021 Conference of the North American Chapter of the Association for Computational
Linguistics: Human Language Technologies, (Online), pp. 4110–4124, 2021.

[42] N. Baker, F. Alexander, T. Bremer, A. Hagberg, Y. Kevrekidis, H. Najm, M. Parashar,

A. Patra, J. Sethian, S. Wild, K. Willcox, and S. Lee, “Workshop Report on Basic Research
Needs for Scientific Machine Learning: Core Technologies for Artificial Intelligence,” Tech.
Rep. 1478744, Feb. 2019.

[43] J. K. Snelson, “Propagation of Travelling Waves on Transmission Lines - Frequency De-

pendent Parameters,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-91,
pp. 85–91, Jan. 1972.

[44] J. Schoene and T. E. McDermott, “Modeling large-scale power systems with frequency-
dependent transmission line models,” in IEEE/PES Transmission and Distribution
Conference and Exposition, pp. 1–7, Apr. 2008.

[45] H. W. Dommel, “Digital Computer Solution of Electromagnetic Transients in Single-and

Multiphase Networks,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-88,
pp. 388–399, Apr. 1969.

[46] “Brief Summary of all 2nd Generation Generic Renewable Energy System Dynamic Models,”
tech. rep., Modeling and Validation Subcommittee Renewable Energy Modeling Working
Group, WECC, Feb. 2021.

[47] “Technical Report: BPS-Connected Inverter-Based Resource Modeling and Studies,” tech.
rep., North American Electric Reliability Corporation, May 2020.

[48] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and

C. S. Woodward, “SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers,”
ACM Transactions on Mathematical Software, vol. 31, pp. 363–396, Sept. 2005.

[49] A. Bojańczyk, “Complexity of Solving Linear Systems in Different Models of Computation,”

SIAM Journal on Numerical Analysis, vol. 21, no. 3, pp. 591–603, 1984.

[50] M. Senyuk, M. Safaraliev, F. Kamalov, and H. Sulieman, “Power System Transient Stabil-
ity Assessment Based on Machine Learning Algorithms and Grid Topology,” Mathematics,
vol. 11, p. 525, Jan. 2023.

[51] C. Ren, Y. Xu, and R. Zhang, “An Interpretable Deep Learning Method for Power Sys-
tem Transient Stability Assessment via Tree Regularization,” IEEE Transactions on Power
Systems, vol. 37, pp. 3359–3369, Sept. 2022.
 134

[52] T. Meridji, G. Joós, and J. Restrepo, “A power system stability assessment framework using
machine-learning,” Electric Power Systems Research, vol. 216, p. 108981, Mar. 2023.

[53] P. Sarajcev, A. Kunac, G. Petrovic, and M. Despalatovic, “Power System Transient Stability
Assessment Using Stacked Autoencoder and Voting Ensemble,” Energies, vol. 14, p. 3148,
Jan. 2021.

[54] S. Kunkolienkar, F. Safdarian, J. Snodgrass, A. Birchfield, and T. Overbye, “A Description

of the Texas A&M University Electric Grid Test Case Repository for Power System Studies,”
in 2024 IEEE Texas Power and Energy Conference (TPEC), pp. 1–6, Feb. 2024.

[55] P. Demetriou, M. Asprou, J. Quiros-Tortos, and E. Kyriakides, “Dynamic IEEE Test Systems
for Transient Analysis,” IEEE Systems Journal, vol. 11, pp. 2108–2117, Dec. 2017.

[56] T. Xu, A. B. Birchfield, and T. J. Overbye, “Modeling, Tuning, and Validating System Dy-
namics in Synthetic Electric Grids,” IEEE Transactions on Power Systems, vol. 33, pp. 6501–
6509, Nov. 2018.

[57] K. S. Shetye, B. Yancey, and T. J. Overbye, “Benchmarking Standard Library and User De-
fined Models of Renewable Generation,” in 2021 North American Power Symposium (NAPS),
pp. 1–6, Nov. 2021.

[58] A. Büttner, A. Plietzsch, M. Anvari, and F. Hellmann, “A framework for synthetic power sys-
tem dynamics,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 33, p. 083120,
Aug. 2023.

[59] R. Kogler, A. Plietzsch, P. Schultz, and F. Hellmann, “Normal Form for Grid-Forming Power
Grid Actors,” PRX Energy, vol. 1, p. 013008, June 2022.

[60] P. Kundur, Power System Stability And Control. New York: McGraw- Hill, 1994.

[61] V. Purba, B. B. Johnson, M. Rodriguez, S. Jafarpour, F. Bullo, and S. V. Dhople, “Reduced-

order Aggregate Model for Parallel-connected Single-phase Inverters,” IEEE Transactions on
Energy Conversion, vol. 34, pp. 824–837, June 2019.

[62] P. Pourbeik, J. J. Sanchez-Gasca, J. Senthil, J. D. Weber, P. S. Zadehkhost, Y. Kazachkov,

S. Tacke, J. Wen, and A. Ellis, “Generic Dynamic Models for Modeling Wind Power Plants
and Other Renewable Technologies in Large-Scale Power System Studies,” IEEE Transactions
on Energy Conversion, vol. 32, pp. 1108–1116, Sept. 2017. Conference Name: IEEE Trans-
actions on Energy Conversion.

[63] S. Mat Zali and J. V. Milanović, “Generic Model of Active Distribution Network for Large
Power System Stability Studies,” IEEE Transactions on Power Systems, vol. 28, pp. 3126–
3133, Aug. 2013. Conference Name: IEEE Transactions on Power Systems.

[64] “Model User Guide for Generic Renewable Energy System Models,” tech. rep., Electric Power
Research Institute, July 2018.

[65] V. Purba, B. B. Johnson, S. Jafarpour, F. Bullo, and S. V. Dhople, “Dynamic Aggregation of

Grid-Tied Three-Phase Inverters,” IEEE Transactions on Power Systems, vol. 35, pp. 1520–
1530, Mar. 2020.
 135

[66] Z. Ma, Z. Wang, Y. Wang, R. Diao, and D. Shi, “Mathematical Representation of WECC
Composite Load Model,” Journal of Modern Power Systems and Clean Energy, vol. 8,
pp. 1015–1023, Sept. 2020.

[67] C. Legaard, T. Schranz, G. Schweiger, J. Drgoňa, B. Falay, C. Gomes, A. Iosifidis, M. Abkar,

and P. Larsen, “Constructing Neural Network Based Models for Simulating Dynamical Sys-
tems,” ACM Computing Surveys, vol. 55, pp. 1–34, Nov. 2023.

[68] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural Ordinary Differential

Equations,” in Advances in Neural Information Processing Systems 31 (NeurIPS), 2018.

[69] S. Bai, J. Z. Kolter, and V. Koltun, “Deep Equilibrium Models,” in Advances in Neural
Information Processing Systems 32 (NeurIPS), 2019.

[70] J. Drgona, A. Tuor, J. Koch, M. Shapiro, J. Bruno, and V. Draguna, “NeuroMANCER:

Neural Modules with Adaptive Nonlinear Constraints and Efficient Regularizations,” 2023.

[71] C. Rackauckas, M. Innes, Y. Ma, J. Bettencourt, L. White, and V. Dixit, “[Link] - A

Julia Library for Neural Differential Equations,” Feb. 2019.

[72] T. Xiao, Y. Chen, S. Huang, T. He, and H. Guan, “Feasibility Study of Neural ODE and
DAE Modules for Power System Dynamic Component Modeling,” IEEE Transactions on
Power Systems, pp. 1–13, 2022.

[73] Y. Zhou and P. Zhang, “Neuro-Reachability of Networked Microgrids,” IEEE Transactions

on Power Systems, vol. 37, pp. 142–152, Jan. 2022.

[74] A. Azmy, I. Erlich, and P. Sowa, “Artificial neural network-based dynamic equivalents for
distribution systems containing active sources,” IEE Proceedings - Generation, Transmission
and Distribution, vol. 151, no. 6, p. 681, 2004.

[75] C. Han, D. Changhong, and L. Dalu, “Recurrent Neural Network-based Dynamic Equivalenc-
ing in Power System,” in 2007 IEEE International Conference on Control and Automation,
pp. 2396–2399, May 2007.

[76] F. Calero, C. A. Cañizares, and K. Bhattacharya, “Aggregated BESS Dynamic Models for
Active Distribution Network Studies,” IEEE Transactions on Smart Grid, vol. 12, pp. 2077–
2088, May 2021.

[77] J. Vorwerk, T. Zufferey, P. Aristidou, and G. Hug, “Using Quantile Forecasts for Dynamic
Equivalents of Active Distribution Grids under Uncertainty,” in 11th Bulk Power Systems
Dynamics and Control Symposium (IREP 2022), (Banff, Canada), p. 13, July 2022.

[78] M. Raissi, P. Perdikaris, and G. Karniadakis, “Physics-informed neural networks: A deep

learning framework for solving forward and inverse problems involving nonlinear partial dif-
ferential equations,” Journal of Computational Physics, vol. 378, pp. 686–707, Feb. 2019.

[79] G. S. Misyris, A. Venzke, and S. Chatzivasileiadis, “Physics-Informed Neural Networks for

Power Systems,” in 2020 IEEE Power Energy Society General Meeting (PESGM), pp. 1–5,
Aug. 2020.
 136

[80] J. Stiasny, S. Chevalier, and S. Chatzivasileiadis, “Learning without Data: Physics-

Informed Neural Networks for Fast Time-Domain Simulation,” in 2021 IEEE International
Conference on Communications, Control, and Computing Technologies for Smart Grids
(SmartGridComm), pp. 438–443, Oct. 2021.

[81] Y. N. Raffoul, “Delay differential equations,” in Advanced Differential Equations, pp. 287–
328, Elsevier, 2023.

[82] P. Eberhard and C. Bischof, “Automatic differentiation of numerical integration algorithms,”

Mathematics of Computation, vol. 68, pp. 717–732, Apr. 1999.

[83] A. Sei and W. W. Symes, “A note on consistency and adjointness for numerical schemes,”
Institute for Mathematics and Its Applications Preprint Series, Aug. 1995.

[84] F. Sapienza, J. Bolibar, F. Schäfer, B. Groenke, A. Pal, V. Boussange, P. Heimbach,

G. Hooker, F. Pérez, P.-O. Persson, and C. Rackauckas, “Differentiable Programming for
Differential Equations: A Review,” June 2024. arXiv:2406.09699 [[Link]], preprint.

[85] E. Huntley, “A note on the application of the matrix Riccati equation to the optimal control of
distributed parameter systems,” IEEE Transactions on Automatic Control, vol. 24, pp. 487–
489, June 1979.

[86] Z. Sirkes and E. Tziperman, “Finite Difference of Adjoint or Adjoint of Finite Difference?,”
AMS Monthly Weather Review, pp. 3373–3378, Dec. 1997.

[87] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ra-

madhan, and A. Edelman, “Universal Differential Equations for Scientific Machine Learning,”
Nov. 2021.

[88] H. Zhang, S. Abhyankar, E. Constantinescu, and M. Anitescu, “Discrete Adjoint Sensitivity

Analysis of Hybrid Dynamical Systems With Switching,” IEEE Transactions on Circuits and
Systems I: Regular Papers, vol. 64, pp. 1247–1259, May 2017.

[89] Y. Cao, S. Li, L. Petzold, and R. Serban, “Adjoint Sensitivity Analysis for Differential-
Algebraic Equations: The Adjoint DAE System and Its Numerical Solution,” SIAM Journal
on Scientific Computing, vol. 24, pp. 1076–1089, Jan. 2003.

[90] A. Griewank and A. Walther, Evaluating Derivatives: Principles and Techniques of

Algorithmic Differentiation. Society for Industrial and Applied Mathematics, second ed.,
Aug. 2008.

[91] I. Hiskens and M. Pai, “Trajectory sensitivity analysis of hybrid systems,” IEEE Transactions
on Circuits and Systems I: Fundamental Theory and Applications, vol. 47, pp. 204–220, Feb.
2000.

[92] X. Pan, W. Liang, J. Guo, X. Sun, and H. Chen, “Parameter Estimation by Considering Both
Voltage Dip and Phase Angle Jump for DFIG-Based WTGs in Distribution Grids,” IEEE
Transactions on Sustainable Energy, vol. 15, pp. 2170–2181, Oct. 2024.

[93] I. M. Granitsas, I. A. Hiskens, J. L. Mathieu, and G. S. Ledva, “Parameter Identifiability

and Estimation of Thermostatically Controlled Loads,” in 2023 IEEE Belgrade PowerTech,
pp. 01–06, June 2023.
 137

[94] R. Nagi, X. Huan, and Y. C. Chen, “Bayesian Inference of Parameters in Power System Dy-
namic Models Using Trajectory Sensitivities,” IEEE Transactions on Power Systems, vol. 37,
pp. 1253–1263, Mar. 2022.

[95] I. Hiskens, “Nonlinear dynamic model evaluation from disturbance measurements,” IEEE
Transactions on Power Systems, vol. 16, pp. 702–710, Nov. 2001.

[96] J. Ma, D. Han, R.-M. He, Z.-Y. Dong, and D. J. Hill, “Reducing Identified Parameters of
Measurement-Based Composite Load Model,” IEEE Transactions on Power Systems, vol. 23,
pp. 76–83, Feb. 2008.

[97] T. B. Nguyen, M. A. Pai, and I. A. Hiskens, “Sensitivity approaches for direct computation
of critical parameters in a power system,” International Journal of Electrical Power & Energy
Systems, vol. 24, pp. 337–343, June 2002.

[98] D. Chatterjee and A. Ghosh, “Transient Stability Assessment of Power Systems Containing
Series and Shunt Compensators,” IEEE Transactions on Power Systems, vol. 22, pp. 1210–
1220, Aug. 2007.

[99] I. Hiskens and J. Alseddiqui, “Sensitivity, Approximation, and Uncertainty in Power System
Dynamic Simulation,” IEEE Transactions on Power Systems, vol. 21, pp. 1808–1820, Nov.
2006.

[100] H. Choi, R. Elliott, D. Venkat, and D. Trudnowski, “Security Constrained Uncertainty In-
terval Estimation using Sensitivity Trajectories in Dynamical Systems,” in 2024 American
Control Conference (ACC), pp. 3411–3416, July 2024.

[101] K. Shubhanga and A. Kulkarni, “Determination of effectiveness of transient stability controls

using reduced number of trajectory sensitivity computations,” IEEE Transactions on Power
Systems, vol. 19, pp. 473–482, Feb. 2004.

[102] S. Q. Yuan and D. Z. Fang, “Robust PSS Parameters Design Using a Trajectory Sensitivity
Approach,” IEEE Transactions on Power Systems, vol. 24, pp. 1011–1018, May 2009.

[103] S.-M. Baek, J.-W. Park, and I. Hiskens, “Optimal Tuning for Linear and Nonlinear Pa-
rameters of Power System Stabilizers in Hybrid System Modeling,” in 2007 IEEE Industry
Applications Annual Meeting, pp. 1665–1672, Sept. 2007.

[104] A. Bidram, M.-e. Hamedani-golshan, and A. Davoudi, “Capacitor Design Considering First
Swing Stability of Distributed Generations,” IEEE Transactions on Power Systems, vol. 27,
pp. 1941–1948, Nov. 2012.

[105] D. Chatterjee and A. Ghosh, “TCSC control design for transient stability improvement of a
multi-machine power system using trajectory sensitivity,” Electric Power Systems Research,
vol. 77, pp. 470–483, Apr. 2007.

[106] A. Zamora-Cárdenas and C. R. Fuerte-Esquivel, “Multi-parameter trajectory sensitivity ap-

proach for location of series-connected controllers to enhance power system transient stabil-
ity,” Electric Power Systems Research, vol. 80, pp. 1096–1103, Sept. 2010.
 138

[107] M. Zima, P. Korba, and G. Andersson, “Power systems voltage emergency control ap-
proach using trajectory sensitivities,” in Proceedings of 2003 IEEE Conference on Control
Applications, 2003. CCA 2003., vol. 1, pp. 189–194 vol.1, June 2003.

[108] G. Hou and V. Vittal, “Trajectory Sensitivity Based Preventive Control of Voltage Instability
Considering Load Uncertainties,” IEEE Transactions on Power Systems, vol. 27, pp. 2280–
2288, Nov. 2012.

[109] Y. Xu, Z. Y. Dong, J. Zhao, Y. Xue, and D. J. Hill, “Trajectory sensitivity analysis on
the equivalent one-machine-infinite-bus of multi-machine systems for preventive transient
stability control,” IET Generation, Transmission & Distribution, vol. 9, no. 3, pp. 276–286,
2015.

[110] T. Nguyen and M. Pai, “Dynamic security-constrained rescheduling of power systems using
trajectory sensitivities,” IEEE Transactions on Power Systems, vol. 18, pp. 848–854, May
2003.

[111] S. Abhyankar, V. Rao, and M. Anitescu, “Dynamic security constrained optimal power flow
using finite difference sensitivities,” in 2014 IEEE PES General Meeting | Conference &
Exposition, pp. 1–5, July 2014. ISSN: 1932-5517.

[112] Y. Ma, V. Dixit, M. J. Innes, X. Guo, and C. Rackauckas, “A Comparison of Automatic

Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation
Solutions,” in IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–9,
Sept. 2021.

[113] M. Innes, A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. B. Shah, and W. Tebbutt, “A

Differentiable Programming System to Bridge Machine Learning and Scientific Computing,”
July 2019. arXiv:1907.07587 [cs], preprint.

[114] W. Moses and V. Churavy, “Instead of Rewriting Foreign Code for Machine Learning, Auto-
matically Synthesize Fast Gradients,” in Advances in Neural Information Processing Systems,
vol. 33, pp. 12472–12485, 2020.

[115] R. W. Kenyon, B. Wang, A. Hoke, J. Tan, and B. M. Hodge, “Comparison of Electromag-

netic Transient and Phasor Dynamic Simulations: Implications for Inverter Dominated Sys-
tems,” in 2023 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference
(ISGT), pp. 1–5, Jan. 2023.

[116] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, vol. 14 of Springer
Series in Computational Mathematics. Springer Berlin Heidelberg, 1996.

[117] D. Charles, “Off the Grid: Computer models that forecast overloaded power lines are holding
back U.S. solar and wind energy projects,” Science, vol. 381, no. 6662, pp. 1042–1045, 2023.

[118] “Value and Limitations of the Positive Sequence Generic Models of Renewable Energy Sys-
tems,” tech. rep., WECC Modeling and Validation Working Group, Dec. 2015.

[119] R. Anantharaman, Y. Ma, S. Gowda, C. Laughman, V. Shah, A. Edelman, and C. Rack-

auckas, “Accelerating Simulation of Stiff Nonlinear Systems using Continuous-Time Echo
State Networks,” Mar. 2021. arXiv:2010.04004 [cs, math], preprint.
 139

[120] A. Azmy and I. Erlich, “Identification of dynamic equivalents for distribution power networks
using recurrent ANNs,” in IEEE PES Power Systems Conference and Exposition, pp. 348–
353, Oct. 2004.
[121] L. El Ghaoui, F. Gu, B. Travacca, A. Askari, and A. Tsai, “Implicit Deep Learning,” SIAM
Journal on Mathematics of Data Science, vol. 3, pp. 930–958, Jan. 2021.
[122] F. Milano, Power System Modelling and Scripting. Springer Berlin Heidelberg, 2010.
[123] R. H. Park, “Two-reaction theory of synchronous machines generalized method of analysis-
part I,” Transactions of the American Institute of Electrical Engineers, vol. 48, no. 3, p. 716,
1929.
[124] A. Pal, A. Edelman, and C. Rackauckas, “Mixing Implicit and Explicit Deep Learning with
Skip DEQs and Infinite Time Neural ODEs (Continuous DEQs),” 2022. arXiv:2201.12240
[[Link]], preprint.
[125] S. Kim, W. Ji, S. Deng, Y. Ma, and C. Rackauckas, “Stiff neural ordinary differential equa-
tions,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 31, p. 093122, Sept.
2021.
[126] D. P. Kingma and J. L. Ba, “Adam: A Method for Stochastic Optimization,” in 3rd
International Conference for Learning Representations, (San Diego, CA), 2015.
[127] “SciML: Open Source Software for Scientific Machine Learning.” url: [Link]
[128] R. W. Kenyon, A. Sajadi, A. Hoke, and B.-M. Hodge, “Open-Source PSCAD Grid-Following
and Grid-Forming Inverters and A Benchmark for Zero-Inertia Power System Simulations,”
in 2021 IEEE Kansas Power and Energy Conference (KPEC), pp. 1–6, Apr. 2021.
[129] S. D’Arco, J. A. Suul, and O. B. Fosso, “A Virtual Synchronous Machine implementation
for distributed control of power converters in SmartGrids,” Electric Power Systems Research,
vol. 122, pp. 180–197, May 2015.
[130] C. Shen, A. P. Appling, P. Gentine, T. Bandai, H. Gupta, A. Tartakovsky, M. Baity-Jesi,
F. Fenicia, D. Kifer, L. Li, X. Liu, W. Ren, Y. Zheng, C. J. Harman, M. Clark, M. Farthing,
D. Feng, P. Kumar, D. Aboelyazeed, F. Rahmani, Y. Song, H. E. Beck, T. Bindas, D. Dwivedi,
K. Fang, M. Höge, C. Rackauckas, B. Mohanty, T. Roy, C. Xu, and K. Lawson, “Differentiable
modelling to unify machine learning and physical models for geosciences,” Nature Reviews
Earth & Environment, vol. 4, pp. 552–567, Aug. 2023.
[131] I. V. Nadal, J. Stiasny, and S. Chatzivasileiadis, “Integrating Physics-Informed Neural Net-
works into Power System Dynamic Simulations,” Apr. 2024.
[132] R. Bellman, “Dynamic Programming,” Science, vol. 153, pp. 34–37, July 1966.
[133] D. Widmann and C. Rackauckas, “DelayDiffEq: Generating Delay Differential Equation
Solvers via Recursive Embedding of Ordinary Differential Equation Solvers,” Aug. 2022.
arXiv:2208.12879 [[Link]], preprint.
[134] A. Bellen, M. Zennaro, A. Bellen, and M. Zennaro, Numerical Methods for Delay Differential
Equations. Numerical Mathematics and Scientific Computation, Oxford University Press,
Jan. 2013.
 140

[135] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” Neural Computation, vol. 9,
pp. 1735–1780, Nov. 1997.

[136] B. Lindemann, T. Müller, H. Vietz, N. Jazdi, and M. Weyrich, “A survey on long short-term
memory networks for time series prediction,” Procedia CIRP, vol. 99, pp. 650–655, Jan. 2021.

[137] A. Pal, “Lux: Explicit Parameterization of Deep Neural Networks in Julia,” Apr. 2023.
doi:10.5281/zenodo.7808904.

[138] G. Steinebach, “Construction of Rosenbrock–Wanner method Rodas5P and numerical bench-

marks within the Julia Differential Equations package,” BIT Numerical Mathematics, vol. 63,
p. 27, Apr. 2023.

[139] J. Stiasny and S. Chatzivasileiadis, “Physics-informed neural networks for time-domain sim-
ulations: Accuracy, computational cost, and flexibility,” Electric Power Systems Research,
vol. 224, p. 109748, Nov. 2023.

[140] G. K. W. Kenway, C. A. Mader, P. He, and J. R. R. A. Martins, “Effective adjoint approaches

for computational fluid dynamics,” Progress in Aerospace Sciences, vol. 110, p. 100542, Oct.
2019.

[141] R. H. Langland and N. L. Baker, “Estimation of observation impact using the NRL atmo-
spheric variational data assimilation adjoint system,” Tellus A: Dynamic Meteorology and
Oceanography, vol. 56A, Jan. 2004.

[142] A. Callejo, O. A. Bauchau, and B. Diskin, “Adjoint method for the sensitivity analysis of
composite beam cross-sections,” Computers & Structures, vol. 213, pp. 100–111, Mar. 2019.

[143] S. S. Gaikwad, L. Hascoet, S. H. K. Narayanan, L. Curry-Logan, R. Greve, and P. Heimbach,

“SICOPOLIS-AD v2: tangent linear and adjoint modeling framework for ice sheet modeling
enabled by automatic differentiation tool Tapenade,” Journal of Open Source Software, vol. 8,
p. 4679, Mar. 2023.

[144] I. Hiskens, “Power System Modeling for Inverse Problems,” IEEE Transactions on Circuits
and Systems I: Regular Papers, vol. 51, pp. 539–551, Mar. 2004.

[145] I. Hiskens and M. Pai, “Power system applications of trajectory sensitivities,” in 2002 IEEE
Power Engineering Society Winter Meeting. Conference Proceedings (Cat. No.02CH37309),
vol. 2, (New York, NY, USA), pp. 1200–1205, IEEE, 2002.

[146] S. G. Johnson, “Notes on Adjoint Methods for 18.335.” Apr. 2021.

[147] R. J. LeVeque, “1. Finite Difference Approximations,” in Finite Difference Methods for
Ordinary and Partial Differential Equations, Other Titles in Applied Mathematics, pp. 3–11,
Society for Industrial and Applied Mathematics, Jan. 2007.

[148] R. R. Barton, “Computing forward difference derivatives in engineering optimization,”

Engineering Optimization, vol. 20, pp. 205–224, Dec. 1992.

[149] R. Mathur, An analytical approach to computing step sizes for finite-difference derivatives.
PhD thesis, May 2012.
 141

[150] 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, pp. 5595–5637,
Jan. 2017.

[151] J. Schmidhuber, “Deep learning in neural networks: An overview,” Neural Networks, vol. 61,
pp. 85–117, Jan. 2015.

[152] F. Milano, “Semi-Implicit Formulation of Differential-Algebraic Equations for Transient Sta-

bility Analysis,” IEEE Transactions on Power Systems, vol. 31, pp. 4534–4543, Nov. 2016.

[153] S. Balay, S. Abhyankar, M. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Con-

stantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. Gropp, V. Hapla,
T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. Knepley, F. Kong, S. Kruger, D. May,
L. McInnes, R. Mills, L. Mitchell, T. Munson, J. Roman, K. Rupp, P. Sanan, J. Sarich,
B. Smith, S. Zampini, H. Zhang, and J. Zhang, “PETSc/TAO Users Manual (Rev. 3.20),”
Tech. Rep. ANL–21/39 Rev. 3-20, 2205494, 185712, Nov. 2023.

[154] H. Zhang, E. M. Constantinescu, and B. F. Smith, “PETSc TSAdjoint: A Discrete Ad-

joint ODE Solver for First-Order and Second-Order Sensitivity Analysis,” SIAM Journal on
Scientific Computing, vol. 44, pp. C1–C24, Feb. 2022.

[155] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet Classification with Deep Convo-
lutional Neural Networks,” in Advances in Neural Information Processing Systems, vol. 25,
2012.

[156] P. Mattson, V. J. Reddi, C. Cheng, C. Coleman, G. Diamos, D. Kanter, P. Micikevicius,

D. Patterson, G. Schmuelling, H. Tang, G.-Y. Wei, and C.-J. Wu, “MLPerf: An Industry
Standard Benchmark Suite for Machine Learning Performance,” IEEE Micro, vol. 40, pp. 8–
16, Mar. 2020.

[157] V. J. Reddi, “Benchmarking AI,” in Machine Learning Systems: Principles and Practices of
Engineering Artificially Intelligent Systems, Harvard University, 2024.

[158] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by Simulated Annealing,”

Science, vol. 220, pp. 671–680, May 1983.

[159] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and

A. Courville, “On the Spectral Bias of Neural Networks,” in Proceedings of the 36th
International Conference on Machine Learning, pp. 5301–5310, PMLR, May 2019.

[160] A. Jordana, J. Carpentier, and L. Righetti, “Learning Dynamical Systems from Noisy Sensor
Measurements using Multiple Shooting,” June 2021. arXiv:2106.11712 [cs, stat], preprint.

[161] C. Lataniotis, Data-driven uncertainty quantification for high-dimensional engineering

problems. PhD thesis, ETH Zurich, Nov. 2019.

[162] A. Kendall and Y. Gal, “What Uncertainties Do We Need in Bayesian Deep Learning for
Computer Vision?,” in Advances in Neural Information Processing Systems, vol. 30, 2017.

[163] S. Stock, J. Stiasny, D. Babazadeh, C. Becker, and S. Chatzivasileiadis, “Bayesian Physics-

Informed Neural Networks for Robust System Identification of Power Systems,” in 2023 IEEE
Belgrade PowerTech, pp. 1–6, June 2023.
 142

[164] J. Thiyagalingam, M. Shankar, G. Fox, and T. Hey, “Scientific machine learning benchmarks,”
Nature Reviews Physics, vol. 4, pp. 413–420, June 2022.

[165] S. A. Desai, Physics-Informed Neural Networks for Data-Efficient Learning. PhD thesis,
University of Oxford, 2021.
 ProQuest Number: 31637257

INFORMATION TO ALL USERS

The quality and completeness of this reproduction is dependent on the quality
and completeness of the copy made available to ProQuest.

Distributed by
ProQuest LLC a part of Clarivate ( 2024 ).
Copyright of the Dissertation is held by the Author unless otherwise noted.

This work is protected against unauthorized copying under Title 17,

United States Code and other applicable copyright laws.

This work may be used in accordance with the terms of the Creative Commons license
or other rights statement, as indicated in the copyright statement or in the metadata
associated with this work. Unless otherwise specified in the copyright statement
or the metadata, all rights are reserved by the copyright holder.

ProQuest LLC
789 East Eisenhower Parkway
Ann Arbor, MI 48108 USA