# Learning Networked Dynamical System Models with
Weak Form and Graph Neural Networks

Yin Yu111Graduate Student Research Assistant, Department of Aerospace
Engineering, Student Member AIAA<EMAIL_ADDRESS>The Pennsylvania State
University, University Park, PA, 16802  Daning Huang222Assistant Professor,
Department of Aerospace Engineering, Member AIAA<EMAIL_ADDRESS>The
Pennsylvania State University, University Park, PA, 16802  Seho
Park333Graduate Student Research Assistant, Department of Mechanical
Engineering<EMAIL_ADDRESS>The Pennsylvania State University, University
Park, PA, 16802  Herschel C. Pangborn444Assistant Professor, Department of
Mechanical Engineering, Member, AIAA<EMAIL_ADDRESS>The Pennsylvania
State University, University Park, PA, 16802

###### Abstract

This paper presents a sequence of two approaches for the data-driven control-
oriented modeling of networked systems, i.e., the systems that involve many
interacting dynamical components. First, a novel deep learning approach named
the weak Latent Dynamics Model (wLDM) is developed for learning generic
nonlinear dynamics with control. Leveraging the weak form, the wLDM enables
more numerically stable and computationally efficient training as well as more
accurate prediction, when compared to conventional methods such as neural
ordinary differential equations. Building upon the wLDM framework, we propose
the weak Graph Koopman Bilinear Form (wGKBF) model, which integrates geometric
deep learning and Koopman theory to learn latent space dynamics for networked
systems, especially for the challenging cases having multiple timescales. The
effectiveness of the wLDM framework and wGKBF model are demonstrated on three
example systems of increasing complexity - a controlled double pendulum, the
stiff Brusselator dynamics, and an electrified aircraft energy system. These
numerical examples show that the wLDM and wGKBF achieve superior predictive
accuracy and training efficiency as compared to baseline models. Parametric
studies provide insights into the effects of hyperparameters in the weak form.
The proposed framework shows the capability to efficiently capture control-
dependent dynamics in these systems, including stiff dynamics and multi-
physics interactions, offering a promising direction for learning control-
oriented models of complex networked systems.

## Nomenclature

Roman Symbols

$a$, $b$ = | Time interval bounds of test functions
---|---
$A$, $B$ = | Brusselator system parameters
${\mathbf{A}}$, ${\mathbf{B}}_{i}$ = | System matrices of bilinear dynamical systems
${\mathbf{C}}$, ${\mathbf{D}}$ = | Weight matrices of integration approximation
${\mathcal{D}}$ = | Dataset
${\mathcal{E}}$ = | Edges of a graph
${\mathbf{f}}$ = | System dynamics
${\mathbf{f}}_{E}$, ${\mathbf{f}}_{P}$, ${\mathbf{f}}_{D}$ = | Encoder, processor, and decoder of neural network
${\mathbf{g}}$ = | Observation function
${\mathcal{G}}$ = | Graph
${\mathbf{h}}$, ${\mathbf{H}}$ = | Node feature vectors and their collection
$J$ = | Optimization objective
${\mathcal{K}}_{t}$ = | Koopman operator
${\mathcal{L}}$ = | Loss function
$L$ = | Number of steps in a trajectory
${\mathbf{L}}$ = | Laplacian matrix of a graph
$N$ = | Time horizon, i.e., number of steps in a time interval of a test function
$N_{E}$, $N_{P}$, $N_{D}$ = | Number of layers for encoder, processor, and decoder
${\mathcal{N}}_{I}$ = | Neighbors of node $I$
$n_{I},m_{I}$ = | Dimensions of states and inputs at node $I$
$n_{N}$ = | Number of nodes
$n_{T}$ = | Number of trajectories
$P^{(\alpha,\beta)}$ = | Jacobi polynomial
$S$ = | Latent dimension
$\tilde{S}$ = | Latent dimension of a node
$T_{k}$ = | The $k$th order Chebyshev polynomial
$t$ = | Time
$\Delta t$ = | Time step
${\mathbf{u}}$, ${\mathbf{U}}$ = | Control input vector and its collection
${\mathcal{V}}$ = | Nodes of a graph
$\mathbf{v}$, $v_{i}$ = | Numerical integration weights and the $i$th element
${\mathbf{w}}$, ${\mathbf{W}}$ = | Latent representation vector and its collection
${\mathbf{x}}$, ${\mathbf{X}}$ = | State vector and its collection
${\mathbf{y}}$, ${\mathbf{Y}}$ = | System observation and its collection
$\hat{{\mathbf{y}}}$, $\hat{{\mathbf{Y}}}$ = | Measured observation and its collection

Greek Symbols

$\alpha$, $\beta$ = | Jacobi polynomial parameters
---|---
$\gamma$ = | Penalty parameter
$\delta_{ij}$ = | Kronecker delta
${\boldsymbol{\Theta}}$ = | Trainable parameters of the model
${\boldsymbol{\lambda}}$ = | Adjoint state vector
$\rho$ = | Weight function for the weak form test functions
${\boldsymbol{\sigma}}$ = | Activation function
$\phi$ = | Test functions

Others

$\hat{\Box}$ = | Measured observation for training
---|---
$\Box_{I}$ = | Index for nodes of a graph
$\Box_{i}$ = | Index for time or element
$n_{\Box}$ = | Dimension of the quantity $\Box$

Acronyms

BDF = | Backward Differentiation Formula
---|---
biDMDc = | bilinear Dynamic Mode Decomposition with control
FCNN = | Fully-Connected Neural Network
GCL = | Graph Convolution Layer
GNN = | Graph Neural Network
IO = | Integral Order
KBF = | Koopman Bilinear Form
LSTM = | Long Short-Term Memory
MP = | Message Passing
MPC = | Model Predictive Control
NODE = | Neural Ordinary Differential Equation
NRMSE = | Normalized Root Mean Squared Error
ODE = | Ordinary Differential Equation
PO = | Polynomial Order
PReLU = | Parametric Rectified Linear Unit
wLDM = | weak Latent Dynamics Model
wGKBF = | weak Graph Koopman Bilinear Form
SINDy = | Sparse Identification of Nonlinear Dynamics

##  1 Introduction

Powertrain electrification is a megatrend in vehicle design, leading to
enhanced performance, new paradigms in safety and sustainability, and
decreased lifetime costs [1, 2, 3, 4, 5]. Electrified vehicle energy systems
are characterized by complex nonlinear dynamics, having tightly coupled
interactions among system components across multiple energy domains (e.g.,
electrical, thermal, mechanical). This is typical of many networked systems,
which consist of many components whose dynamics may affect each other.
Networked systems usually demonstrate two unique features: sparsity and
heterogeneity [6]. Sparsity refers to a coupling structure in which each
component only interacts with a small number of other components. This is a
feature that can be leveraged to simplify the modeling and control of
networked systems. Heterogeneity implies that the components may be governed
by different dynamics and interact with a different number of other
components. Importantly, the dynamics across all the components may span a
wide range of timescales, thus producing a stiff system as a whole.
Heterogeneity, together with the potential for stiffness, typically pose a
challenge to the modeling and control of networked systems.

Control methods for nonlinear systems have been explored and demonstrated
extensively in the literature [7, 8]. A notable example for the networked
systems is hierarchical Model Predictive Control (MPC) using graph-based
models [9, 10], which is particularly effective at coordinating decision-
making across energy systems with multiple timescales. The performance and
robustness of such methods are highly sensitive to model error, yet an
accurate model is usually difficult and time-consuming to obtain, especially
for complex systems such as the electrified vehicle energy systems. For
example, even if physics-based equations are employed to model the energy
system, there may still be unknown functional terms and parameters, e.g.,
electrical efficiencies and heat transfer coefficients, that need to be
estimated empirically. In addition, the physics-based model may have nonlinear
terms that are inconvenient for controller design [11].

This paper aims to develop methods for creating data-driven models of
networked systems that are suitable for use in model-based control design. The
key questions are twofold: (1) What is the appropriate form for such a
control-oriented model, and (2) how to identify the model from data?

##### Choice of Model Form

A generic answer to the first question is that a control-oriented model should
be sufficiently compact, i.e., low-dimensional, for computational efficiency
and hardware deployment. Ideally the model should possess a structure that a
control method can leverage [7, 8], with classical examples being transfer
functions and linear time-invariant state space models. Relaxing the linearity
requirement, another example is control-affine models that are only linear in
the control inputs [8]. A control-oriented modeling approach undergoing recent
revival is the bilinear model [12], which has a special, mildly nonlinear
control-affine form whose autonomous component is linear in states and control
component is bilinear in both states and control. One advantage of bilinear
models over generic control-affine models is the availability of various
specialized, computationally efficient control methods [13, 14]. A second
advantage is the theoretical foundation in Koopman theory [15, 16, 17, 18],
particularly the Koopman Bilinear Form (KBF) [19, 20, 21], which provides
rigorous conditions for the existence of bilinear models for the
representation of a broad range of nonlinear dynamics [22, 23, 24]. In Koopman
theory, nonlinear dynamics are lifted to a higher-dimensional space where the
dynamics are governed by a linear operator. Notable data-driven methods
include extended DMD with control (eDMDc) [25] and bilinear DMD with control
(biDMDc) [19]. These conventional Koopman methods typically use a library of
feature functions, e.g., polynomials, to construct the mapping between
physical and Koopman states. Such library-based mapping poses significant
challenges, particularly in maintaining predictive accuracy and
generalizability, as well as scalability to high-dimensional physical states.
Recent approaches leverage deep neural networks to learn this mapping [26,
27], and the resulting models show improvements over conventional Koopman
methods that use library-based mapping in model compactness, predictive
accuracy, and control performance [28, 29].

Despite the success of data-driven methods for identifying Koopman operators,
the extension of KBF to networked systems remains underexplored due to the
challenge in the data representation of networked systems. It is natural to
use graphs to describe the networked systems, where the nodes represent the
components and the edges represent the interactions. However, traditional
machine learning algorithms, such as convolutional neural networks, are
defined for structured data, e.g., vectors and matrices, and struggle with
unstructured data, e.g., those defined on graphs [30]. The unstructured data
do not conform to a tabular structure such as matrices, and conventional tools
such as linear algebra cannot be directly applied. Graph Neural Networks
(GNNs) [31, 32] have emerged as a powerful tool for modeling data defined on
graphs, learning from both node-level and edge-level features to capture local
and global patterns. They handle unstructured data that is common in real-
world applications. In particular, incorporating domain knowledge and system
properties into graph structures has been found to better characterize the
underlying dynamics and interactions, leading to more accurate and efficient
dynamical models [33, 34]. GNN-based models have been successfully applied to
complex physical simulations, such as fluid dynamics [35, 36], weather systems
[37, 38], robotic control systems [39, 40], power systems [41, 34], etc.
Hence, the extension of KBF with GNNs is worth exploring, as it holds
potential to produce a computationally efficient control-oriented model that
encodes interactions among the components of a complex networked system while
efficiently handling sparse connectivity and heterogeneous dynamics.

##### Learning Dynamics from Data

After selecting an appropriate model form, the next step is to learn the
system dynamics from available data. The algorithms for learning dynamics are
broadly categorized into discrete-time and continuous-time formulations.
Discrete-time methods, such as autoregressive models [42] and recurrent neural
networks [43], map the system state from one time step to the next. These
methods are easy to implement but may have highly nonlinear model forms that
are inconvenient for controller design. Furthermore, it is either difficult or
computationally inefficient to adapt the model to different time step sizes.
This limitation poses significant challenges in various control scenarios,
such as hierarchical MPC, where the ability to handle multiple time step sizes
is essential [44]. In contrast, continuous-time methods learn a model to map
the states to their time derivatives, offering a more natural and often more
compact representation of underlying dynamical processes, especially when the
underlying physics are governed by ordinary differential equations (ODEs). In
control applications, the ODE may be discretized differently to maximize the
computational efficiency to resolve the dynamics at different timescales.

One major category of continuous-time methods is the “derivative-based”
formulation, which directly fits a model to match the time derivative of the
states [45]. This formulation usually results in a linear or nonlinear
constraint-free regression problem. A representative and widely used method in
this category is Sparse Identification of Nonlinear Dynamics (SINDy) [46].
SINDy builds upon the idea of sparse regression, seeking to find a minimal set
of nonlinear functions that best describe the system dynamics via an ODE. The
SINDy algorithm is relatively easy to implement as it only requires an
iterative process of linear regressions. However, one limitation of SINDy, and
derivative-based formulations in general, is sensitivity to the error of
estimated state time derivatives, where inaccurate derivative estimates can
result in the identification of incorrect or suboptimal governing equations.
The weak SINDy (wSINDy) method [47] addresses this issue by re-formulating the
linear regression problem using the so-called weak form. The weak form uses a
set of test functions to replace time derivatives with numerical integrals,
thus bypassing the need for the derivative estimates that may be numerically
unstable. The wSINDy has been shown to significantly improve model performance
in identifying system equations from data by avoiding the amplification of
errors resulting from numerical differentiation or noisy measurements [48].
Despite this improvement, wSINDy still heavily relies on a fixed library of
functions, usually polynomials, that suffer from the curse of dimensionality.
This characteristic fundamentally limits its capacity to model complex, high-
dimensional systems, such as the networked systems.

Deep learning-based models can automatically discover low-dimensional latent
representations of high-dimensional systems, overcoming the limitations of
fixed function libraries and enabling more efficient learning and prediction
of complex, nonlinear behaviors [49]. Learning dynamics in a latent space
offers several key advantages. First, it significantly reduces the
dimensionality of the problem, leading to improved computational efficiency,
particularly for systems with high-dimensional state spaces or intricate
interactions. This reduction allows for faster training and inference, making
real-time applications more feasible [50]. Second, for high-dimensional
dynamical system, by projecting the system dynamics onto a lower-dimensional
manifold, latent space representations can capture the essential features and
underlying structure of the system, filtering out noise and irrelevant details
[51]. This focus on core dynamics can lead to models with improved
generalizability and predictive accuracy. There is indeed an extension of
wSINDy to a latent dynamics formulation using dimension reduction techniques,
such as autoencoders [52]. However, in this approach, the dimension reduction
and dynamics learning are performed in a sequential manner, and the dynamics
still employ a fixed library of functions, which may suffer from lack of
expressiveness. Furthermore, the latent dynamics extension [52] did not
account for control inputs.

Hence, an alternative continuous-time formulation that is more promising for
networked systems is Neural Ordinary Differential Equations (NODEs) [53] and
their graph-based variants [54, 55, 56, 57]. NODEs are examples of deep
learning models that learn continuous-time dynamics in a latent space, and
graph-based NODEs extend this idea by incorporating graph structures to
capture interactions among the components of a networked system. The NODE
methods belong to another category of continuous-time methods, the
“trajectory-based” formulation, which fits a model so that the trajectories
generated by the model match with those in the dataset. This formulation
results in an ODE-constrained problem that is usually highly nonlinear and
non-convex and thus more difficult to solve than the regression problem in the
derivative-based formulation. Furthermore, the gradient computation of NODEs,
as needed in the learning process, can be numerically unstable and thus
computationally expensive and memory-intensive [58, 59]. This drawback is
particularly pronounced in the case of stiff systems [60], which may result in
completely incorrect models or even the failure of learning. Hence the
bottleneck of stiffness is highly undesirable for the learning of dynamics for
networked systems, which are often stiff due to including multiple timescales.

##### Contributions

To bridge the gaps in the literature, we propose a novel and systematic
approach for data-driven control-oriented modeling of networked systems. The
major contributions of this paper are as follows:

  1. 1.

We introduce the weak Latent Dynamics Model (wLDM), a framework for learning
generic nonlinear dynamics. This framework extends the NODE method from the
trajectory-based formulation with the weak form concept from wSINDy, i.e., the
derivative-based formulation, to eliminate the bottleneck of stiffness in the
learning of NODEs and enhance the accuracy of the learned dynamical models.

  2. 2.

Building upon the general wLDM framework, we propose the weak Graph Koopman
Bilinear Form (wGKBF) model as a control-oriented model for networked systems.
The wGKBF model captures the complex dynamics of networked systems while
maintaining a relatively simple bilinear structure, leveraging the
expressiveness of deep learning for sparse and heterogeneous network dynamics.

  3. 3.

We demonstrate the effectiveness of the proposed wLDM method, including its
wGKBF variant, through three example systems of increasing complexity: a
controlled double pendulum, a stiff Brusselator dynamics model, and an
electrified aircraft energy system. The results highlight the capability of
the proposed framework to produce accurate and computationally efficient
control-oriented models for complex networked systems.

The rest of the paper is organized as follows. Section 2 states the general
problem of learning dynamics from data and discusses the existing algorithms
and their limitations. Section 3 presents the proposed wLDM framework for
solving the general problem. Section 4 adapts the wLDM framework to develop
wGKBF for the control-oriented modeling of networked systems. Section 5
demonstrates the efficacy of wLDM and wGKBF models via three numerical
examples. Lastly, the conclusions are provided in Section 6.

##  2 Problem Statement and Preliminaries

In this section, we formally state the problem of inferring the parameters of
a nonlinear dynamical system from data. Then, we briefly discuss existing
algorithms for solving such problems, namely the direct and indirect methods,
and their limitations, which serves as the basis for the proposed algorithms.

###  2.1 Problem Statement

Consider a nonlinear dynamical system

$\displaystyle\dot{{\mathbf{x}}}$
$\displaystyle={\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})$
(1a) $\displaystyle{\mathbf{y}}$
$\displaystyle={\mathbf{g}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})$
(1b) $\displaystyle{\mathbf{x}}(t_{1})$
$\displaystyle={\mathbf{x}}_{0}({\boldsymbol{\Theta}})$ (1c)

where ${\mathbf{x}}\in{\mathbb{R}}^{n_{x}}$ is the state,
${\mathbf{u}}\in{\mathbb{R}}^{n_{u}}$ the control,
${\mathbf{y}}\in{\mathbb{R}}^{n_{y}}$ the observation, and
${\boldsymbol{\Theta}}\in{\mathbb{R}}^{n_{\theta}}$ the model parameters.
Functions ${\mathbf{f}}$ and ${\mathbf{g}}$ can depend on time $t$, however we
omit this for simplicity of exposition. The initial condition may be partially
or fully unknown, and Eq. (1c) provides an estimate. The goal is to infer
parameters ${\boldsymbol{\Theta}}$ of the potentially nonlinear models
${\mathbf{f}}$ and ${\mathbf{g}}$ from a dataset
${\mathcal{D}}=\\{({\mathbf{U}}^{(i)},\hat{\mathbf{Y}}^{(i)})\\}_{i=1}^{n_{T}}$,
where $n_{T}$ is the total number of trajectories. The $i$th trajectory
consists of a sequence of control inputs
${\mathbf{U}}^{(i)}=[{\mathbf{u}}_{1}^{(i)},{\mathbf{u}}_{2}^{(i)},\dots,{\mathbf{u}}_{L^{(i)}}^{(i)}]\in{\mathbb{R}}^{n_{u}\times
L^{(i)}}$ at time stations $\\{t_{1},t_{2},\cdots,t_{L}\\}$ and corresponding
measurements
$\hat{\mathbf{Y}}^{(i)}=[\hat{\mathbf{y}}_{1}^{(i)},\hat{\mathbf{y}}_{2}^{(i)},\dots,\hat{\mathbf{y}}_{L^{(i)}}^{(i)}]\in{\mathbb{R}}^{n_{y}\times
L^{(i)}}$, where $L^{(i)}$ is the length of the $i$th trajectory. Without loss
of generality, we formulate the optimization problem for learning the dynamics
using a single trajectory:

$\displaystyle\min_{\boldsymbol{\Theta}}$ $\displaystyle\
J({\boldsymbol{\Theta}};{\mathcal{D}})=\int_{t_{1}}^{t_{L}}\left\lVert{\mathbf{g}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})-\hat{\mathbf{y}}\right\rVert^{2}dt$
(2a) s.t. $\displaystyle\
\dot{{\mathbf{x}}}={\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}}),\quad{\mathbf{x}}(t_{1})={\mathbf{x}}_{0}({\boldsymbol{\Theta}})$
(2b)

This is the trajectory-based formulation as categorized in the introduction.

###  2.2 Available Algorithms

To solve Eq. (2), the continuous-time optimization problem is transformed into
a discrete-time approximation through discretization methods, converting the
infinite-dimensional problem into a finite-dimensional one. There are two main
categories of discretization methods [61]: Direct methods, also known as
“discretize-then-optimize”, and indirect methods, or “optimize-then-
discretize”. These two methods are briefly discussed and compared next.

####  2.2.1 Direct Methods

In the direct approach, the differential constraints are discretized, and the
continuous-time dynamics and the optimization problem are transformed into a
discrete form. This discretization process typically involves approximating
the system dynamics at a series of discrete time points using numerical
integration methods. For illustrative purposes, applying a simple first-order
discretization to Eq. (2) yields a discretized form,

$\displaystyle\min_{\boldsymbol{\Theta}}$ $\displaystyle\
J_{d}({\boldsymbol{\Theta}};{\mathcal{D}})=\sum_{i=1}^{L}\left\lVert{\mathbf{g}}({\mathbf{x}}_{i},{\mathbf{u}}_{i};{\boldsymbol{\Theta}})-\hat{\mathbf{y}}_{i}\right\rVert^{2}\Delta
t$ (3a) s.t. $\displaystyle\ \frac{{\mathbf{x}}_{i}-{\mathbf{x}}_{i-1}}{\Delta
t}={\mathbf{f}}({\mathbf{x}}_{i},{\mathbf{u}}_{i};{\boldsymbol{\Theta}}),\
i=2,\dots,L,\quad{\mathbf{x}}(t_{1})={\mathbf{x}}_{0}({\boldsymbol{\Theta}})$
(3b)

where $\Delta t$ is the time step between consecutive measurements, and
${\mathbf{x}}_{i}$ and ${\mathbf{u}}_{i}$ represent the state and control
input at the $i$-th time step, respectively. By discretizing the continuous-
time problem, the direct approach transforms the dynamical system into a set
of algebraic equations, which can be solved using standard numerical
optimization techniques. The optimization process seeks the parameters
${\boldsymbol{\Theta}}$ that minimize the objective function
$J_{d}({\boldsymbol{\Theta}};{\mathcal{D}})$ while satisfying the discretized
dynamics constraints.

In practice, more sophisticated discretization schemes, such as pseudospectral
methods [62], could be used to reduce the discretization error. However, due
to the inevitable discretization errors in Eq. (3b), the direct method can
suffer from the accumulation of approximation errors over time, which may lead
to suboptimal solutions. The accuracy and computational efficiency of the
optimization process depends heavily on the choice of the discretization
scheme and the time step size.

####  2.2.2 Indirect Methods

Conversely, the indirect approach starts with the necessary optimality
conditions of the continuous-time optimization problem before any
discretization. The optimality conditions are typically obtained by applying
Pontryagin’s Maximum Principle (PMP) [63], which introduces the adjoint state
${\boldsymbol{\lambda}}(t)\in{\mathbb{R}}^{n_{x}}$ and the Hamiltonian
$H({\mathbf{x}},{\mathbf{u}},{\boldsymbol{\lambda}},t)={\boldsymbol{\lambda}}^{T}{\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})+l({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})$,
with
$l({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})=\left\lVert{\mathbf{g}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})-\hat{\mathbf{y}}\right\rVert^{2}$.
PMP states that for an optimal solution, the following conditions must hold:

$\displaystyle\dot{{\mathbf{x}}}$ $\displaystyle=\frac{\partial
H}{\partial{\boldsymbol{\lambda}}}={\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}}),\quad{\mathbf{x}}(t_{1})={\mathbf{x}}_{0}$
(4a) $\displaystyle\dot{{\boldsymbol{\lambda}}}$
$\displaystyle=-\frac{\partial
H}{\partial{\mathbf{x}}}=-{\boldsymbol{\lambda}}^{T}\frac{\partial{\mathbf{f}}}{\partial{\mathbf{x}}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})-\frac{\partial
l}{\partial{\mathbf{x}}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}}),\quad{\boldsymbol{\lambda}}(t_{L})=\frac{\partial
J}{\partial{\mathbf{x}}(t_{L})}$ (4b) $\displaystyle 0$
$\displaystyle=\frac{\partial
H}{\partial{\mathbf{u}}}={\boldsymbol{\lambda}}^{T}\frac{\partial{\mathbf{f}}}{\partial{\mathbf{u}}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})+\frac{\partial
l}{\partial{\mathbf{u}}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})$
(4c)

where Eq. (4b) governs the evolution of the adjoint state, and is solved
backward in time with the terminal condition of $\frac{\partial
J}{\partial{\mathbf{x}}(t_{L})}$; the value of $l$ at any time may be
interpolated from the given discrete measurements. Equation (4) defines a two
point boundary value problem.

There are two popular approaches for solving the boundary value problem: The
shooting method and the collocation method [64]. The shooting method
iteratively solves the problem by guessing the initial adjoint state and
updating it based on the mismatch between the final state and the desired
terminal condition. The collocation method discretizes the continuous-time
optimality conditions using numerical schemes and solves the resulting
nonlinear algebraic equation. One of the main advantages of the indirect
approach is that it inherently handles the continuous nature of the dynamics
and adjoint equations, leading to more accurate gradients and solutions as
compared to the direct approach. However, the indirect approach often suffers
from numerical difficulties during optimization, such as the need for good
initial guesses for the adjoint state and system parameters, and the
sensitivity to the choice of numerical schemes for discretizing the optimality
conditions.

####  2.2.3 Applications in Machine Learning

In the context of machine learning, when the trajectory optimization problem
is chosen to learn dynamics from data, the indirect method is usually selected
due to its guarantee optimality guarantees. In this case, the learning
algorithm falls under the category of NODEs. The training of NODEs has been
implemented using shooting methods [53]. While NODEs provide a continuous-time
modeling framework and convert the constrained optimization problem into an
unconstrained one, training them can be computationally expensive and memory-
intensive for solving the adjoint ODE. Particularly for stiff systems,
solution of the adjoint ODE may be numerically unstable and can even lead to
failure in convergence of model learning [60]. Moreover, the numerical
integration process can introduce approximation errors and lead to numerical
instabilities. These challenges in existing algorithms motivate the
exploration of alternative approaches, such as the weak form, as discussed in
the following section.

##  3 Formulation of Weak Latent Dynamics Model

In this section, we introduce the weak form approach and propose the weak
Latent Dynamics Model (wLDM), a general deep learning architecture that learns
the dynamics of a system in a latent space. We detail the key components of
the wLDM and discuss its advantages over conventional methods.

###  3.1 Weak Form for Learning Dynamics

The original formulation of weak form for learning dynamics in Ref. [47] is
developed for derivative-based formulation assuming that the states are known;
this results in a linear regression problem, possibly with regularization
terms. We extend the approach to the NODE methods, i.e., trajectory-based
formulation, where the states are unknown. This results in a nonlinear
regression problem, where the mapping between the physical and latent states
as well as the latent dynamics are learned together in an iterative manner.
Essentially, the NODE methods are reformulated using the weak form to provide
a computationally efficient and numerically stable framework for
identification of the nonlinear dynamics.

####  3.1.1 Representation of Dynamics by Weak Form

Consider a set of $C^{1}$-continuous test functions
$\\{\phi_{i}(t)\\}_{i=1}^{n_{\phi}}$. A test function $\phi(t)$ is defined on
a time interval $[a,b]$, such that $\phi(t)=0$ for $t\notin[a,b]$ and $\phi$
satisfies the boundary conditions:

$\phi(a)=\phi(b)=0$ (5)

The ODE for the dynamics in Eq. (1) can be converted to the weak form by
multiplying both sides of the equation by a test function $\phi(t)$ and
integrating over the time interval $[a,b]$:

$\int_{a}^{b}\phi(t)\dot{{\mathbf{x}}}(t)dt=\int_{a}^{b}\phi(t){\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})dt$
(6)

Using integration by parts and the boundary conditions in Eq. (5), we can
eliminate the time derivative of ${\mathbf{x}}$ on the left-hand side of Eq.
(6):

$\int_{a}^{b}\phi(t)\dot{{\mathbf{x}}}(t)dt=\phi(t){\mathbf{x}}(t)\Big{|}_{a}^{b}-\int_{a}^{b}\dot{\phi}_{i}(t){\mathbf{x}}(t)dt=-\int_{a}^{b}\dot{\phi}_{i}(t){\mathbf{x}}(t)dt$
(7)

Substituting Eq. (7) into Eq. (6), we obtain the weak form of the dynamics:

$-\int_{a}^{b}\dot{\phi}(t){\mathbf{x}}(t)dt=\int_{a}^{b}\phi(t){\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})dt$
(8)

The key difference of Eq. (8) when compared to the original dynamics in Eq.
(1) is that $\dot{{\mathbf{x}}}$ no longer appears. This would have to be
estimated from the data of ${\mathbf{x}}$ and the estimation error becomes
part of the model error that lead to incorrect learning of dynamics. In the
weak form of Eq. (8), the time derivative is transferred to the test function
$\phi$ that is chosen by the user and the exact expression of $\dot{\phi}$ is
known. Hence, the use of a test function bypasses the need to estimate
$\dot{{\mathbf{x}}}$ and eliminates the model error due to derivative
estimation.

The weak form in Eq. (8) allows us to reformulate the optimization problem in
Eq. (2) as:

$\displaystyle\min_{{\boldsymbol{\Theta}}}$ $\displaystyle\
J({\boldsymbol{\Theta}};{\mathcal{D}})=\int_{t_{1}}^{t_{L}}\left\lVert{\mathbf{g}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})-\hat{\mathbf{y}}\right\rVert^{2}dt$
(9a) s.t. $\displaystyle\
-\int_{a_{i}}^{b_{i}}\dot{\phi}_{i}(t){\mathbf{x}}(t)dt=\int_{a_{i}}^{b_{i}}\phi_{i}(t){\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})dt,\quad
i=1,\dots,n_{\phi}$ (9b)

As the number of test functions $n_{\phi}$ approaches infinity and the union
of intervals $[a_{i},b_{i}]$ cover $[t_{1},t_{L}]$, the weak form constraints
in Eq. (9b) converge to the strong form of the dynamics, resulting in exact
satisfaction of the system dynamics. However, when working with a finite
number of test functions, Eq. (9b) only “weakly” satisfies the true dynamics,
hence the name “weak form”. The accuracy of the approximation depends on a
number of factors that are discussed later in Sec. 3.1.3.

####  3.1.2 Discretization and Integration Matrices

Next, we discretize the integrals in Eq. (9) based on the dataset
${\mathcal{D}}$, and convert the optimization to a conventional equality-
constrained problem. For any integrable function $\eta(t)$, given data pairs
$\\{(t_{k},\eta(t_{k}))\\}_{k=1}^{N}$, $a\leq t_{1}<t_{2}<\cdots<t_{N}\leq b$,
we can perform numerical integration as follows:

$\int_{a}^{b}\eta(t)dt\approx\sum_{i=1}^{N}v_{k}\eta(t_{k})\equiv{\boldsymbol{\eta}}^{T}\mathbf{v}$
(10)

where $v_{k}$ are weights at time stations $t_{k}$, and

$\displaystyle{\boldsymbol{\eta}}=[\eta(t_{1}),\eta(t_{2}),\cdots,\eta(t_{N})]^{T},\quad\mathbf{v}=[v_{1},v_{2},\cdots,v_{N}]^{T}$
(11)

Here the number of steps $N$ in the interval of the integral will be referred
to as the “time horizon” for conciseness.

The objective function in Eq. (9a) can be approximated using numerical
integration and the weighted Frobenius norm $\left\lVert\cdot\right\rVert_{v}$
as follows:

$\displaystyle J({\boldsymbol{\Theta}};{\mathcal{D}})$
$\displaystyle=\int_{t_{1}}^{t_{L}}\left\lVert{\mathbf{g}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})-\hat{\mathbf{y}}\right\rVert^{2}dt$
(12)
$\displaystyle\approx\sum_{k=1}^{N}v_{k}\left\lVert{\mathbf{y}}_{k}-\hat{\mathbf{y}}_{k}\right\rVert^{2}\equiv\left\lVert{\mathbf{Y}}-\hat{\mathbf{Y}}\right\rVert_{v}^{2}$

where
${\mathbf{y}}_{k}={\mathbf{g}}({\mathbf{x}}_{k},{\mathbf{u}}_{k};{\boldsymbol{\Theta}})$.
The weighted Frobenius norm is defined as:

$\left\lVert{\mathbf{Y}}-\hat{\mathbf{Y}}\right\rVert_{v}^{2}=\mathrm{tr}(({\mathbf{Y}}-\hat{\mathbf{Y}})\mathbf{V}({\mathbf{Y}}-\hat{\mathbf{Y}})^{T})$
(13)

where
${\mathbf{Y}}=[{\mathbf{y}}_{1},{\mathbf{y}}_{2},\dots,{\mathbf{y}}_{N}]$,
$\hat{\mathbf{Y}}=[\hat{\mathbf{y}}_{1},\hat{\mathbf{y}}_{2},\dots,\hat{\mathbf{y}}_{N}]$,
$\mathbf{V}=\mathrm{Diag}(v)$, and $\mathrm{tr}(\cdot)$ denotes the matrix
trace.

Similarly, for the constraints in Eq. (9b), the left-hand side of the $i$th
constraint is approximated as

$-\int_{a}^{b}\dot{\phi}_{i}{\mathbf{x}}dt\approx-\sum_{k=1}^{N}v_{k}\dot{\phi}_{i}(t_{k}){\mathbf{x}}_{k}\equiv{\mathbf{X}}{\mathbf{d}}_{i}$
(14)

where
${\mathbf{X}}=[{\mathbf{x}}_{1},{\mathbf{x}}_{2},\cdots,{\mathbf{x}}_{N}]$,
${\mathbf{d}}_{i}=[v_{1}\dot{\phi}_{i}(t_{1}),v_{2}\dot{\phi}_{i}(t_{2}),\cdots,v_{N}\dot{\phi}_{i}(t_{N})]^{T}$,
and the right hand side is

$\int_{a}^{b}\phi_{i}{\mathbf{f}}({\mathbf{x}},{\mathbf{u}};{\boldsymbol{\Theta}})dt\approx-\sum_{k=1}^{K}v_{k}\phi_{i}(t_{k}){\mathbf{f}}_{k}\equiv{\mathbf{F}}{\mathbf{c}}_{i}$
(15)

where
${\mathbf{f}}_{k}={\mathbf{f}}({\mathbf{x}}_{k},{\mathbf{u}}_{k};{\boldsymbol{\Theta}})$,
${\mathbf{F}}=[{\mathbf{f}}_{1},{\mathbf{f}}_{2},\cdots,{\mathbf{f}}_{N}]$ and
${\mathbf{c}}_{i}=[v_{1}{\phi}_{i}(t_{1}),v_{2}{\phi}_{i}(t_{2}),\cdots,v_{N}{\phi}_{i}(t_{N})]^{T}$.

The total of $n_{\phi}$ integral constraints can be written compactly as

${\mathbf{X}}{\mathbf{D}}={\mathbf{F}}({\mathbf{X}},{\mathbf{U}};{\boldsymbol{\Theta}}){\mathbf{C}}$
(16)

where
${\mathbf{D}}=[{\mathbf{d}}_{1},{\mathbf{d}}_{2},\cdots,{\mathbf{d}}_{n_{\phi}}]$,
${\mathbf{C}}=[{\mathbf{c}}_{1},{\mathbf{c}}_{2},\cdots,{\mathbf{c}}_{n_{\phi}}]$,
and the dependence of ${\mathbf{F}}$ on states, inputs, and parameters is
explicitly written.

The discretized optimization problem with algebraic equality constraints takes
the form:

$\displaystyle\min_{\boldsymbol{\Theta}}$ $\displaystyle\
\left\lVert{\mathbf{Y}}-\hat{\mathbf{Y}}\right\rVert_{v}^{2}$ (17a) s.t.
$\displaystyle\
{\mathbf{X}}{\mathbf{D}}={\mathbf{F}}({\mathbf{X}},{\mathbf{U}};{\boldsymbol{\Theta}}){\mathbf{C}}$
(17b)

To solve the constrained optimization problem in Equation (17), we introduce a
penalty for constraint violations, resulting in the following loss function:

${\mathcal{L}}({\boldsymbol{\Theta}})=\left\lVert{\mathbf{Y}}-\hat{\mathbf{Y}}\right\rVert_{v}^{2}+\gamma\left\lVert{\mathbf{X}}{\mathbf{D}}-{\mathbf{F}}({\mathbf{X}},{\mathbf{U}};{\boldsymbol{\Theta}}){\mathbf{C}}\right\rVert^{2},$
(18)

where $\gamma$ is the penalty parameter. More advanced methods such as the
interior point method may be used for more accurate enforcement of the
constraints.

####  3.1.3 Choice of Test Functions

The choice of test functions $\phi_{i}(t)$ is crucial to the performance of
the weak form. In Ref. [47] the test functions are from the space of piece-
wise polynomials that are non-negative, unimodal, and compactly supported in
$[0,T]$ with a certain number of continuous derivatives depending on the
polynomial degrees.

In this work, we generalize the approach and propose using test functions
based on Jacobi polynomials $P_{n}^{(\alpha,\beta)}(\tau)$, which are
orthogonal on the interval $\tau\in[-1,1]$ with respect to the weight function
$\rho^{(\alpha,\beta)}(\tau)=(1-\tau)^{\alpha}(1+\tau)^{\beta}$. The
orthogonality property of Jacobi polynomials is given by:

$\int_{-1}^{1}P_{n}^{(\alpha,\beta)}(\tau)P_{m}^{(\alpha,\beta)}(\tau)\rho^{(\alpha,\beta)}(\tau)d\tau=K_{n}^{(\alpha,\beta)}\delta_{nm}$
(19)

where $K_{n}^{(\alpha,\beta)}$ is a normalization constant and $\delta_{nm}$
is the Dirac delta. We define the test functions $\phi_{i}(t)$ using Jacobi
polynomials and the weight function $\rho^{(1,1)}(\tau)$ as follows:

$\phi_{i}(t)=P_{i}^{(1,1)}(\tau)\rho^{(1,1)}(\tau),\quad\tau=\frac{2t-(b+a)}{b-a}$
(20)

The variable transformation $\tau$ maps the time interval $[a,b]$ to the
standard interval $[-1,1]$. The weight function
$\rho^{(1,1)}(\tau)=(1-\tau^{2})$ ensures that the boundary conditions in Eq.
(5) are satisfied, as $\phi_{i}(a_{i})=\phi_{i}(b_{i})=0$. The order of the
Jacobi polynomial also becomes a tunable hyperparameter of the weak form.

The use of Jacobi polynomials as test functions offers several potential
advantages over the original choice in Ref. [47]. First, the orthogonality
property of Jacobi polynomials can lead to better convergence in representing
the system trajectories as compared to non-orthogonal basis functions. Second,
Jacobi polynomials are smooth functions, which can help capture the underlying
dynamics more accurately. Third, the weight function $\rho^{(1,1)}(\tau)$
automatically enforces the boundary conditions in Eq. (5), simplifying the
implementation of the weak form. Furthermore, the choice of test functions
based on Jacobi polynomials allows for the accurate approximation of the
integral constraints using quadrature rules. As shown in Ref. [47], the
trapezoidal rule can discretize the weak derivative relation in Eq. (8) to
order $p+1$ accuracy when the test function $\phi$ has roots
$\phi(a)=\phi(b)=0$ of multiplicity $p$ at the endpoints of its support. This
leads to more accurate representation of the dynamics and potentially a better
model.

The choice of the time horizon $N$, polynomial order (PO), and integration
order (IO) in the weak form presents a delicate balance between approximation
accuracy and computational efficiency. The time horizon $N$ determines the
length of the trajectory segment considered by each constraint; an intuition
is that a longer segment captures more dynamics within the time interval and
results in a more accurate characterization of dynamics by the constraint. The
PO affects the expressiveness of the test functions $\phi_{i}(t)$ used in the
weak form, with higher-order polynomials providing more flexibility, and hence
potentially more accuracy, in representing complex dynamics. The IO, which
corresponds to the accuracy of the numerical integration scheme, can be
improved using higher-order methods, such as Simpson’s rule or higher-order
Newton-Cotes formulas. Balancing these hyperparameters is essential for
accurately and efficiently learning dynamical systems using the weak form. In
particular, while larger values of $N$ and PO may result in a more accurate
dynamic model representation, but the increased $N$ and PO would increase the
dimensions of weight matrices ${\mathbf{C}}$ and ${\mathbf{D}}$ in the
integral constraints, Eq. (17b), and thus increase the computational cost in
the evaluation of these constraints during the model training. In general,
cross-validation or other model selection techniques can be used to tune these
hyperparameters for the best trade-off between accuracy and computational
efficiency. The effects of $N$, PO, and IO on the weak form’s performance will
be investigated through numerical examples in Sec. 5.1, providing empirical
evidence of their impact on model accuracy and guidance in the choice of these
parameters.

###  3.2 Weak Latent Dynamics Model

Building upon the weak form approach introduced in the previous section, we
propose the weak Latent Dynamics Model (wLDM), a deep learning-based general
architecture that learns the dynamics of a system in a latent space. The wLDM
leverages the advantages of the weak form to provide an efficient and stable
framework for learning complex, nonlinear dynamics from data. The wLDM
consists of three main components: An encoder, a processor, and a decoder, as
illustrated in Fig. 1. These components work together to map the observed data
into a latent space, learn the dynamics in this latent space, and reconstruct
the original system states from the latent representations. The components of
the wLDM are detailed as follows:

  1. 1.

Encoder: The encoder ${\mathbf{f}}_{E}$ maps the observations ${\mathbf{y}}$
and control inputs ${\mathbf{u}}$ into a $S$-dimensional latent space
representation ${\mathbf{w}}\in{\mathbb{R}}^{S}$ as:

${\mathbf{w}}={\mathbf{f}}_{E}({\mathbf{y}},{\mathbf{u}};{\boldsymbol{\Theta}}_{E})$
(21)

where ${\boldsymbol{\Theta}}_{E}$ denotes the learnable parameters of the
encoder. This latent representation is intended to capture the essential
features and dynamics of the system in a compact and computationally amenable
form. The choice of the encoder architecture depends on the nature of the data
and the problem at hand. For instance, fully connected neural networks (FCNNs)
can be used for simple, structured data, while graph neural networks (GNNs)
may be more suitable for data with complex, graph-like structures. The encoder
consists of a neural network of $N_{E}$ layers.

  2. 2.

Processor: The processor ${\mathbf{f}}_{P}$ is a neural network of $N_{P}$
layers that learns the dynamics in the latent space and produces the time
derivative of the latent states, $\dot{{\mathbf{w}}}\in{\mathbb{R}}^{S}$ as

$\dot{{\mathbf{w}}}={\mathbf{f}}_{P}({\mathbf{w}},{\mathbf{u}};{\boldsymbol{\Theta}}_{P})$
(22)

where ${\boldsymbol{\Theta}}_{P}$ represents the learnable parameters of the
processor. The processor plays a crucial role in the wLDM, as it learns the
underlying dynamics of the system in the latent space. By operating in a
possibly lower-dimensional space, the processor can learn more efficient and
stable representations of the dynamics as compared to learning directly in the
original state space.

  3. 3.

Decoder: The decoder ${\mathbf{f}}_{D}$ is a neural network of $N_{D}$ layers,
which reconstructs the observations from the latent representation
${\mathbf{w}}$,

${\mathbf{y}}={\mathbf{f}}_{D}({\mathbf{w}};{\boldsymbol{\Theta}}_{D})$ (23)

where ${\boldsymbol{\Theta}}_{D}$ represents the learnable parameters of the
decoder.

Figure 1: Structural diagram of the wLDM architecture.

Each component of the wLDM is typically composed of multiple layers of neural
networks, with nonlinear activation functions applied between layers. A common
choice for the activation function is the Parametric Rectified Linear Unit
(PReLU), which allows for learning the optimal slope of the negative part of
the activation function during training. The choice of latent space
dimensionality $S$ is a crucial hyperparameter that balances the trade-off
between model complexity and information loss. A smaller $S$ leads to a more
compact representation but may result in a loss of information, while a larger
$S$ can capture more complex dynamics but may increase computational cost and
the risk of overfitting. In practice, the optimal value of $S$ can be
determined through cross-validation or based on prior knowledge of the
system’s intrinsic dimensionality. The wLDM is trained by minimizing the loss
function in Eq. (18) using standard optimization algorithms, such as Adam,
with an exponential decay scheduling of the learning rate to ensure stable
convergence. Another set of critical hyperparameters in the wLDM are the
choice of intervals $[a_{i},b_{i}]$, i.e., their lengths and distribution.
Empirically, the intervals are chosen to share the same length, characterized
by the number of time steps within the interval $N$; the intervals are
distributed evenly and closely so that neighboring intervals overlap by at
least $N/4$.

The trained model can then make predictions of system dynamics given initial
conditions and control inputs, as summarized in Algorithm 1. Note that the
algorithm leverages the decoding of the predicted observations at each time
step. By encoding the predicted observations into the latent space using the
encoder ${\mathbf{f}}_{E}$ and then decoding the updated latent states back to
the original space using the decoder ${\mathbf{f}}_{D}$, the algorithm ensures
that the predicted trajectories remain within the valid and meaningful region
of the state space captured by the encoder. This decoder step helps to
maintain the stability of the predictions, reduces the influence of noise, and
improves the generalization capability of the model. The benefits of this
approach have been empirically observed, with the decoder-based predictions
consistently leading to more accurate results as compared to evolving the
dynamics solely in the latent space without intermediate encoding and decoding
steps.

Algorithm 1 Predict System Dynamics Using Trained Model

0: Trained wMLP with learned parameters ${\boldsymbol{\Theta}}$

0: Initial observation $\hat{\mathbf{y}}_{0}$, control inputs
${\mathbf{u}}(t)$, prediction horizon $T$

1: Encode initial observation:
${\mathbf{w}}_{0}={\mathbf{f}}_{E}(\hat{\mathbf{y}}_{0})$

2: Initialize predicted trajectory: ${\mathbf{y}}(0)=\hat{\mathbf{y}}_{0}$

3: Initialize time step: $t=0$

4: while $t<T$ do

5: Encode current observation:
${\mathbf{w}}_{t}={\mathbf{f}}_{E}({\mathbf{y}}_{t},{\mathbf{u}}_{t})$

6: Compute latent state derivative:
$\dot{{\mathbf{w}}}_{t}={\mathbf{f}}_{P}({\mathbf{w}}_{t},{\mathbf{u}}_{t};{\boldsymbol{\Theta}}_{P})$

7: Update latent state using ODE solver: ${\mathbf{w}}_{t+\Delta
t}=\text{ODESolver}({\mathbf{w}}_{t},\dot{{\mathbf{w}}}_{t},\Delta t)$

8: Decode latent state: ${\mathbf{y}}_{t+\Delta
t}={\mathbf{f}}_{D}({\mathbf{w}}_{t+\Delta t})$

9: Update time step: $t=t+\Delta t$

10: end while

11: return Predicted trajectory ${\mathbf{y}}(t)$

The key advantage offered by wLDM over conventional methods such as NODEs is
that the former does not rely on the solution of adjoint ODE for gradient
computation during the training process, which can be computationally
expensive and memory-intensive, especially for long sequences, high-
dimensional systems, or stiff systems. Using the weak form, the wLDM can be
trained with significantly less computational cost, which will be demonstrated
in the results section. In addition, the wLDM architecture, as illustrated in
Fig. 1, is flexible in accommodating different types of neural network layers
for the encoder, processor, and decoder, depending on the specific problem and
the structure of the data. For instance, FCNNs can be used for simple,
structured data, while GNNs may be more suitable for data with complex, graph-
like structures.

##  4 Control-Oriented Modeling of Networked Systems

In this section, we introduce the weak Graph Koopman Bilinear Form (wGKBF),
which is a specific adaptation of the wLDM framework tailored to networked
systems. We start by stating the model problem and the corresponding KBF model
for networked systems, subsequently present GNNs as the key component for
incorporating graph topology into the wLDM, and lastly present the full wGKBF
model under the wLDM framework.

###  4.1 Model Problem

In this study, we consider a networked system whose components are sparsely
connected and each may have different numbers of states and inputs. While this
system is physically heterogeneous, a homogeneous graph will be used to
characterize the topology, i.e., the nodes and the edges. The heterogeneous
information on the nodes, e.g., the states and inputs, will be homogenized
later in the wGKBF model.

First, the notations for graphs are briefly introduced. Let
${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ be a graph with a set of $n_{N}$
nodes, ${\mathcal{V}}\subset{\mathbb{N}}^{+}$, and a set of edges
${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$. The graph may
contain directed edges between nodes $I$ and $J$, i.e.,
$(I,J)\in{\mathcal{E}}$ and $(J,I)\notin{\mathcal{E}}$. By convention, $(I,J)$
means $I$ pointing to $J$. Note that capitalized indices are used to
differentiate from indices for time steps. For node $I$, define its neighbors
as ${\mathcal{N}}_{I}=\\{K|(K,I)\in{\mathcal{E}}\\}$. The number of neighbors,
$\left\lvert{\mathcal{N}}_{I}\right\rvert$, is referred to as the degree of
node $I$, $d_{I}$. As will be used later, the connectivity of a graph is
characterized by the adjacency matrix
${\mathbf{S}}\in{\mathbb{R}}^{n_{N}\times n_{N}}$, with ${\mathbf{S}}_{IJ}=1$
if $(I,J)\in{\mathcal{E}}$ and 0 otherwise. Furthermore, the (normalized)
Laplacian matrix is defined as
${\mathbf{L}}={\mathbf{I}}-{\mathbf{D}}^{-1/2}{\mathbf{S}}{\mathbf{D}}^{-1/2}$,
where ${\mathbf{D}}$ is the diagonal matrix with ${\mathbf{D}}_{II}=d_{I}$.

Next, the dynamics of a networked system on graph ${\mathcal{G}}$ are defined.
Suppose node $I$ has states ${\mathbf{x}}_{I}\in{\mathbb{R}}^{n_{I}}$ and
inputs ${\mathbf{u}}_{I}\in{\mathbb{R}}^{m_{I}}$, and the state and input
dimensions can be different among the nodes. The dynamics at node $I$ are
governed by the following system of ODEs:

$\dot{{\mathbf{x}}}_{I}={\mathbf{f}}_{I}({\mathbf{x}}_{I},{\mathbf{u}}_{I})+\sum_{J\in{\mathcal{N}}_{I}}{\mathbf{f}}_{IJ}({\mathbf{x}}_{I},{\mathbf{x}}_{J},{\mathbf{u}}_{I})$
(24)

where the first term accounts for the self-interaction dynamics, with respect
to ${\mathbf{x}}_{I}$ and ${\mathbf{u}}_{I}$, and the summation accounts for
pair-wise interactions due to the neighbors.

The full system may be denoted compactly in the generic form Eq. (1a), where
the collections of states and inputs of the networked dynamics are,
respectively,

${\mathbf{x}}=[{\mathbf{x}}_{1},{\mathbf{x}}_{2},\cdots,{\mathbf{x}}_{n_{N}}]\in{\mathbb{R}}^{n_{x}},\quad{\mathbf{u}}=[{\mathbf{u}}_{1},{\mathbf{u}}_{2},\cdots,{\mathbf{u}}_{m_{N}}]\in{\mathbb{R}}^{n_{u}}$

with $n_{x}=\sum_{I=1}^{n_{N}}n_{I}$ and $n_{u}=\sum_{I=1}^{n_{N}}m_{I}$.

Koopman theory and its extension to control [65, 20] state that there exists a
pair of nonlinear mappings: (1) a lifting, or effectively an encoder,
${\mathbf{w}}={\mathbf{f}}_{E}({\mathbf{x}},{\mathbf{u}})$ that maps physical
states and inputs to the Koopman states ${\mathbf{w}}$, and (2) an inverse
map, or effectively a decoder, ${\mathbf{x}}={\mathbf{f}}_{D}({\mathbf{w}})$
that recovers physical states from Koopman states. In the Koopman space, the
dynamics evolve in a linear manner,
$\dot{{\mathbf{w}}}={\mathbf{A}}{\mathbf{w}}$, where ${\mathbf{A}}$ is a
system matrix to be identified from data or the governing equations. More
details on the Koopman theory are provided in Appendix A.

In this study, we consider control-affine dynamics, for which the simpler KBF
model [20] may be applied,

$\dot{{\mathbf{w}}}={\mathbf{A}}{\mathbf{w}}+\sum_{k=1}^{n_{u}}{\mathbf{B}}_{k}{\mathbf{w}}u_{k}$
(25)

In the general case where the states and inputs are nonlinearly mixed, the
Koopman model may be obtained via, e.g., the KIC framework [65]; however, such
models might not have convenient structures like the KBF that can be leveraged
for controller design.

To learn the KBF model, two elements need to be determined: (1) The nonlinear
mappings ${\mathbf{f}}_{E}$ and ${\mathbf{f}}_{D}$, and (2) the bilinear
system matrices ${\mathbf{A}}$ and ${\mathbf{B}}_{k}$. Intuitively, the
Koopman states for a node $I$ should depend on only the physical states of
node $I$ and its neighbors, since distant nodes should have negligible impact
on the node $I$. To preserve such sparsity induced by the graph topology, GNNs
will be employed to learn the nonlinear mappings, as will be discussed in Sec.
4.2. Subsequently, noting that the Koopman and physical states are effectively
the states and observations in the general nonlinear dynamics, respectively,
the graph KBF model can be directly learned via the wLDM framework; this will
be discussed in Sec. 4.3.

###  4.2 Graph Neural Networks

The core model structure in GNNs that addresses the graph topology is the
message passing (MP) mechanism, which may consist of multiple consecutive MP
steps [32]. At the $j$th MP step on a graph ${\mathcal{G}}$, a vector of
features ${\mathbf{h}}_{I}^{(j)}\in{\mathbb{R}}^{n_{D}^{(j)}}$ at a node $I$
is updated via the following rules,

$\displaystyle\text{Aggregate:}\quad{\mathbf{m}}_{I}^{(j)}={\boldsymbol{\eta}}_{ag}\left(\\{{\mathbf{h}}_{J}^{(j)}\
|\ J\in{\mathcal{N}}_{I}\\}\right),$ (26a)
$\displaystyle\text{Update:}\quad{\mathbf{h}}_{I}^{(j+1)}={\boldsymbol{\eta}}_{up}\left({\mathbf{h}}_{I}^{(j)},{\mathbf{m}}_{I}^{(j)}\right),$
(26b)

where ${\boldsymbol{\eta}}_{ag}$ and ${\boldsymbol{\eta}}_{up}$ are nonlinear
functions, e.g., neural networks, to be learned from data, and
${\mathbf{m}}_{I}^{(j)}$ denotes the information aggregated from the
neighborhood ${\mathcal{N}}_{I}$. Hence the new vector of features at node
$I$, ${\mathbf{h}}_{I}^{(j+1)}\in{\mathbb{R}}^{n_{D}^{(j+1)}}$, is a function
of the features from its neighbors. In essence, one MP step facilitates the
exchange of information between directly connected nodes. By stacking multiple
MP steps, information can be propagated across wider neighborhoods, allowing
each node’s feature vector to be influenced by its more distant neighbors.

In this study, the MP mechanism is implemented using Graph Convolution Layers
(GCLs) [66]. GCLs generalize convolution from multi-dimensional data arrays to
data on graphs and serve as an effective filter that extracts localized
features from graph data. From the point of view of graph Fourier Transforms,
GCLs perform filtering on the spectrum of the Laplacian matrix ${\mathbf{L}}$
of graph ${\mathcal{G}}$. An efficient implementation of the GCL is the
ChebConv network [66], which performs the MP aggregation and updating over all
nodes simultaneously. Collecting the vectors of features at all nodes as
${\mathbf{H}}^{(j)}=[{\mathbf{h}}_{1}^{(j)},{\mathbf{h}}_{2}^{(j)},\cdots,{\mathbf{h}}_{n_{N}}^{(j)}]^{T}\in{\mathbb{R}}^{n_{N}\times
n_{D}^{(j)}}$, we define the $K$th order ChebConv-based graph convolution as

${\mathbf{H}}^{(j+1)}={\boldsymbol{\sigma}}\left(\sum_{k=0}^{K}T_{k}(\tilde{{\mathbf{L}}}){\mathbf{H}}^{(j)}{\boldsymbol{\Theta}}_{k}^{(j)}\right),$
(27)

where ${\boldsymbol{\sigma}}$ is a nonlinear element-wise activation function,
${\boldsymbol{\Theta}}_{k}^{(j)}\in{\mathbb{R}}^{n_{D}^{(j)}\times
n_{D}^{(j+1)}}$ are learnable parameters, and in the summation
$T_{k}(\tilde{{\mathbf{L}}})$ is the $k$th order Chebyshev polynomial of the
scaled Laplacian $\tilde{{\mathbf{L}}}$. The scaled Laplacian is defined as
$\tilde{{\mathbf{L}}}=(2/\lambda_{\max}){\mathbf{L}}-{\mathbf{I}}$ with
$\lambda_{\max}$ the largest eigenvalue of ${\mathbf{L}}$, so that the range
of eigenvalues of $\tilde{{\mathbf{L}}}$ is $[-1,1]$ and matches with the
domain of Chebyshev polynomials. The matrix $T_{k}(\tilde{{\mathbf{L}}})$
contains up to $K$th power of $\tilde{{\mathbf{L}}}$ and hence one $K$th order
ChebConv in Eq. (27) performs $K$ MP steps.

###  4.3 Weak Graph Koopman Bilinear Form

Building upon the principles of Koopman theory (Sec. 4.1), and the GNN
architecture (Sec. 4.2), we now introduce the weak Graph Koopman Bilinear Form
(wGKBF) model, which incorporates graph representations and KBF to capture the
dynamics of control-affine networked systems. The wGKBF model extends the wLDM
framework by employing GCLs in the encoder and decoder, while the processor
takes the form of a KBF. The details are explained below,

  1. 1.

Encoder: The encoder embeds the physical states at each node into a
$\tilde{S}$-dimensional vector of Koopman states so that the entirety of the
physical states are embedded into a Koopman space of $S=n_{N}\times\tilde{S}$
dimensions. The encoder starts with a simple FCNN step to map the physical
states at each node to a feature vector of dimension $\tilde{S}$. For node
$I$,

${\mathbf{h}}_{I}^{(0)}={\mathbf{f}}_{E,I}({\mathbf{x}}_{I},{\boldsymbol{\Theta}}_{E,I})$
(28)

where ${\boldsymbol{\Theta}}_{E,I}$ denotes learnable parameters. Depending on
the applications, the FCNNs at different nodes may be different due to the
heterogeneity in state dimensions; in this case the FCNN step is to ensure the
uniformity of dimensions in the GCL layers. Subsequently, the initial
collection of features is formed as

${\mathbf{H}}^{(0)}=[{\mathbf{h}}_{1},{\mathbf{h}}_{2},\cdots,{\mathbf{h}}_{n_{N}}]^{T}\in{\mathbb{R}}^{n_{N}\times\tilde{S}}$
(29)

and passed through $N_{E}$ GCL layers to produce
${\mathbf{H}}^{(N_{E})}=[{\mathbf{h}}_{1}^{(N_{E})},{\mathbf{h}}_{2}^{(N_{E})},\cdots,{\mathbf{h}}_{n_{N}}^{(N_{E})}]^{T}$,
which is assembled into the vector of Koopman states,

${\mathbf{w}}=[{\mathbf{h}}_{1}^{(N_{E})T},{\mathbf{h}}_{2}^{(N_{E})T},\cdots,{\mathbf{h}}_{n_{N}}^{(N_{E})T}]^{T}\in{\mathbb{R}}^{S}$
(30)

Combining the FCNN and GCL steps, the encoder is denoted

${\mathbf{w}}={\mathbf{f}}_{E}({\mathbf{x}},{\mathcal{G}};{\boldsymbol{\Theta}}_{E}),$
(31)

where ${\boldsymbol{\Theta}}_{E}$ denotes all the learnable parameters in the
encoder.

  2. 2.

Processor: The processor assumes KBF dynamics from Eq. (46) in the latent
space, where the KBF matrices ${\mathbf{A}}\in{\mathbb{R}}^{S\times S}$ and
${\mathbf{B}}_{i}\in{\mathbb{R}}^{S\times S}$ are trainable together with
other network parameters.

  3. 3.

Decoder: The decoder recovers the physical states at each node from the
Koopman space, which is an inversion of the encoder steps. First, given a new
Koopman state vector ${\mathbf{w}}$, it is reorganized into a collection of
features at all nodes
${\mathbf{H}}^{(0)}\in{\mathbb{R}}^{n_{N}\times\tilde{S}}$, which is
subsequently passed through $N_{D}$ GCL layers, arriving at
${\mathbf{H}}^{(N_{D})}=[{\mathbf{h}}_{1}^{(N_{D})},{\mathbf{h}}_{2}^{(N_{D})},\cdots,{\mathbf{h}}_{n_{N}}^{(N_{D})}]^{T}$.
Lastly, a simple FCNN is applied to map feature vector back to the physical
states at each node. For node $I$,

${\mathbf{x}}_{I}={\mathbf{f}}_{D,I}({\mathbf{h}}_{I}^{(N_{D})},{\boldsymbol{\Theta}}_{D,I})$
(32)

Again, these FCNNs at different nodes may be different due to the
heterogeneity in state dimensions. Combining the GCL and FCNN steps, the
decoder is denoted

${\mathbf{x}}={\mathbf{f}}_{D}({\mathbf{w}},{\mathcal{G}};{\boldsymbol{\Theta}}_{D}),$
(33)

where ${\boldsymbol{\Theta}}_{D}$ denotes all the learnable parameters in the
decoder.

When compared to the general wLDM framework, the encoder in wGKBF does not
treat the inputs. Instead, the inputs are provided directly to the processor
following the KBF formulation. Finally, when the physical states at each node
share the same dimension and physical meaning, the FCNN steps in the encoder
and decoder may be skipped.

The incorporation of graph-based representations and the KBF in the wGKBF
model allows for the effective modeling of networked systems, capturing the
intricate interactions and dependencies among the system components. The GCL
layers enable the model to learn localized features from the graph structure
while the KBF processor provides a global bilinearization of the dynamics in
the latent space. By leveraging the expressive power of GNNs and the linearity
of the Koopman operator, the wGKBF model is well suited for accurate
prediction of the dynamics of complex networked systems, as illustrated
through numerical examples next.

##  5 Results and Discussion

To demonstrate the effectiveness of the proposed approach, we investigate
three example nonlinear systems. The first is a lightly-damped double pendulum
system with dimension 4, which serves to validate the wLDM framework and study
the effects of varying weak form hyperparameters through an ablation study.

The second example addresses the Brusselator dynamics, a classic model problem
that can become stiff depending on the choice of parameters and poses
significant challenges to numerical solvers. By applying the wLDM to this
system and comparing to NODEs, we showcase the proposed framework’s capability
to tackle stiff systems of ODEs, demonstrating its robustness and efficiency
in handling such challenging scenarios.

In the third example, we focus on an electrified aircraft system and employ
the wGKBF model with control inputs. This example showcases the model’s
capability in a realistic, multidisciplinary graph-based system, where the
interactions between various components are crucial to capture for accurate
modeling and control. The wGKBF effectively captures the complex dynamics of
the electrified aircraft energy system with a relatively simple bilinear
dynamics model, which may be amenable to the design of accurate yet fast
controllers in this application domain.

###  5.1 Double Pendulum

####  5.1.1 Problem Setup

In the first numerical example, we consider a damped and controlled double
pendulum problem to show the feasibility of the proposed algorithm. The double
pendulum is illustrated in Fig. 2, and described as:

$\displaystyle\ddot{\theta}_{1}$
$\displaystyle=(M_{2}L_{1}\dot{\theta}_{1}^{2}\sin\delta\cos\delta+M_{2}g\sin\theta_{2}\cos\delta+M_{2}L_{2}\dot{\theta}_{2}^{2}\sin\delta-\bar{M}g\sin\theta_{1}+u_{1})/(L_{1}\rho)-\dot{\theta}_{1}$
(34) $\displaystyle\ddot{\theta}_{2}$
$\displaystyle=(-M_{2}L_{2}\dot{\theta}_{2}^{2}\sin\delta\cos\delta+\bar{M}g\sin\theta_{1}\cos\delta-\bar{M}L_{1}\dot{\theta}_{1}^{2}\sin\delta-\bar{M}g\sin\theta_{1}+u_{2})/(L_{2}\rho)-\dot{\theta}_{2}$
(35)

where $\delta=\theta_{2}-\theta_{1}$, $\bar{M}=M_{1}+M_{2}$, and
$\rho=\bar{M}-M_{2}\cos^{2}\delta$. The masses and lengths of the two
pendulums are $M_{1}=1\text{ kg}$ (upper), $M_{2}=1\text{ kg}$ (lower),
$L_{1}=1\text{ m}$, and $L_{2}=1\text{ m}$. Damping terms are added to both
pendulums so as to create a stable isolated equilibrium point in the system.
To generate the dataset for training and testing, we simulate the double
pendulum system with initial conditions randomly selected for
$\theta_{i}\in[-10^{\circ},10^{\circ}]$ and
$\dot{\theta}_{i}\in[-10^{\circ}/s,10^{\circ}/s]$ and control
$u_{i}\in[-0.25,0.25]\textrm{ N}$. In this example, the control is held
constant for each trajectory to reduce the variability in the subsequent error
analysis. The more practical case of time-varying control is considered in the
third example. A total of 320 trajectories are generated for model training,
followed by another 100 trajectories for testing. Each trajectory is simulated
with a 4th order Runge-Kutta method for 20s with a step size of 0.01 seconds.
All trajectories were normalized to $[0,1]$.

Figure 2: Illustration of the double pendulum.

To evaluate the effectiveness of the proposed weak form, we compare the
performance of the wLDM to two baseline methods: Long Short-Term Memory (LSTM)
networks [67, 68] and Neural Ordinary Differential Equations (NODEs) [53]. The
LSTM network is a type of Recurrrent Neural Network, with its ability to
capture long-term dependencies and sequences making it particularly suited for
sequential data analysis. As a discrete-time model, LSTM networks predict
system observations at fixed time steps. In this case, a standard LSTM network
implementation processes $L$ steps of observation and control
$({\mathbf{Y}},{\mathbf{U}})$ through two LSTM layers, followed by a linear
output layer that predicts the observation at the next step
$\tilde{{\mathbf{y}}}$. On the other hand, NODEs are chosen as a continuous-
time baseline model. Including NODEs as a baseline allows for a direct
comparison between the proposed wLDM and a state-of-the-art continuous-time
approach. The NODE is implemented with multiple layers of FCNN and PReLU
activation functions, which integrates a predictive trajectory based on the
initial condition and control signals and is then trained with the procedures
outlined in Sec. 2.2.2.

The hyperparameters of the nominal wLDM, with a total of 1,481 trainable
parameters, are outlined in Table 5 and follow the architecture specified in
Sec. 3.2. The nominal model serves as the basis of comparison for the
subsequent benchmark against different methods and wLDM’s having different
parameters. To ensure a fair comparison, both baseline methods are implemented
to match the size of the nominal wLDM. The LSTM comprises a single LSTM layer
and an output layer of FCNN (both with dimension $S=32$) followed by PReLU
activation, totaling 5,191 parameters. The NODE shares the same architecture
as the nominal wLDM, as outlined in Table 5, resulting in an identical
parameter count.

Throughout this and the following examples, all models are implemented using
the open-source machine learning framework, PyTorch [69]. The training and
assessment of all the models are performed using a single NVIDIA RTX3090 GPU.

Table 5: Detailed hyperparameters of the nominal wLDM. Component | Layer Type | Parameters
---|---|---
Encoder | FCNN |  $N_{E}=1$, $S=32$
| Activation | PReLU
Processor | FCNN |  $N_{P}=1$, $S=32$
| Activation | PReLU
Decoder | FCNN |  $N_{D}=1$, $S=32$
| Activation | PReLU

####  5.1.2 Numerical Results

To compare the model predictions, we train each model on the dataset generated
from the double pendulum system with the specified range of parameters. The
trained models are then used to predict the future states of the system given
initial conditions and control inputs that are unseen in the training dataset,
as outlined in Sec. 3.2. The predictive accuracy of the models is quantified
by Normalized Root Mean Square Error (NRMSE). For a trajectory of $L$ time
steps, the NRMSE is defined as,

$\text{NRMSE}=\frac{1}{L}\sum_{l=1}^{L}\frac{\sqrt{\frac{1}{n_{y}}\sum_{i=1}^{n_{y}}(\tilde{y}_{l,i}-y_{l,i})^{2}}}{\max_{i}(y_{l,i})-\min_{i}(y_{l,i})},$
(36)

where $[y_{l,1},\cdots,y_{l,n_{y}}]$ denotes the $n_{y}$ elements of the
observations at the $l$th step of the true trajectory.

Figure 3: Comparison of prediction performance of the wLDM and the two
baseline methods for the damped double pendulum dynamics with controls.

Figure 3 compares the model performance of the nominal wLDM against the
baseline models for a representative test case. The predictions from the wLDM
and NODE closely match the ground truth over the entire prediction horizon,
capturing both the transient dynamics and equilibrium points accurately. There
is virtually no deviation in the prediction of the wLDM and NODE, with the
wLDM achieving the lowest NRMSE of $1.15\times 10^{-3}$. Predictions from the
LSTM network show notable deviation across all 4 states in both transient and
steady-state behavior, leading to a higher NRMSE of $2.64\times 10^{-2}$.
These deviations, particularly pronounced in long-term predictions, are likely
exacerbated by cumulative errors. A comparison of the predictive NRMSE for all
three models across the entire test dataset, and training time, is presented
in Table 6. The wLDM achieves the best results in both mean and standard
deviation of NRMSE, demonstrating a modest improvement over the NODE and a
substantial advantage over the LSTM network. The training time of the wLDM is
significantly lower than the baseline methods. This showcases the exceptional
accuracy and training efficiency of the wLDM approach.

Table 6: Summary of model performance for the double pendulum. Model | NRMSE | Training
---|---|---
Mean | Std. | Min | Max | Time (s)
wLDM | 9.41E-4 | 5.05E-4 | 1.53E-4 | 2.87E-3 | 2940
NODE | 1.11E-3 | 8.03E-4 | 1.08E-4 | 4.98E-3 | 25480
LSTM | 2.28E-2 | 1.26E-3 | 1.96E-2 | 2.64E-2 | 64150

####  5.1.3 Parametric Study

To delve deeper into the characteristics of training with a weak form
approach, we performed parametric studies with the wLDM and the double
pendulum dynamics, varying the time horizon $N$, the order of Jacobi
polynomials PO, and integration order IO. To ensure the fair comparison among
the models of different combinations of $N$, PO and IO, the latent space
dimension $S$ is increased from 32 to 128. This is because as PO increases,
the wLDM may learn more complex dynamical behaviors from data and would
require a model having more parameters to characterize such increased
complexity.

Figure 4: Comparison of predictive performance for models trained with
increasing time window horizon $N$.

As a first parametric study, we fixed the PO and IO to isolate and understand
the effect of incrementally increasing the time horizon $N$, i.e., the number
of steps in the time interval of a test function, on the model’s predictive
performance. The results in Fig. 4 demonstrate the impact of the time horizon
$N$ on the predictive performance of the wLDM. As $N$ increases, we observe a
consistent decrease in the predictive NRMSE, supporting the intuition
discussed in Sec. 3.1.3 that using a longer trajectory segment makes the
integral constraints more accurate representations of dynamics and leading to
a more accurate model. This improvement in performance is particularly
noticeable when $N$ is increased from smaller values (e.g., from 3 to 11),
indicating that even a modest increase in the time horizon can lead to
significant gains in accuracy. However, the gains in predictive performance
become less significant as $N$ grows larger, evident from the diminishing
returns in NRMSE reduction beyond $N=27$ in both plots, suggesting a trade-off
between model size and computational efficiency. Furthermore, comparing the
left plot ($\text{PO}=\text{IO}=1$) and the right plot
($\text{PO}=\text{IO}=2$), we observe that the model with higher PO and IO
achieves lower NRMSE values for the same $N$. This suggests that higher-order
polynomials and integration schemes can help capture the system dynamics more
effectively, potentially reducing the required $N$ for a given level of
performance. The findings concerning the time interval $N$ are empirical and
worth more in-depth theoretical analysis, which can be a valuable future
research direction.

Figure 5: Trained wLDM performance with varying PO and IO at a fixed $N=61$
with $95\%$ confidence intervals computed using a bootstrap method.

The second parametric study investigates the effects of increasing the PO and
IO on the model’s predictive performance. To allow for the testing of higher
integration orders, which require more time steps, the time horizon $N$ is
chosen to be 61. The results, shown in Fig. 5, indicate that higher-order
polynomials in the weak form offer more flexibility in representing complex
dynamics, and higher-order integration schemes provide better approximations
of the integral constraints. As PO increases from 1 to 6, the mean predictive
NRMSE decreases consistently across all integration orders, suggesting that
with a sufficient time horizon $N$, increasing the expressiveness of the test
functions through higher-order polynomials leads to more accurate models.

Moreover, for a given PO, increasing the IO generally results in lower NRMSE
values and tighter confidence intervals, indicating that higher-order
integration schemes lead to more accurate and confident predictions by better
approximating the integral constraints. The confidence intervals are computed
using a bootstrap method, where a subset of the RMSE values is randomly
sampled with replacement, and the mean RMSE is calculated for each bootstrap
sample. The tighter confidence intervals associated with higher IO values
demonstrate that the model’s predictions are more consistent and reliable when
using higher-order integration schemes. However, the gains in performance
become less pronounced at higher orders, suggesting a trade-off between the
increased computational cost in training (due to more constraints associated
with higher-order polynomials) and the resulting performance of the model. The
optimal choice of PO and IO should balance the desired level of accuracy with
the computational cost associated with higher orders.

###  5.2 Brusselator Dynamics

####  5.2.1 Problem Setup

In this numerical example, we investigate the Brusselator system, a classic
model that originates from chemical kinetics and consists of a pair of
nonlinear ODEs [60]. The governing equations of the Brusselator model are:

$\displaystyle\dot{x}_{1}$ $\displaystyle=A+x_{1}^{2}x_{2}-(B+1)x_{1}$ (37)
$\displaystyle\dot{x}_{2}$ $\displaystyle=Bx_{1}-x_{1}^{2}x_{2}$ (38)

where $A$ and $B$ are constants that determine the stiffness of the system.
This model has tunable numerical stiffness so it is suitable for assessing the
capabilities of different methods to learn stiff systems.

To assess the predictive accuracy and training efficiency of the wLDM and NODE
in handling stiff systems, we generate training data with $A=1$ while varying
$B$ between 2 and 5. The increasing stiffness of the system is quantified by
computing the stiffness ratio $S_{R}$, defined as the ratio of the maximum to
the minimum absolute eigenvalues of the Jacobian matrix:

$S_{R}=\frac{\max(|\lambda_{i}|)}{\min(|\lambda_{i}|)}$ (39)

where $\lambda_{i}$ are the eigenvalues of the Jacobian matrix of the
Brusselator system, given by:

$J=\begin{bmatrix}2x_{1}x_{2}-(B+1)&x_{1}^{2}\\\
B-2x_{1}x_{2}&-x_{1}^{2}\end{bmatrix}$ (40)

The stiffness ratio is computed for each trajectory, confirming the increase
in stiffness as $B$ increases, with a minimum of 4.33 at $B=2$ and a maximum
of 261.21 at $B=5$.

Initial conditions for the trajectories were randomly chosen from within the
range $[0,2]$ for both $x_{1}$ and $x_{2}$. Each trajectory was evolved over a
time period of $t\in[0,20]$ seconds using the backward differentiation formula
with a timestep of $\Delta t=0.2$ s. To adequately train the models, 200
trajectories were created for each value of $B$, complemented by an additional
set of 100 trajectories for model testing.

The nominal wLDM used in this example has a latent space dimension of $S=128$
and $N_{P}=2$ processor layers to ensure sufficient capacity for handling the
more challenging stiff dynamics. The model was trained with increasing time
horizons, from $N=61$ to $N=241$, to capture both the long-term dynamics and
short-term spikes accurately. This approach proved particularly useful in the
prediction of highly stiff cases, such as $B=5$; in the case of $B=5$, the
system exhibits typical phenomenon in stiff dynamics, where the response is
smooth most of the time but sudden large variations appear within a short time
window. The PO and IO are both set to be 4.

The baseline NODE model was implemented with a similar architecture, matching
the number of layers and hidden dimensions to ensure a fair comparison. The
models were trained until the relative change in the loss function fell below
a tolerance of $10^{-4}$ or a maximum of 10,000 iterations was reached.

####  5.2.2 Numerical Results

(a) $B=2$

(b) $B=3$

(c) $B=4$

(d) $B=5$

Figure 6: Predictive comparison of the wLDM and baseline NODE model on
Brusselator dynamics with increasing stiffness.

Figure 6 compares the predictive accuracy of the wLDM and the baseline NODE
model on Brusselator dynamics with increasing stiffness. As $B$ increases, the
system undergoes a transition from decaying oscillations converging to a limit
cycle ($B=2$) to highly stiff dynamics with sharp, almost discontinuous
transitions ($B=5$). These abrupt changes in the system’s behavior pose
significant challenges to the models in terms of both training and prediction.

The comparison of the NRMSE of the wLDM and NODE are shown in Table 7. The
NRMSE values highlight the wLDM’s consistently better performance as compared
to NODE across all stiffness levels, and the performance gap becomes more
pronounced as $B$ increases. For $B=4$, NODE struggles to capture the sharp
peaks accurately, resulting in a higher NRMSE. In the highly stiff case of
$B=5$, NODE fails to capture the stiffness of the state $x_{1}$, leading to a
substantial increase in the NRMSE.

Table 7: Summary of model performance for the Brusselator dynamics. Model | $B=2$ | $B=3$ | $B=4$ | $B=5$
---|---|---|---|---
wLDM | 0.003 | 0.001 | 0.011 | 0.013
NODE | 0.011 | 0.012 | 0.041 | 0.159

The phase portrait analysis in Fig. 6 provides valuable insights into the
models’ ability to capture the qualitative behavior of the Brusselator system.
In the case of $B=2$, both the wLDM and NODE successfully reproduce the spiral
convergence towards a stable limit cycle, indicating their capability to learn
the overall dynamics in a less stiff scenario. At $B=3$, the wLDM can
accurately capture the elongated and asymmetric shape of the limit cycle,
closely matching the ground truth. In contrast, NODE shows deviations,
particularly in the outer regions of the limit cycle.

In the highly stiff scenarios ($B=4$ and $B=5$), the advantages of the wLDM
become even more pronounced. This model proves more adept at reproducing the
complex loops, spirals, and scattered trajectories characteristic of the
Brusselator system in these regimes. This can be attributed to the use of
higher-order polynomials and integration schemes in the weak form, which
provide more flexibility in representing stiff dynamics with rapid
transitions. On the other hand, NODE struggles to match the ground truth in
these highly stiff cases, especially in the regions of stiff dynamics with
rapid transitions. This is likely due to the challenges associated with the
numerical integration of stiff ODEs, which can lead to instability and
inaccuracies in the learned dynamics.

Figure 7: Comparison of training efficiency.

Figure 7 compares the training efficiency of the wLDM and NODE models by
plotting the training time per iteration against the number of iterations
required for convergence. The results demonstrate that the wLDM requires
significantly less computational resources, with each iteration taking less
than 1s and convergence achieved in around 4500 iterations. In contrast, NODE
models require longer training times per iteration and more iterations to
converge, resulting in a much larger total training time as compared to the
wLDM. Particularly, the training time per iteration of the NODE models
increases as the stiffness increases, but that of the wLDM is insensitive to
the stiffness. While the NODE models require more iterations to converge with
increased stiffness, there is an exception of $B=5$, noted specially using a
hollow marker in Fig. 7. In this case, early convergence is achieved but with
poor prediction performance; this phenomenon is also observed in other studies
[60].

Overall, this example demonstrates the superior performance of the wLDM as
compared to the baseline NODE model in learning stiff ODEs, consistently
achieving better predictive accuracy, capturing intricate details and complex
behaviors, and faithfully representing qualitative features in the phase
portrait analysis. Moreover, the wLDM exhibits significantly faster training
times and converges within a consistent range of iterations, making it a more
computationally efficient and practical choice for modeling stiff ODEs.

###  5.3 Electrified Aircraft System

####  5.3.1 Problem Setup

Finally, for a more realistic and complex example, we consider the electrified
aircraft system shown in Fig. 8. This electro-thermo-mechanical system
consists of a battery delivering power to an inverter and motor, a simplified
1-D model of the longitudinal vehicle dynamics, and a single-phase thermal
management system responsible for cooling the motor, consisting of a cold
plate, heat exchanger, and tank [44].

Figure 8: Schematic of the electrified aircraft system. Figure 9: Graph-based
model structure of the electrified aircraft system.

The system dynamics are modeled using the modular graph-based modeling
approach from [70, 71], producing the oriented graph shown in Fig. 9 where
each node represents a dynamic state or exogenous coupling term and edges
indicate which states directly exchange energy. The full-order system, shown
in Fig. 9, has 24 dynamic states and 3 control inputs, specifically the
inverter duty cycle, primary coolant flow rate driven by the pump, and
secondary flow rate through the other side of the heat exchanger.

In this example, the physics-based model [44] is assumed to be the true model
that generates the training data, and the goal is to demonstrate the
capability of wGKBF model to learn the nonlinear dynamics of the energy
system. A future step is to use experimental measurements of real hardware
system to train the wGKBF model.

We apply the wGKBF model as a method of data-driven control-oriented modeling
that incorporates the graph structure of the system. To train the wGKBF model,
we generate 75 trajectories by simulating the system with randomly perturbed
initial conditions, constrained within the predefined operational boundaries.
An additional 25 trajectories are generated for testing the trained model’s
predictive performance. The system is subjected to chirp control signals with
initial and final frequencies randomly sampled from 1 to 4 Hz, ensuring a rich
and complex excitation of dynamics across the system. The data are again
generated using the backward differentiation formula, with a simulation time
of 2 seconds and $\Delta t=0.001$. The resulting time series data is then
normalized to the range [0, 1] before being used for model training and
testing.

We benchmark the predictive performance of wGKBF against the conventional
implementation of KBF, i.e., the biDMDc approach [19]. The biDMDc model
utilizes past system states, which is a form of time-delay embedding that
incorporates historical data to enhance predictive accuracy. We use 10 steps
of past physical states of the system and a truncation order of 20 for the DMD
modes, determined after a parametric sweep to achieve the best possible
predictive accuracy. To ensure a fair comparison, the wGKBF model also
incorporates the same 10 steps of past physical states as node features in the
graph. Theoretically, the graph-based KBF could be trained without the weak
form, e.g., using the conventional graph NODE methods [54, 55, 56, 57].
However, these methods resulted in intractable computational cost in the
learning of the electrified aircraft dynamics due to the high numerical
stiffness. Hence in this study only the wGKBF and biDMDc results are reported
and compared.

Table 8: Detailed hyperparameters of the nominal wGKBF model. Component | Layer Type | Parameters
---|---|---
Encoder | GCL |  $N_{E}=2$, $n_{N}=24$, $\tilde{S}=16$
| Activation | PReLU
Processor | KBF (${\mathbf{A}},{\mathbf{B}}_{i})$ |  $S=384$, $i=1,2,3$
Decoder | GCL |  $N_{D}=2$, $\tilde{S}=16$
| Activation | PReLU

The hyperparameters of the nominal wGKBF model are summarized in Table 8. In
this case, the states at all nodes share the same dimension of 1, hence the
FCNN steps in the encoder and decoder are removed. The encoder and decoder
both consist of two GCL layers with PReLU activation functions. The latent
space dimension is set to $\tilde{S}=16$, determined through a parametric
sweep to optimize model prediction accuracy while maintaining a tractable
number of parameters in the KBF processor. The number of nodes in the graph,
$n_{N}=24$, matches the number of system states. Consequently, the KBF
processor has a state dimension of $S=384$, with three input matrices
${\mathbf{B}}_{i}$ corresponding to the three control inputs.

Figure 10: Comparison of wGKBF and biDMDc prediction for the electrified
aircraft system.

####  5.3.2 Numerical Results

Figure 10 compares the performance of wGKBF and biDMDc in predicting the
simulated aircraft system dynamics with chirp control signals. The true system
dynamics, represented by the black dashed lines, exhibit a mix of steady,
slowly varying, and rapidly fluctuating behaviors across different state
variables. The thermal states, such as cold plate fluid and wall temperatures
($x_{1}$ and $x_{2}$), show slow response characteristics, while some
electrical states, including inverter voltages and currents ($x_{16}$,
$x_{17}$, $x_{19}$, and $x_{20}$), display fast dynamics with significant
variations. The mechanical states, like motor speed ($x_{22}$) and vehicle
speed ($x_{24}$), demonstrate intermediate dynamics.

The wGKBF model exhibits superior performance in capturing the true dynamics
of the electrified aircraft system, achieving an NRMSE of 0.026, significantly
lower than the NRMSE of biDMDc of 0.475. The wGKBF’s ability to accurately
predict the rapid fluctuations in electrical states and the subtle variations
in thermal states highlights its effectiveness in handling the multi-domain
and multi-timescale dynamics of the electro-thermo-mechanical system.

This enhanced performance can be attributed to several aspects of the wGKBF
model. The graph-based representation allows the wGKBF to explicitly model the
interactions and dependencies among the various components of the aircraft
system. By incorporating the system’s inherent structure into the model
architecture, the wGKBF can more effectively capture localized features and
propagate information across the graph, leading to a more accurate and
interpretable representation of the system dynamics. Furthermore, the deep
learning nature of the wGKBF model endows it with a high level of
expressiveness, allowing it to learn complex nonlinear mappings between the
input and output spaces. The combination of GCLs and KBF in the wGKBF
architecture enables the model to capture the dynamics of multiple timescales
of the system, resulting in improved predictive accuracy compared to the
biDMDc baseline.

Figure 11: Comparison of predictive performance of wGKBF and biDMDc across all test cases. Table 9: Summary of model performance for the electrified aircraft system. Model | NRMSE | Training
---|---|---
Mean | Std. | Min | Max | Time (s)
wGKBF | 3.38E-2 | 1.12E-2 | 6.47E-2 | 1.84E-2 | 16680
biDMDc | 5.82E-1 | 6.23E-2 | 6.97E-1 | 4.73E-1 | 724

The predictive performance of the two models across all 25 test cases is shown
in Fig. 11, and summarized in Table 9. The wGKBF maintains a consistently low
NRMSE with a mean value of $3.38\times 10^{-2}$ and a small standard deviation
of $1.12\times 10^{-2}$. In contrast, the biDMDc model, while takes
significantly less time to train, shows greater variability in prediction
error, with higher NRMSE values and a larger standard deviation. The
substantial gap in performance between the two models highlights the
effectiveness of the wGKBF approach in leveraging the graph topology and deep
learning to capture the complex nonlinear dynamics of the electrified aircraft
system. By incorporating the system’s structure and learning expressive
mappings, the wGKBF model achieves significantly improved predictive accuracy
and consistency, showcasing its superior capability in modeling the complex
electro-thermo-mechanical systems.

##  6 Conclusions

Accurate and computationally efficient control-oriented models are essential
for the design and implementation of advanced control strategies to optimize
system performance, safety, and efficiency of networked systems, such as
electrified aircraft energy systems. The development of such models requires
answering two questions: What model form to use and how to train the model?
This paper addresses the second question by introducing a novel deep learning
approach, the weak Latent Dynamics Model (wLDM), for learning generic
nonlinear dynamics with control. Leveraging the weak form, the wLDM enables
more numerically stable and computationally efficient training as well as more
accurate predictive model, when compared to traditional approaches. Addressing
the first question, we build upon the wLDM and propose the weak Graph Koopman
Bilinear Form (wGKBF) model that integrates graph-based representations and
Koopman theory to learn control-oriented bilinear dynamics for networked
systems.

The effectiveness of the wLDM and its wGKBF variant is demonstrated through
three example systems of increasing complexity:

  1. 1.

A controlled double pendulum: The numerical experiments show that the wLDM
consistently outperforms baseline models in terms of predictive accuracy and
training efficiency. Parametric studies provide valuable insights into the
role of hyperparameters in the weak form. The time horizon $N$ is found to
have a significant impact on predictive performance, with larger values
leading to more accurate models at the cost of increased training time. The
order of Jacobi polynomials (PO) and integration order (IO) also play crucial
roles, with higher orders capturing more complex dynamics but exhibiting
diminishing returns in performance gains.

  2. 2.

The Brusselator dynamics: A systematic study is performed to demonstrate the
capability of the wLDM to learn stiff dynamics, which are challenging for
conventional learning algorithms such as NODEs. It is found that the wLDM
outperforms the NODEs in both predictive accuracy and training efficiency. In
particular, the NODEs take substantially longer time to train as the stiffness
increases while the training time of the wLDM is shown to be insensitive to
the stiffness. This comparison highlights the potential of the wLDM for
learning complex stiff dynamics, such as the systems having disparate
timescales.

  3. 3.

An electrified aircraft energy system: The wLDM is extended to wGKBF by
incorporating graph-based representations and the Koopman bilinear form. The
wGKBF model effectively captures the intricate couplings and dependencies
among the various components of the electrified aircraft energy system and
outperforms the existing bilinear system identification method in predictive
accuracy by an order of magnitude.

Collectively, the wLDM and wGKBF model offer a powerful and flexible framework
for learning the dynamics of complex networked systems in aerospace
applications. By leveraging the weak form, graph-based representations, and
Koopman theory, these models enable accurate, efficient, and control-oriented
modeling of challenging systems, opening up new possibilities for advanced
system modeling and control. Future work will focus on leveraging these models
for the design of advanced control strategies, such as hierarchical MPC and
reinforcement learning-based approaches, as well as system analysis tasks,
including stability assessment, sensitivity studies, and optimization of
system configurations. Additionally, investigating the scalability and
generalizability of the proposed framework to larger-scale systems with
higher-dimensional state spaces and more intricate network structures will be
a key direction for further research.

##  Appendix A Koopman Theory and Bilinear Form

Koopman theory provides a convenient tool for the analysis and control of
nonlinear dynamical systems. When the system is autonomous and has an isolated
stable equilibrium point, a global linearization can be achieved in the entire
attracting basin of the system [15]. With such linearization, classical linear
system theory can be leveraged to characterize the nonlinear dynamics. For a
controlled system, the Koopman operator can be extended to the bilinear form
[20], allowing for the efficient representation and control of nonlinear
systems with inputs. For the more general controlled case, see Ref. [65].

Consider a control-affine system of the form

$\dot{{\mathbf{x}}}={\mathbf{f}}_{0}({\mathbf{x}})+\sum_{i=1}^{n_{u}}{\mathbf{f}}_{i}({\mathbf{x}})u_{i}\>,\quad{\mathbf{x}}(0)={\mathbf{x}}_{0}\>,$
(41)

where ${\mathbf{x}}\in{\mathbb{X}}\subseteq{\mathbb{R}}^{n_{x}}$ and
${\mathbf{u}}=[u_{1}\dots u_{n_{u}}]^{T}\in{\mathbb{R}}^{n_{u}}$ are states
and control, respectively,
${\mathbf{f}}_{0}:{\mathbb{X}}\rightarrow{\mathbb{X}}$ captures the autonomous
dynamics, and ${\mathbf{f}}_{i}:{\mathbb{X}}\rightarrow{\mathbb{X}}$ captures
the influence of control actions on the system.

In the absence of control inputs [17], i.e., when ${\mathbf{u}}=0$, the system
generates a flow ${\mathbf{F}}_{t}({\mathbf{x}}_{0})={\mathbf{x}}(t)$ from an
initial condition ${\mathbf{x}}_{0}$. The continuous time Koopman operator
${\mathcal{K}}_{t}:{\mathcal{F}}\rightarrow{\mathcal{F}}$ is an infinite-
dimensional linear operator such that
${\mathcal{K}}_{t}z=z\circ{\mathbf{F}}_{t}$ for all $z\in{\mathcal{F}}$, where
$z:{\mathbb{X}}\rightarrow{\mathbb{C}}$ is a complex-valued observable
function of the state vector ${\mathbf{x}}$, ${\mathcal{F}}$ is the function
space of all possible observables, and $\circ$ denotes function composition.
As a linear operator, ${\mathcal{K}}_{t}$ admits eigenpairs
$(\lambda,\varphi)$ such that

${\mathcal{K}}_{t}\varphi=\varphi\circ{\mathbf{F}}_{t}=e^{\lambda
t}\varphi\>,$ (42)

where $\lambda\in{\mathbb{C}}$ and $\varphi\in{\mathcal{F}}$ are the Koopman
eigenvalue and Koopman eigenfunction, respectively.

The infinitesimal generator of ${\mathcal{K}}_{t}$ associated with
${\mathbf{f}}_{0}$, referred to as the Koopman generator, is defined as
${\mathcal{L}}_{{\mathbf{f}}_{0}}=\lim_{t\rightarrow
0}\frac{{\mathcal{K}}_{t}-I}{t}$, where $I$ is the identity operator, and
turns out to be the Lie derivative
${\mathcal{L}}_{{\mathbf{f}}_{0}}={\mathbf{f}}_{0}\cdot\nabla$, with eigenpair
$(\lambda,\varphi)$,

$\dot{\varphi}={\mathcal{L}}_{{\mathbf{f}}_{0}}\varphi=\lambda\varphi\>.$ (43)

Given a set of eigenpairs $\\{(\lambda_{i},\varphi_{i})\\}_{i=1}^{n}$, the
Koopman Canonical Transform (KCT) [16] of the control-affine system in Eq.
(41) is given as

$\dot{{\boldsymbol{\varphi}}}={\boldsymbol{\Lambda}}{\boldsymbol{\varphi}}+\sum_{i=1}^{n_{u}}{\mathcal{L}}_{{\mathbf{f}}_{i}}{\boldsymbol{\varphi}}u_{i}\>,$
(44)

where
${\boldsymbol{\Lambda}}=\mathrm{Diag}([\lambda_{1},\cdots,\lambda_{n}])$,
${\boldsymbol{\varphi}}=[\varphi_{1},\cdots,\varphi_{n}]$, and Lie derivatives
for the control terms are
${\mathcal{L}}_{{\mathbf{f}}_{i}}={\mathbf{f}}_{i}\cdot\nabla$.

Suppose the set of eigenfunctions is sufficiently large, such that
${\boldsymbol{\varphi}}$ span an invariant space for
${\mathcal{L}}_{{\mathbf{f}}_{i}}$, i.e., each
${\mathcal{L}}_{{\mathbf{f}}_{i}}$ can be represented using an $n\times n$
matrix ${\mathbf{D}}_{i}$ such that
${\mathcal{L}}_{{\mathbf{f}}_{i}}{\boldsymbol{\varphi}}={\mathbf{D}}_{i}{\boldsymbol{\varphi}}$.
Then the KCT can be brought to a bilinear form [20],

$\dot{{\boldsymbol{\varphi}}}={\boldsymbol{\Lambda}}{\boldsymbol{\varphi}}+\sum_{i=1}^{n_{u}}{\mathbf{D}}_{i}{\boldsymbol{\varphi}}u_{i}\>.$
(45)

Often it is difficult to directly obtain the eigenfunctions of
${\mathcal{K}}_{t}$. Instead, by introducing lifted coordinates via a mapping
${\mathbf{w}}={\boldsymbol{\phi}}({\mathbf{x}})$ such that
${\boldsymbol{\varphi}}={\boldsymbol{\Psi}}^{H}{\mathbf{w}}$, and substituting
it into the KCT in bilinear form from Eq. (45), we reach the commonly used
Koopman Bilinear Form (KBF) [20, 21],

$\dot{{\mathbf{w}}}={\mathbf{A}}{\mathbf{w}}+\sum_{i=1}^{n_{u}}{\mathbf{B}}_{i}{\mathbf{w}}u_{i}\>,$
(46)

where an eigendecomposition
${\mathbf{A}}={\boldsymbol{\Phi}}{\boldsymbol{\Lambda}}{\boldsymbol{\Psi}}^{H}$
reproduces the Koopman eigenvalues ${\boldsymbol{\Lambda}}$ and
${\mathbf{B}}_{i}={\boldsymbol{\Phi}}{\mathbf{D}}_{i}{\boldsymbol{\Psi}}^{H}$.
The original states are recovered from an inverse mapping as
${\mathbf{x}}={\boldsymbol{\phi}}^{-1}({\mathbf{w}})\equiv{\boldsymbol{\psi}}({\mathbf{w}})$.

The KBF provides a global bilinearization of the control-affine system in Eq.
(41), enabling the application of linear and bilinear control techniques to
the nonlinear aircraft energy system. By incorporating the KBF into the wLDM
framework, we aim to learn a control-oriented representation of the aircraft
energy system that captures its nonlinear dynamics while maintaining a
structure amenable to efficient control design.

The condition for bilinearization into the form Eq. (46) is restated as
follows:

###### Theorem 1

[20] Suppose a set of the Koopman eigenfunctions
$\\{\varphi_{1},\varphi_{2},\dots,\varphi_{n}\\}$ of the autonomous system
forms a invariant subspace of ${\mathcal{L}}_{{\mathbf{f}}_{i}},i=1,\dots,m$.
Then Eq. (41) is bilinearizable with an $n$ dimensional state space.

## Acknowledgements

This work was partially supported by a seed grant from The Pennsylvania State
University Institute for Computational and Data Sciences (ICDS) . The research
of YY and DH was partially supported by National Science Foundation under the
award DMS-2229435.

## References

  * Boglietti et al. [2009] Boglietti, A., Cavagnino, A., Tenconi, A., and Vaschetto, S., “The Safety Critical Electric Machines and Drives in the More Electric Aircraft: A Survey,” _Industrial Electronics Conference_ , 2009, pp. 2587–2594.
  * Roberts and Decker [2013] Roberts, R. A., and Decker, D. D., “Energy Optimization of an Aircraft Focused on Component Sizing and Control Architecture Interactions,” _11th International Energy Conversion Engineering Conference_ , 2013, p. 3805.
  * Karimi [2007] Karimi, K. J., “Future Aircraft Power Systems-Integration Challenges,” , 2007.
  * Gerstler and Bunker [2008] Gerstler, W. D., and Bunker, R. S., “Aircraft Engine Thermal Management: The Impact of Aviation Electric Power Demands,” _ASME Global Gas Turbine News_ , 2008.
  * Park and Pangborn [2021] Park, S., and Pangborn, H. C., “Power and Thermal Management with Battery Degradation for Hybrid Electric Vehicles,” _IEEE Conference on Control Technology and Applications_ , 2021, pp. 832–838.
  * Newman [2018] Newman, M., _Networks: An Introduction_ , Oxford University Press, 2018.
  * Slotine and Li [1991] Slotine, J., and Li, W., _Applied Nonlinear Control_ , Prentice-Hall International Editions, Prentice-Hall, 1991.
  * Khalil [2002] Khalil, H. K., _Control of Nonlinear Systems_ , Prentice Hall, New York, NY, 2002.
  * Aksland et al. [2023] Aksland, C. T., Tannous, P. J., Wagenmaker, M. J., Pangborn, H. C., and Alleyne, A. G., “Hierarchical Predictive Control of an Unmanned Aerial Vehicle Integrated Power, Propulsion, and Thermal Management System,” _IEEE Transactions on Control Systems Technology_ , Vol. 31, No. 3, 2023.  10.1109/TCST.2022.3220913.
  * Koeln et al. [2019] Koeln, J. P., Pangborn, H. C., Williams, M. A., Kawamura, M. L., and Alleyne, A. G., “Hierarchical Control of Aircraft Electro-Thermal Systems,” _IEEE Transactions on Control Systems Technology_ , Vol. 28, No. 4, 2019, pp. 1–15.  10.1109/tcst.2019.2905221.
  * Park and Pangborn [2023] Park, S., and Pangborn, H. C., “Lifted Graph-based Modeling for Linear Predictive Control of Nonlinear Energy Systems,” _2023 IEEE Conference on Control Technology and Applications (CCTA)_ , 2023, pp. 926–933.  10.1109/CCTA54093.2023.10252310.
  * Bruder et al. [2021] Bruder, D., Fu, X., and Vasudevan, R., “Advantages of bilinear Koopman realizations for the modeling and control of systems with unknown dynamics,” _IEEE Robotics and Automation Letters_ , Vol. 6, No. 3, 2021, pp. 4369–4376.
  * Kanai and Yamakita [2023] Kanai, M., and Yamakita, M., “Lifted Bilinear Model-based Linear Model Predictive Control with Scalability,” _IFAC-PapersOnLine_ , Vol. 56, No. 2, 2023, pp. 9405–9410.  https://doi.org/10.1016/j.ifacol.2023.10.232, 22nd IFAC World Congress.
  * Lian et al. [2023] Lian, Y., Jiang, Y., Opila, D. F., and Jones, C. N., “A Proximal-Point Lagrangian-Based Parallelizable Nonconvex Solver for Bilinear Model Predictive Control,” _IEEE Transactions on Control Systems Technology_ , Vol. 31, No. 5, 2023, pp. 2245–2260.  10.1109/TCST.2023.3291537.
  * Budišić et al. [2012] Budišić, M., Mohr, R., and Mezić, I., “Applied Koopmanism,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , Vol. 22, No. 4, 2012, p. 047510.
  * Surana and Banaszuk [2016] Surana, A., and Banaszuk, A., “Linear Observer Synthesis for Nonlinear Systems Using Koopman Operator Framework,” _IFAC-PapersOnLine_ , Vol. 49, No. 18, 2016, pp. 716–723.  10th IFAC Symposium on Nonlinear Control Systems NOLCOS 2016.
  * Mauroy and Mezić [2016] Mauroy, A., and Mezić, I., “Global Stability Analysis Using the Eigenfunctions of the Koopman Operator,” _IEEE Transactions on Automatic Control_ , Vol. 61, No. 11, 2016, pp. 3356–3369.
  * Mauroy et al. [2020] Mauroy, A., Susuki, Y., and Mezić, I., _Koopman operator in systems and control_ , Springer, 2020.
  * Goldschmidt et al. [2021] Goldschmidt, A., Kaiser, E., Dubois, J. L., Brunton, S. L., and Kutz, J. N., “Bilinear dynamic mode decomposition for quantum control,” _New Journal of Physics_ , Vol. 23, No. 3, 2021, p. 033035.
  * Goswami and Paley [2022] Goswami, D., and Paley, D. A., “Bilinearization, Reachability, and Optimal Control of Control-Affine Nonlinear Systems: A Koopman Spectral Approach,” _IEEE Transactions on Automatic Control_ , Vol. 67, 2022, pp. 2715–2728.
  * Jiang et al. [2024] Jiang, X., Li, Y., and Huang, D., “Modularized Bilinear Koopman Operator for Modeling and Predicting Transients of Microgrids,” _IEEE Transactions on Smart Grid_ , 2024.  10.1109/TSG.2024.3399076.
  * Abraham and Murphey [2019] Abraham, I., and Murphey, T. D., “Active learning of dynamics for data-driven control using Koopman operators,” _IEEE Transactions on Robotics_ , Vol. 35, No. 5, 2019, pp. 1071–1083.
  * Rosenfeld and Kamalapurkar [2021] Rosenfeld, J. A., and Kamalapurkar, R., “Dynamic mode decomposition with control Liouville operators,” _IFAC-PapersOnLine_ , Vol. 54, No. 9, 2021, pp. 707–712.
  * Schulze and Mitsos [2022] Schulze, J. C., and Mitsos, A., “Data-Driven Nonlinear Model Reduction Using Koopman Theory: Integrated Control Form and NMPC Case Study,” _IEEE Control Systems Letters_ , Vol. 6, 2022, pp. 2978–2983.
  * Proctor et al. [2016] Proctor, J. L., Brunton, S. L., and Kutz, J. N., “Dynamic Mode Decomposition with Control,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 15, No. 1, 2016, pp. 142–161.  10.1137/15M1013857.
  * Li et al. [2017] Li, Q., Dietrich, F., Bollt, E. M., and Kevrekidis, I. G., “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , Vol. 27, No. 10, 2017.
  * Huang et al. [2023a] Huang, D., Prasetyo, M. B., Yu, Y., and Geng, J., “Learning Koopman Operators with Control Using Bi-level Optimization,” _IEEE CDC 2023_ , 2023a.  10.1109/CDC49753.2023.10383799.
  * Folkestad et al. [2022] Folkestad, C., Wei, S. X., and Burdick, J. W., “KoopNet: Joint learning of Koopman bilinear models and function dictionaries with application to quadrotor trajectory tracking,” _2022 International Conference on Robotics and Automation (ICRA)_ , IEEE, 2022, pp. 1344–1350.
  * Shi and Meng [2022] Shi, H., and Meng, M. Q.-H., “Deep Koopman operator with control for nonlinear systems,” _IEEE Robotics and Automation Letters_ , Vol. 7, No. 3, 2022, pp. 7700–7707.
  * Battaglia et al. [2018] Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., et al., “Relational inductive biases, deep learning, and graph networks,” _arXiv preprint arXiv:1806.01261_ , 2018.
  * Scarselli et al. [2009] Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M., and Monfardini, G., “The graph neural network model,” _IEEE Transactions on Neural Networks_ , Vol. 20, No. 1, 2009, pp. 61–80.  10.1109/TNN.2008.2005605.
  * Hamilton [2020] Hamilton, W. L., _Graph representation learning_ , Morgan & Claypool Publishers, 2020.
  * Yu et al. [2024a] Yu, Y., Harlim, J., Huang, D., and Li, Y., “Learning Coarse-Grained Dynamics on Graph,” _arXiv preprint arXiv:2405.09324_ , 2024a.
  * Yu et al. [2024b] Yu, Y., Jiang, X., Huang, D., Li, Y., Yue, M., and Zhao, T., “PIDGeuN: Graph Neural Network-Enabled Transient Dynamics Prediction of Networked Microgrids Through Full-Field Measurement,” _IEEE Access_ , Vol. 12, 2024b, pp. 49464–49475.  10.1109/ACCESS.2024.3384457.
  * Yu et al. [2022] Yu, Y., Major, D. A., Huang, D., and Schmitz, S., “Reduced-Order Modeling of Ship Airwakes with Atmospheric Turbulence Effects using Dynamic Graph Networks,” _AIAA SCITECH 2022 Forum_ , 2022, p. 2533.
  * Massegur Sampietro and Da Ronch [2024] Massegur Sampietro, D., and Da Ronch, A., “Gradient-Guided Graph Convolutional Multi-Mesh Frameworks for Aircraft Aerodynamics Modelling,” _AIAA SCITECH 2024 Forum_ , 2024, p. 0456.
  * Li et al. [2023] Li, P., Yu, Y., Huang, D., Wang, Z.-H., and Sharma, A., “Regional Heatwave Prediction Using Graph Neural Network and Weather Station Data,” _Geophysical Research Letters_ , Vol. 50, No. 7, 2023, p. e2023GL103405.  https://doi.org/10.1029/2023GL103405, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2023GL103405, e2023GL103405 2023GL103405.
  * Yu et al. [2024c] Yu, Y., Li, P., Huang, D., and Sharma, A., “Street-level temperature estimation using Graph Neural Networks: Performance, feature embedding and interpretability,” _Authorea Preprints_ , 2024c.
  * Sanchez-Gonzalez et al. [2018] Sanchez-Gonzalez, A., Heess, N., Springenberg, J. T., Merel, J., Riedmiller, M., Hadsell, R., and Battaglia, P., “Graph networks as learnable physics engines for inference and control,” _International Conference on Machine Learning_ , PMLR, 2018, pp. 4470–4479.
  * Li et al. [2019] Li, Y., Wu, J., Zhu, J.-Y., Tenenbaum, J. B., Torralba, A., and Tedrake, R., “Propagation networks for model-based control under partial observation,” _2019 International Conference on Robotics and Automation (ICRA)_ , IEEE, 2019, pp. 1205–1211.
  * Karimi et al. [2021] Karimi, A. M., Wu, Y., Koyuturk, M., and French, R. H., “Spatiotemporal Graph Neural Network for Performance Prediction of Photovoltaic Power Systems,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 35, No. 17, 2021, pp. 15323–15330.
  * Lütkepohl [2005] Lütkepohl, H., _VAR Processes with Parameter Constraints_ , Springer Berlin Heidelberg, Berlin, Heidelberg, 2005, Chap. 3, pp. 193–231.  10.1007/978-3-540-27752-1_5.
  * Hochreiter and Schmidhuber [1997a] Hochreiter, S., and Schmidhuber, J., “Long short-term memory,” _Neural computation_ , Vol. 9, No. 8, 1997a, pp. 1735–1780.
  * Yu et al. [2023] Yu, Y., Park, S., Huang, D., and Pangborn, H., “Towards Automated Timescale Decomposition for Graph-based Hierarchical Control of Electrified Aircraft Systems,” _AIAA AVIATION 2023 Forum_ , 2023, p. 4509.
  * Brunton and Kutz [2022] Brunton, S. L., and Kutz, J. N., _Data-driven science and engineering: Machine learning, dynamical systems, and control_ , Cambridge University Press, 2022.
  * Brunton et al. [2016] Brunton, S. L., Brunton, B. W., Proctor, J. L., and Kutz, J. N., “Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control,” _PloS one_ , Vol. 11, No. 2, 2016, pp. 1–19.
  * Messenger and Bortz [2021] Messenger, D. A., and Bortz, D. M., “Weak SINDy for partial differential equations,” _Journal of Computational Physics_ , Vol. 443, 2021, p. 110525.  10.1016/j.jcp.2021.110525, URL http://dx.doi.org/10.1016/j.jcp.2021.110525.
  * Kaptanoglu et al. [2023] Kaptanoglu, A. A., Zhang, L., Nicolaou, Z. G., Fasel, U., and Brunton, S. L., “Benchmarking sparse system identification with low-dimensional chaos,” _Nonlinear Dynamics_ , Vol. 111, No. 14, 2023, pp. 13143–13164.  10.1007/s11071-023-08525-4.
  * Bengio et al. [2013] Bengio, Y., Courville, A., and Vincent, P., “Representation learning: A review and new perspectives,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , Vol. 35, No. 8, 2013, pp. 1798–1828.
  * Hinton and Salakhutdinov [2006] Hinton, G. E., and Salakhutdinov, R. R., “Reducing the dimensionality of data with neural networks,” _Science_ , Vol. 313, No. 5786, 2006, pp. 504–507.
  * Tenenbaum et al. [2000] Tenenbaum, J. B., De Silva, V., and Langford, J. C., “A global geometric framework for nonlinear dimensionality reduction,” _Science_ , Vol. 290, No. 5500, 2000, pp. 2319–2323.
  * Tran et al. [2024] Tran, A., He, X., Messenger, D. A., Choi, Y., and Bortz, D. M., “Weak-form latent space dynamics identification,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 427, 2024, p. 116998.  https://doi.org/10.1016/j.cma.2024.116998, URL https://www.sciencedirect.com/science/article/pii/S0045782524002548.
  * Chen et al. [2018] Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D., “Neural ordinary differential equations,” _arXiv preprint arXiv:1806.07366_ , 2018.
  * Poli et al. [2019] Poli, M., Massaroli, S., Park, J., Yamashita, A., Asama, H., and Park, J., “Graph neural ordinary differential equations,” _arXiv preprint arXiv:1911.07532_ , 2019.
  * Zang and Wang [2020] Zang, C., and Wang, F., “Neural dynamics on complex networks,” _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_, 2020, pp. 892–902.
  * Luo et al. [2023] Luo, X., Yuan, J., Huang, Z., Jiang, H., Qin, Y., Ju, W., Zhang, M., and Sun, Y., “Hope: High-order graph ode for modeling interacting dynamics,” _International Conference on Machine Learning_ , PMLR, 2023, pp. 23124–23139.
  * Huang et al. [2023b] Huang, Z., Zhao, W., Gao, J., Hu, Z., Luo, X., Cao, Y., Chen, Y., Sun, Y., and Wang, W., “TANGO: Time-Reversal Latent GraphODE for Multi-Agent Dynamical Systems,” , 2023b.
  * Gholami et al. [2019] Gholami, A., Keutzer, K., and Biros, G., “ANODE: Unconditionally accurate memory-efficient gradients for neural ODEs,” _arXiv preprint arXiv:1902.10298_ , 2019.
  * Zhuang et al. [2020] Zhuang, J., Dvornek, N., Li, X., Tatikonda, S., Papademetris, X., and Duncan, J., “Adaptive checkpoint adjoint method for gradient estimation in neural ODE,” _arXiv preprint arXiv:2006.02493_ , 2020.
  * Kim et al. [2021] Kim, S., Ji, W., Deng, S., Ma, Y., and Rackauckas, C., “Stiff neural ordinary differential equations,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , Vol. 31, No. 9, 2021.
  * Ross and Fahroo [2002] Ross, I. M., and Fahroo, F., “A perspective on methods for trajectory optimization,” _AIAA/AAS Astrodynamics Specialist Conference and Exhibit_ , 2002, p. 4727.
  * Patterson and Rao [2014] Patterson, M. A., and Rao, A. V., “GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming,” _ACM Transactions on Mathematical Software (TOMS)_ , Vol. 41, No. 1, 2014, pp. 1–37.
  * Pontryagin [2018] Pontryagin, L. S., _Mathematical theory of optimal processes_ , Routledge, 2018.
  * Ross [2015] Ross, I. M., _A Primer on Pontryagin’s Principle in Optimal Control_ , 2nd ed., Collegiate Publishers, 2015.
  * Proctor et al. [2018] Proctor, J. L., Brunton, S. L., and Kutz, J. N., “Generalizing Koopman theory to allow for inputs and control,” _SIAM Journal on Applied Dynamical Systems_ , Vol. 17, No. 1, 2018, pp. 909–930.
  * Defferrard et al. [2017] Defferrard, M., Bresson, X., and Vandergheynst, P., “Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering,” _arXiv:1606.09375_ , 2017.
  * Hochreiter and Schmidhuber [1997b] Hochreiter, S., and Schmidhuber, J., “Long Short-Term Memory,” _Neural Computation_ , Vol. 9, No. 8, 1997b, pp. 1735–1780.  10.1162/neco.1997.9.8.1735.
  * Sutskever et al. [2014] Sutskever, I., Vinyals, O., and Le, Q. V., “Sequence to Sequence Learning with Neural Networks,” _Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2_ , MIT Press, Cambridge, MA, USA, 2014, p. 3104–3112.
  * Paszke et al. [2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S., “PyTorch: An imperative style, high-performance deep learning library,” _Advances in Neural Information Processing Systems 32_ , edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Curran Associates, Inc., 2019, pp. 8024–8035.
  * Aksland et al. [2019] Aksland, C. T., Bixel, T. W., Raymond, L. C., Rottmayer, M. A., and Alleyne, A. G., “Graph-Based Electro-Mechanical Modeling of a Hybrid Unmanned Aerial Vehicle for Real-Time Applications,” _American Control Conference_ , American Automatic Control Council, 2019, pp. 4253–4259.
  * Pangborn et al. [2018] Pangborn, H. C., Koeln, J. P., Williams, M. A., and Alleyne, A. G., “Experimental Validation of Graph-Based Hierarchical Control for Thermal Management,” _Journal of Dynamic Systems, Measurement, and Control_ , Vol. 140, No. 10, 2018, pp. 1–16.