# Physics-informed neural network for quantum control of NMR registers

Priya Batra<EMAIL_ADDRESS>T. S. Mahesh
<EMAIL_ADDRESS>Physics Department and NMR Research Center, Indian
Institute of Science Education and Research, Pune 411008, India

###### Abstract

Classical and quantum machine learning are being increasingly applied to
various tasks in quantum information technologies. Here, we present an
experimental demonstration of quantum control using a physics-informed neural
network (PINN). PINN’s salient feature is how it encodes the entire control
sequence in terms of its network parameters. This feature enables the control
sequence to be later adopted to any hardware with optimal time discretization,
which contrasts with conventional methods involving a priory time
discretization. Here, we discuss two important quantum information tasks: gate
synthesis and state preparation. First, we demonstrate quantum gate synthesis
by designing a two-qubit CNOT gate and experimentally implementing it on a
heteronuclear two-spin NMR register. Second, we demonstrate quantum state
preparation by designing a control sequence to efficiently transfer the
thermal state into the long-lived singlet state and experimentally implement
it on a homonuclear two-spin NMR register. We present a detailed numerical
analysis of the PINN control sequences regarding bandwidth, discretization
levels, control field errors, and external noise.

##  I Introduction

Efficient and robust quantum control is one of the critical needs for the
practical realization of quantum technologies [1, 2]. It deals with the
problem of efficient state-independent unitary synthesis or optimal transfer
of one quantum state to another. Given its importance in a wide variety of
problems, the development of quantum control methodologies is an active area
of research, and numerous optimization algorithms have been designed and
demonstrated [3, 4, 5, 6, 7].

The use of machine learning protocols for different quantum problems has been
a point of significant attention recently [8, 9, 10]. Physics-informed neural
networks (PINNs) are machine-learning algorithms that involve built-in
equations pertaining to certain physical laws [11]. After successfully
implementing a neural network for solving Navier Stokes’s equation [12, 13],
PINN started to gain attention for solving various other types of differential
equations [11, 14]. It is said to be a universal function approximator [15].
Recently, it has been proposed to solve Schrödinger’s equation [16] as well as
the open quantum system Lindbladian equation [17].

An important strength of the PINN algorithm is that it utilizes smooth basis
functions to convert network parameters into the control function. Therefore,
it generates smooth controls with low bandwidth that are easier to implement
in experimental hardware [16]. Hence, PINN avoids the need for any pre-
discretization prevalent in other quantum control algorithms [3, 18, 19].

This work demonstrates PINN-based quantum control for both quantum gate
synthesis and quantum state preparation. First, we utilize PINN to design a
two-qubit CNOT gate and experimentally implement it on a heteronuclear two-
qubit NMR register. We then utilize PINN to design a control sequence to
efficiently transfer the thermal magnetization into the long-lived singlet
state (LLS) and experimentally implement the sequence on a homonuclear two-
qubit NMR register. In both cases, we make detailed numerical analyses of the
control sequences in terms of bandwidth, discretization levels, external
noise, and control errors.

In the next section, we explain the theoretical ideas behind PINN. In Sec. III
and IV, we describe the PINN-based CNOT gate and PINN-based LLS preparation,
respectively. Finally, we conclude in Sec. V.

##  II Physics-informed neural network (PINN)

Consider the total Hamiltonian $H(t)$ consisting of two components, a constant
part $H_{0}$ and a time-dependent part such that

$\displaystyle H(t)=H_{0}+\sum_{k=0}^{M}u_{k,x}(t)I_{k,x}+u_{k,y}(t)I_{k,y},$
(1)

where $u_{k,x(y)}(t)$ are the $x(y)$ time-dependent control amplitudes on the
$k$th nuclear species and $I_{k,x(y)}=\sum_{m}I_{k,x(y)}^{(m)}$ are the
collective spin operators summed over spins of the same isotope. The task here
is to find $u_{k,x(y)}(t)$ to achieve a target unitary $U_{t}$ or a target
state $\rho_{t}$ efficiently. For this purpose, we utilize a deep learning-
based physics-informed neural network (PINN) [11, 14]. We implemented our
neural network with JAX’s Haiku library [20, 21].

Figure 1: Deep neural network and the flowchart describing different steps of
PINN. Here $x_{i}^{l}$ and $b_{i}^{l}$ denote the network nodes and their
biases for the of $l$th layer, $w^{l}_{i,j}$ are the weights connecting $i$th
node of $l$th layer with $j$th node of $(l+1)$th layer. The first layer (input
layer) has a single node which is the time parameter $t$ itself, i.e.,
$x^{1}_{1}=t$. The last layer (output layer) encodes the time-dependent
control parameters $u_{x}(t)=x^{4}_{1}$ and $u_{y}(t)=x^{4}_{2}$ for a single
isotope. The 2nd and 3rd layers are the hidden layers of the (1,3,3,2) neural
network explicitly shown here.

We employ a feed-forward network, as shown in Fig. 1, with the time $t$ as
input and control amplitudes $u_{k,x}(t)$ as output of the network. The node
values are calculated according to

$\displaystyle x^{l+1}_{j}=\Gamma\left(\sum_{i}w^{l}_{i,j}\cdot
x^{l}_{i}+b^{l+1}_{j}\right),$ (2)

where the activation function $\Gamma$ is a nonlinear function that helps
efficiently explore the parameter space. We have used
$\Gamma(\alpha)=\tanh{\alpha}$ as the activation function.

It is important to notice that the control parameters are functions of time
and network parameters. We need to iteratively optimize the network to get the
best controls for any desired target unitary $U_{t}$ or state $\rho_{t}$. In
every iteration, we inject the Hamiltonian $H(t)$ into Schrödinger equation
(in $\hbar=1$ unit)

$\displaystyle\frac{dU(t)}{dt}=-iH(t)U(t),$ (3)

and solve for the overall unitary $U(T)$ using JAX’s ordinary differential
equation solver [22]. Here we use $U(0)=\mathbbm{1}.$ Starting with a set of
random network parameters, we calculate the gate-fidelity

$\displaystyle F_{U}$
$\displaystyle=\absolutevalue{\Tr(U_{t}^{\dagger}U(T))}^{2}.$ (4)

In the case of state-to-state transfer, we use the von Neumann equation

$\displaystyle\frac{d\rho(t)}{dt}=-i[H(t),\rho(t)]$ (5)

and solve for the final achieved state $\rho(T)$ by inputting a definite
initial state $\rho(0)=\rho_{i}$. Here, we use the state-fidelity

$\displaystyle F_{\rho}$ $\displaystyle=\Tr(\rho_{t}\rho(T)).$ (6)

In either case, the goal is to maximize fidelity by optimizing the network. If
the fidelity is less than a threshold value $F_{th}$, then we calculate
gradients $\nabla_{w}F$ and $\nabla_{b}F$ for each weight and bias parameter
using the auto differentiation (autodiff) algorithm in the Optax module [23].
The updated parameters for the next iteration are calculated by adding the
gradients scaled by a step parameter. We optimize the step parameter using the
ADAM optimizer. The entire process is iterated, as shown in Fig. 1 until the
fidelity is maximized.

##  III Gate synthesis: generating CNOT

Single-qubit gates and the two-qubit CNOT gate form a universal set of gates
for quantum computation [24]. As a proof-of-principle demonstration, we employ
PINN to prepare the CNOT gate and implement it on the two-qubit NMR register
diethylfluoromalonate (DEFM) (shown in the inset of Fig. 2) dissolved in
solvent CDCl3. Here 1H and 19F nuclei constitute the two qubits. All
experiments are carried out on a high-resolution 500 MHz NMR spectrometer at
an ambient temperature of 300 K. In on-resonant doubly rotating frame
precessing at the individual Larmor frequencies, the internal Hamiltonian for
the system is

$\displaystyle H_{0}=2\pi JI_{1z}I_{2z},$ (7)

where $J=48.2$ Hz is the scalar coupling and $I_{1z}$, $I_{2z}$ are
respectively the $z$-components of F and H spin operators. The control RF
Hamiltonian is of the following form

$\displaystyle
H_{\mathrm{RF}}(t)=u_{1x}(t)I_{1x}+u_{2x}(t)I_{2x}+u_{1y}(t)I_{1y}+u_{2y}(t)I_{2y},$
(8)

with $u_{1(2)x(y)}$ being the $x(y)$ control amplitudes for qubits $1(2)$. The
total Hamiltonian is thus

$\displaystyle H(t)=H_{0}+H_{\mathrm{RF}}(t).$ (9)

Figure 2: The $x$ and $y$ amplitudes of the PINN-CNOT on 19F (a) and 1H (b).
The reference NMR spectra and the spectra after CNOT for 19F (c) and 1H (d).
The $x$ and $y$ amplitudes of GRAPE-CNOT on 19F (e) and 1H (f). The Fourier
transforms of the PINN and GRAPE sequences on 19F (g) and 1H (h) showing
relatively lower bandwidth of the PINN sequence.

The PINN-generated RF controls for CNOT with 19F as the control qubit and 1H
as the target qubit are shown in Fig. 2 (a,b). Here, we used (1,40,40,4) PINN
with total time set to $T=20$ ms. The experimental 19F and 1H spectra after
the CNOT gate, along with the reference spectra (with $\pi/2$ detection
pulse), are shown in Fig. 2 (c,d).

One major benefit of PINN is that the shape of the control profile is encoded
as a nonlinear function of the network parameters, unlike discrete jumps in
controls obtained by other algorithms such as GRAPE. This has two advantages.
First, the absence of rapid oscillations means smooth control profiles with
minimum bandwidth. To illustrate this, we compare the bandwidth with a GRAPE
control sequence. Fig. 2 (e,f) shows a GRAPE-CNOT sequence (with $2^{7}$
segments), and Fig. 2 (g,h) compares the Fourier transforms of PINN and GRAPE
sequences. It is clear that the PINN sequence has a much lower bandwidth
compared to the GRAPE sequence. The second advantage is that PINN does not
require a priory time-discretization. Once optimized, the neural network can
output the control sequence with any time discretization. This is highly
useful for generating common control profiles and adopting them appropriately
for multiple platforms with different time resolutions. To illustrate this
fact, we now discretize the PINN-CNOT sequence with varying segments and
analyze its infidelity numerically and experimentally. The results shown in
Fig. 3 confirm that very low discretization, with smaller than $2^{5}$
segments, can not faithfully capture the control profiles, thereby leading to
higher infidelities. At $2^{7}$ segments, the experimental infidelity is close
to $10^{-2}$, indicating a good performance of the CNOT gate. However, for
discretizations better than $2^{7}$ segments, the experimental infidelities
show no significant improvement. In other control methods, such as GRAPE, one
needs to run the whole optimization again to improve discretization. The PINN
algorithm avoids under or over-discretizations by capturing the control
profiles in the network parameters, which can output any suitable
discretization depending on the hardware requirement.

Figure 3: The infidelity ($1-F_{U}$) of the PINN-CNOT sequence with varying
discretization size. The experimental infidelity hardly improves for beyond
$2^{7}$ segments, while the simulated infidelity saturates after $2^{1}0$
segments.

##  IV State-to-state transfer: preparation of LLS

Long-lived singlet (LLS) order corresponds to the population difference
between the antisymmetric singlet state
$\ket{S_{0}}=\frac{\ket{01}-\ket{10}}{\sqrt{2}}$ and the symmetric triplet
state $\ket{T_{0}}=\frac{\ket{01}+\ket{10}}{\sqrt{2}}$ in a pair of spin-1/2
nuclei. For the isotropic Hamiltonian $-\mathbf{I}_{1}\cdot\mathbf{I}_{2}$,
the state $\ket{T_{0}}$ forms a triply-degenerate eigenspace together with
other symmetric triplet states $\ket{T_{+}}=\ket{00}$ and
$\ket{T_{-}}=\ket{11}$, while $\ket{S_{0}}$ being the nondegenerate ground
eigenstate. In the absence of symmetry breaking interactions, such as the
chemical shift interaction, it is established that LLS can far outlive the
single-spin relaxation time-constants $T_{1}$ and $T_{2}$ [25]. The survival
of quantum coherence in LLS for such long durations has opened a wide range of
applications from high-precision spectroscopy [26], quantum computing [27],
magnetic resonance imaging [28], etc. Numerous methods have been developed to
transfer thermal state to LLS efficiently [29, 30, 31, 32, 33, 19, 34]. Here
we demonstrate the use of PINN to prepare LLS starting from thermal
magnetization efficiently.

###  IV.1 Experimental preparation of LLS by PINN

We prepare LLS on the two protons of the 2-qubit sample 2,3,6-trichlorophenol
(TCP) (shown in the inset of Fig. 4 (c)) dissolved in (CD3)2SO. The internal
Hamiltonian of the two-qubit system is given by

$\displaystyle H_{0}=2\pi JI_{1z}I_{2z}-\pi\Delta I_{1z}+\pi\Delta I_{2z},$
(10)

where coupling between two qubits is $J=8.75$ Hz and chemical shift difference
is $\Delta=127.5$ Hz. The external RF field has the following Hamiltonian

$\displaystyle H_{RF}(t)=u_{x}(t)(I_{1x}+I_{2x})+u_{y}(t)(I_{1y}+I_{2y}).$
(11)

Thus, the total Hamiltonian is

$\displaystyle H(t)=H_{0}+H_{RF}(t).$ (12)

Figure 4: (a) LLS spectrum (b) Pulse sequence, $x$ and $y$ control amplitudes
for ${}^{1}H$ as a function of time (c) Decay of singlet magnetization under
WALTZ 16 spin-lock. The T1 values of the two protons are 5.5 s and 5.6 s.
Thus, LLS outlives the single-spin non-equilibrium magnetizations by a factor
of around 5. (d) The progression of the expectation value of all four
eigenstates of the singlet-triplet basis. (e) Effect of discretization on
continuous time pulse in simulation as well as experiments.

Here, we used (1,60,60,60,2) PINN with three hidden layers. The PINN sequence
thus generated is shown in Fig. 4 (c). Since the LLS order is not directly
observable, we insert a delay $1/(4\Delta)$ followed by a $\pi/2$-pulse to
convert it to the observable single-quantum magnetization. The resulting
characteristic anti-phase LLS spectrum is shown 4 (a). To measure the lifetime
of LLS, we stored it under the isotropic Hamiltonian realized by the WALTZ-16
spin-lock for up to $44$ sec. The decay during the LLS storage plotted in Fig
4 (b) yields the time-constant $T_{\mathrm{LLS}}=25.5$ which is about five
times the single-spin $T_{1}$ values $5.5$ s and $5.6$ s for the two protons.

To understand the workings of the PINN sequence, we analyze the evolution of
all four eigenstates of the singlet-triplet basis. The progressions of
expectation values of these states during the PINN sequence are plotted in
Fig. 4 (d). The initial thermal state has descending populations from
$\ket{00}=\ket{T_{+}}$ to
$\\{\ket{01},\ket{10}\\}=\frac{1}{\sqrt{2}}(\ket{T_{0}}\pm\ket{S_{0}})$ to
$\ket{11}=\ket{T_{-}}$. Interestingly, after an initial oscillation, the PINN
sequence keeps all four expectation values close to zero until the last 10 ms,
during which oscillations pick up again and population imbalance between
$\ket{S_{0}}$ and $\ket{T_{0}}$, corresponding to high LLS content, is
achieved.

Like in the previous case, here also we analyze the effect of discretization
level on the performance of the PINN-LLS sequence. After the neural network is
optimized, we can discretize the control sequence to any number of segments,
which we have varied from $2^{0}$ to $2^{15}$, and analyze their performance
both by simulations and experiments. The results are shown in Fig. 4 (e).
While the performance was seriously reduced for up to $2^{5}$ due to under-
discretization, the sequence gave good fidelity, and a high LLS order was
prepared at $2^{7}$ segments, which hardly improved thereafter.

###  IV.2 Robustness against noise

It is important for any quantum control method to consider practical
limitations such as external noise. For example, one of the potential
applications of LLS is in magnetic resonance imaging (MRI), where the
longevity of LLS and hyperpolarization help improve the resolution, contrast,
and sensitivity [35]. In such applications, the nuclear spin pair in an
endogenous molecule experiences high external noise, which the state-
preparation sequence should account for. To this end, instead of Eq. 5, we
solve the Lindbladian equation,

$\displaystyle\frac{d\rho(t)}{dt}$ $\displaystyle=-i\,[H(t),\rho(t)]$
$\displaystyle~{}~{}+\gamma\norm{H_{0}}\sum_{n}\biggl{(}V_{n}\,\rho(t)\,V_{n}^{\dagger}-\frac{1}{2}\\{\rho(t),V_{n}^{\dagger}V_{n}\\}\biggr{)},$
(13)

where $\rho(t)=U(t)\rho(0)U^{\dagger}(t)$, $V_{n}$ are the noise operators,
and we have assumed uniform noise coefficient $\gamma\norm{H_{0}}$ scaled by
the norm $\norm{H_{0}}$ of $H_{0}$.

With the above adaptation, we now construct the PINN sequences using the same
Hamiltonian parameters as that of TCP and two noise types: (i) local noise
operators $V_{1}=I_{1x}$, $V_{2}=I_{1y}$,$V_{3}=I_{2x}$, $V_{4}=I_{2y}$, and
(ii) global noise operators $V_{1}=I_{1x}+I_{2x}$, $V_{2}=I_{1y}+I_{2y}$. Fig.
5 (a) and (c) plot the fidelity of LLS vs the noise parameter $\gamma$ for
local and global noises, respectively. We find that the global noise is more
effective in lowering the fidelity. For local noise, we find fairly good
fidelity for up to $\gamma=0.07$, while for global noise, the fidelity drops
drastically beyond 0.04. In addition to making the control sequence robust
against noise, the above Lindbladian approach makes the control sequence
robust against deviations in the control field, either due to miscalibration
or field inhomogeneity. Fig. 5 (b) and (d) plot the fidelity of LLS vs control
field deviation $\Delta u/u$ for various values of noise coefficient $\gamma$.
We find significantly improved robustness in both cases with a small value of
$\gamma$, while we see adverse effects for larger values of $\gamma$. In the
case of global noise, we also observed an oscillatory profile for larger
values of $\gamma$.

Figure 5: (a,c) LLS preparation fidelity vs noise coefficient $\gamma$ for
local (a) and global (c) noise types. (b,d) LLS preparation fidelity vs
control field deviation $\Delta u/u$ for local (b) and global (d) noise types.

##  V Summary and outlook

Elaborate, precise, and robust control over quantum dynamics is important in
all quantum technology applications. Given the ever-growing capabilities of
machine learning, it is a natural pick for the complex challenge of generating
quantum controls. Generally, a quantum control sequence involves time-
dependent Hamiltonians corresponding to various control fields, which makes
finding exact propagators challenging. Most numerical methods for quantum
control design rely on piece-wise constant approximation renders the output
control sequence coarse, higher-bandwidth, and hardware-dependent.

In this work, we employ a physics-informed neural network (PINN) that
overcomes the above problem by encoding the entire shape of the control fields
in terms of nonlinear functions characterized by neural network parameters. As
a result, PINN generates smooth controls, eliminating the need for prior
discretization, and therefore, the control sequence can be adapted to any
hardware with its optimum timing resolution.

Here, we used PINN to design a quantum gate such as CNOT and prepare a quantum
state such as the long-lived singlet (LLS) order in NMR. Using two-qubit NMR
registers, we experimentally confirmed the good performance of CNOT gate and
LLS preparation. We have numerically analyzed the PINN-generated sequences
regarding bandwidth and discretization level. We have also described a
Lindbladian-based approach for generating robust quantum control sequences
against a general noise source. This approach also makes the control sequence
robust against control field errors such as miscalibration or field
inhomogeneities. We hope this work contributes towards advancing the
versatility and effectiveness of quantum control protocols based on machine
learning.

## Acknowledgments

P.B. acknowledges support from the Prime Minister’s Research Fellowship (PMRF)
of the Government of India. The DST/ICPS/QuST/2019/Q67 funding is
acknowledged. We also thank the National Mission on Interdisciplinary Cyber-
Physical Systems for funding from the DST, Government of India, through the
I-HUB Quantum Technology Foundation, IISER-Pune.

## References

  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: a european community view, New Journal of Physics 20, 080201 (2018).
  * Dong and Petersen [2023] D. Dong and I. R. Petersen,  _Learning and Robust Control in Quantum Technology_ (Springer Nature, 2023).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms, J. Magn. Reson. 172, 296 (2005).
  * Nielsen  _et al._ [2007] N. C. Nielsen, C. Kehlet, S. J. Glaser, and N. Khaneja, Optimal control methods in NMR spectroscopy, eMagRes (2007).
  * Caneva  _et al._ [2011] T. Caneva, T. Calarco, and S. Montangero, Chopped random-basis quantum optimization, Phys. Rev. A 84, 022326 (2011).
  * Koch  _et al._ [2022] C. P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Quantum optimal control in quantum technologies. Strategic report on current status, visions and goals for research in Europe, EPJ Quantum Technol. 9, 19 (2022).
  * Mahesh  _et al._ [2023] T. S. Mahesh, P. Batra, and M. H. Ram, Quantum Optimal Control: Practical Aspects and Diverse Methods, J Indian Inst Sci 103, 591 (2023).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C.-H. Wang, A. G. R. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to Machine Learning for physicists, Phys. Rep. A High-Bias, Low-Variance Introduction to Machine Learning for Physicists, 810, 1 (2019).
  * Krenn  _et al._ [2023] M. Krenn, J. Landgraf, T. Foesel, and F. Marquardt, Artificial intelligence and machine learning for quantum technologies, Phys. Rev. A 107, 010101 (2023), publisher: American Physical Society.
  * Batra  _et al._ [2021] P. Batra, A. Singh, and T. S. Mahesh, Efficient Characterization of Quantum Evolutions via a Recommender System, Quantum 5, 598 (2021).
  * Raissi  _et al._ [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations (2017), arXiv:1711.10561.
  * Raissi  _et al._ [2018] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden Fluid Mechanics: A Navier-Stokes Informed Deep Learning Framework for Assimilating Flow Visualization Data (2018), arxiv:1808.04327 [physics, stat] .
  * Raissi  _et al._ [2020] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367, 1026 (2020).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378, 686 (2019).
  * Hornik  _et al._ [1989] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural Networks 2, 359 (1989).
  * Güngördü and Kestner [2022] U. Güngördü and J. P. Kestner, Robust quantum gates using smooth pulses and physics-informed neural networks, Phys. Rev. Res. 4, 023155 (2022).
  * Norambuena  _et al._ [2024] A. Norambuena, M. Mattheakis, F. J. González, and R. Coto, Physics-informed neural networks for quantum control, Phys. Rev. Lett. 132, 010801 (2024).
  * Maximov  _et al._ [2008] I. I. Maximov, Z. Tošner, and N. C. Nielsen, Optimal control design of NMR and dynamic nuclear polarization experiments using monotonically convergent algorithms, The Journal of Chemical Physics 128, 184505 (2008), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2903458/15413228/184505_1_online.pdf .
  * Batra  _et al._ [2020] P. Batra, V. R. Krithika, and T. S. Mahesh, Push-pull optimization of quantum controls, Phys. Rev. Res. 2, 013314 (2020).
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Hennigan  _et al._ [2020] T. Hennigan, T. Cai, T. Norman, and I. Babuschkin, Haiku: Sonnet for JAX (2020).
  * Babuschkin  _et al._ [2020] I. Babuschkin, K. Baumli, A. Bell, S. Bhupatiraju, J. Bruce, P. Buchlovsky, D. Budden, T. Cai, A. Clark, I. Danihelka, A. Dedieu, C. Fantacci, J. Godwin, C. Jones, R. Hemsley, T. Hennigan, M. Hessel, S. Hou, S. Kapturowski, T. Keck, I. Kemaev, M. King, M. Kunesch, L. Martens, H. Merzic, V. Mikulik, T. Norman, G. Papamakarios, J. Quan, R. Ring, F. Ruiz, A. Sanchez, R. Schneider, E. Sezener, S. Spencer, S. Srinivasan, W. Stokowiec, L. Wang, G. Zhou, and F. Viola, The DeepMind JAX Ecosystem (2020).
  * Baydin  _et al._ [2017] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: A survey, J. Mach. Learn. Res. 18, 5595 (2017).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Carravetta  _et al._ [2004] M. Carravetta, O. G. Johannessen, and M. H. Levitt, Beyond the ${T}_{1}$ limit: Singlet nuclear spin states in low magnetic fields, Phys. Rev. Lett. 92, 153003 (2004).
  * Pileio [2020] G. Pileio,  _Long-lived Nuclear Spin Order: Theory and Applications_ (The Royal Society of Chemistry, 2020).
  * Roy and Mahesh [2010] S. S. Roy and T. S. Mahesh, Initialization of nmr quantum registers using long-lived singlet states, Phys. Rev. A 82, 052302 (2010).
  * Vasos  _et al._ [2009] P. R. Vasos, A. Comment, R. Sarkar, P. Ahuja, S. Jannin, J.-P. Ansermet, J. A. Konter, P. Hautle, B. van den Brandt, and G. Bodenhausen, Long-lived states to sustain hyperpolarized magnetization, Proceedings of the National Academy of Sciences 106, 18469 (2009), https://www.pnas.org/doi/pdf/10.1073/pnas.0908123106 .
  * Carravetta and Levitt [2004] M. Carravetta and M. H. Levitt, Long-lived nuclear spin states in high-field solution nmr, Journal of the American Chemical Society 126, 6228 (2004), pMID: 15149209, https://doi.org/10.1021/ja0490931 .
  * Feng  _et al._ [2012] Y. Feng, R. M. Davis, and W. S. Warren, Accessing long-lived nuclear singlet states between chemically equivalent spins without breaking symmetry, Nature physics 8, 831 (2012).
  * DeVience  _et al._ [2013] S. J. DeVience, R. L. Walsworth, and M. S. Rosen, Preparation of nuclear spin singlet states using spin-lock induced crossing, Phys. Rev. Lett. 111, 173002 (2013).
  * Pravdivtsev  _et al._ [2016] A. N. Pravdivtsev, A. S. Kiryutin, A. V. Yurkovskaya, H.-M. Vieth, and K. L. Ivanov, Robust conversion of singlet spin order in coupled spin-1/2 pairs by adiabatically ramped rf-fields, Journal of Magnetic Resonance 273, 56 (2016).
  * Khurana and Mahesh [2017] D. Khurana and T. Mahesh, Bang-bang optimal control of large spin systems: Enhancement of 13c–13c singlet-order at natural abundance, Journal of Magnetic Resonance 284, 8 (2017).
  * Suresh  _et al._ [2023] A. Suresh, V. Varma, P. Batra, and T. S. Mahesh, Counterdiabatic driving for long-lived singlet state preparation, The Journal of Chemical Physics 159, 024202 (2023), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0159448/18041352/024202_1_5.0159448.pdf .
  * Laustsen  _et al._ [2012] C. Laustsen, G. Pileio, M. C. D. Tayler, L. J. Brown, R. C. D. Brown, M. H. Levitt, and J. H. Ardenkjaer-Larsen, Hyperpolarized singlet nmr on a small animal imaging system, Magnetic Resonance in Medicine 68, 1262 (2012), https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.24430 .


# Learning of Causal Observable Functions for Koopman-DFL Lifting
Linearization of Nonlinear Controlled Systems and Its Application to
Excavation Automation

Nicholas Stearns Selby and H. Harry Asada, Manuscript received: March 1, 2021;
Revised May 20, 2021; Accepted June 14, 2021. This paper was recommended for
publication by Editor Youngjin Choi upon evaluation of the Associate Editor
and Reviewers’ comments. This material is based upon work supported by
National Science Foundation Grant NSF-CMMI 2021625. The authors are with the
School of Engineering, Massachusetts Institute of Technology, Cambridge, MA,
USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Object Identifier (DOI):
10.1109/LRA.2021.3092256.© 2021 IEEE. Personal use of this material is
permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for
advertising or promotional purposes, creating new collective works, for resale
or redistribution to servers or lists, or reuse of any copyrighted component
of this work in other works.

###### Abstract

Effective and causal observable functions for low-order lifting linearization
of nonlinear controlled systems are learned from data by using neural
networks. While Koopman operator theory allows us to represent a nonlinear
system as a linear system in an infinite-dimensional space of observables,
exact linearization is guaranteed only for autonomous systems with no input,
and finding effective observable functions for approximation with a low-order
linear system remains an open question. Dual-Faceted Linearization uses a set
of effective observables for low-order lifting linearization, but the method
requires knowledge of the physical structure of the nonlinear system. Here, a
data-driven method is presented for generating a set of nonlinear observable
functions that can accurately approximate a nonlinear control system to a low-
order linear control system. A caveat in using data of measured variables as
observables is that the measured variables may contain input to the system,
which incurs a causality contradiction when lifting the system, i.e. taking
derivatives of the observables. The current work presents a method for
eliminating such anti-causal components of the observables and lifting the
system using only causal observables. The method is applied to excavation
automation, a complex nonlinear dynamical system, to obtain a low-order lifted
linear model for control design.

###### Index Terms:

Robotics and Automation in Construction; Reinforcement Learning

##  I Introduction

There is a growing need in the construction and mining industries for
excavation automation. Various technologies are being developed for operating
excavators autonomously with increased productivity and fuel efficiency [1].
The recent and projected growth of the global construction industry [2] and
the dangers of the excavation work environment [3] are major drivers behind
the development of intelligent excavators for performing earth-moving tasks.

Excavation is a highly nonlinear process where soil and rocks interact with
the bucket of an excavator in a complex manner (see Fig. 1). While
terramechanics models have been studied for many decades, their validity is
limited due to the difficulty of identifying the numerous parameters of
mechanistic models. Data-driven methods have recently been introduced to
autonomous excavation for capturing complex nonlinearities [4, 5, 6, 3, 7],
yet the nonlinear models are still too complex to use, in particular, for
real-time control.

Figure 1: Learned lifting linearization of autonomous excavation.

Lifting linearization is a methodology for representing a nonlinear dynamical
system with a linear dynamic model in a high-dimensional space. Underpinned by
Koopman operator theory, nonlinear systems represented with supernumerary
state variables behave more linearly in the lifted space. The method has
recently been applied to various robotics and automation challenges, including
active learning [8], soft robotics [9], human-robot interaction [10], power
systems [11], and mission planning [12]. More broadly, deep learning has
proven a valuable tool for lifting linearization techniques [13, 14, 15, 16].

The original Koopman Operator has two major limitations:

  1. 1.

The theory is applicable only to dynamical systems with no exogenous input,
i.e. autonomous systems, and

  2. 2.

Exact linearization requires an infinite-dimensional space, except for a
restricted class of systems.

Any extension to non-autonomous, finite-dimensional systems is no longer
exact, but an approximation. Various methods for truncating the system with a
finite-dimensional space have been reported. Among others, the
eigendecomposition of the lifted system allows us to represent the system at
desirable accuracy and granularity while providing useful insights into the
system [17]. Furthermore, the extended Dynamic Mode Decomposition (eDMD) is
completely data-driven, providing a practical tool for complex nonlinear
system representation [18]. These methods, however, need a set of observables,
i.e. output variables, which are nonlinear functions of independent state
variables. It is still an open question how to find an effective set of
observable functions.

One of the key challenges in the lifting linearization of nonlinear systems
with exogenous input is causality. If observable functions are functions of
both state variables and input variables, we cannot use such observables for
lifting the system. Lifting entails computing time derivatives of the
observables and, thereby, the dynamic equations inevitably include the time
derivative of input. In discrete formulation, this means the use of future
input. Including these input terms, the time-evolution of the observables
turns out not to be causal. If one measures a set of candidates of observables
from a nonlinear system that is subject to control inputs and uses the
measured variables for lifting the system, they may end up with a non-causal
dynamical system.

In Dual-Faceted Linearization (DFL), another approach to lifting
linearization, the causality issue is analyzed based on physical system
modeling theory [19]. In DFL, the propagation of inputs across the nonlinear
dynamical system can be tracked, and their effect on all observables, called
auxiliary variables, can be localized. Assuming that inputs are linearly
involved in observables, a method has been established for eliminating the
input-dependent component from each observable and lifting the dynamics using
the remaining input-free observables. In the Koopman-based lifting
linearization, too, it is assumed that the observables are input-affine in
order to eliminate input-dependent components from observable functions so a
causal dynamic model can be obtained [20].

Lifting linearization is a powerful methodology for tackling a broad spectrum
of nonlinear problems, in particular, excavation process modeling and control.
However, two critical challenges have not yet been fully solved:

  1. 1.

Finding an effective set of observables to approximate a nonlinear system in a
low-dimensional lifted space

  2. 2.

Finding causal observables uncorrelated with inputs

The objective of the current work is to solve these two challenges. We present
a low-dimensional, causal, lifting linear model obtained from experimental
data. Neural networks are used to find effective observables through learning.

In the following, we summarize a basic formulation of lifting linearization in
§ II. We present the learning method for obtaining an effective set of
observables in § III. First, we deal with nonlinear controlled systems where
all measured observables are not affected by inputs. Then, the method is
extended for physical observables that may be functions of inputs. Simple
numerical examples are discussed for validating the proposed method in § IV,
and we apply it to excavation process modeling in § IV-B.

##  II Background

This section summarizes background knowledge for readability. More details can
be found in [21, 19].

###  II-A Koopman Operator Theory

First proposed in 1931 by Koopman [21], Koopman operator theory originally
modeled autonomous systems by mapping nonlinear dynamics onto an infinite-
dimensional linear space. Later techniques expanded the use of the operator
for nonautonomous systems [22] and developed methods for approximating the
infinite-dimensional mapping with a computationally feasible, finite-
dimensional space [23, 24, 25].

Let the discrete-time dynamics of a nonlinear, autonomous system with state
$x_{t}\in\mathbb{R}^{l}$ at time $t$ be given by $x_{t+1}=f(x_{t})$.
Furthermore, define a vector of nonlinear observables of the state,
$\eta_{t}=g(x_{t})\in\mathbb{R}^{m}$.

The Koopman operator, $\mathcal{K}$, is linear and infinite-dimensional and
applies to observable functions: $\mathcal{K}_{f}g=g\circ f$, where $\circ$
represents the composition operator.

###  II-B Dual-Faceted Linearization (DFL)

Despite the use of Koopman operators to provide a lifting linearization for
autonomous systems, the theory provides no method by which to select an
effective set of observables. Because it is infeasible to compute the
infinite-dimensional space with finite computational resources, the choice of
which observables to use is very important. DFL [19] uses a particular class
of observables that are determined based on physical modeling theory and bond
graphs [26]. Those observables, called auxiliary variables, are physically
meaningful, and may be measured physically. Furthermore, causality analysis of
the method allows us to examine how exogenous inputs propagate the system and
influence specific auxiliary variables. Using those variables with no input
influence, one can obtain a lifted system that is causal. Alternatively,
input-dependent variables can be “laundered” into causal variables.

Consider the discrete-time dynamics of a nonlinear, nonautonomous system with
input $u_{t}\in\mathbb{R}^{n}$ at time $t$ given by:

$x_{t+1}=f(x_{t},u_{t})$ (1)

Assuming that the system is a lumped-parameter system with integral causality,
we can choose outputs of all the nonlinear elements involved in the system as
observables $\eta_{t}=g(x_{t})\in\mathbb{R}^{m}$ to augment the system state
and construct a linear representation of the system dynamics:

$x_{t+1}=A_{x}x_{t}+A_{\eta}\eta_{t}+B_{x}u_{t}$ (2)

where $A_{x}\in\mathbb{R}^{l\times l}$, $A_{\eta}\in\mathbb{R}^{l\times m}$,
and $B_{x}\in\mathbb{R}^{l\times n}$ are fixed matrix coefficients determined
by the physical structure of the system. This part of the state evolution is
exact.

We approximate the $\eta$-dynamics using a second equation:

$\eta_{t+1}=H_{x}x_{t}+H_{\eta}\eta_{t}+H_{u}u_{t}+r_{\eta_{t+1}}$ (3)

where $H_{x}\in\mathbb{R}^{m\times l}$, $H_{\eta}\in\mathbb{R}^{m\times m}$,
and $H_{u}\in\mathbb{R}^{m\times n}$ are fixed matrix coefficients and
$r_{\eta_{t+1}}\in\mathbb{R}^{m}$ is a residual. Unlike $A_{x}$, $A_{\eta}$,
and $B_{x}$, which are determined from the physical structure of the system,
$H_{x}$, $H_{\eta}$, and $H_{u}$ must be regressed from data. For brevity,
define coefficient matrices
$A\triangleq\left(A_{x},A_{\eta},B_{x}\right)\in\mathbb{R}^{l\times p}$ and
$H\triangleq\left(H_{x},H_{\eta},H_{u}\right)\in\mathbb{R}^{m\times p}$ and
datum vector
$\xi_{t}\triangleq\left(x_{t}^{\intercal},\eta_{t}^{\intercal},u_{t}^{\intercal}\right)^{\intercal}\in\mathbb{R}^{p}$
where $p=l+m+n$. Apply a negative discrete-time shift operator $T_{-1}$ to (3)
to optimize $H$ to minimize the mean squared error of predicting $\eta_{t+1}$:

$\begin{array}[]{rl}H^{\mathrm{o}}&=\operatorname*{arg\,min}_{H}\mathrm{E}\left[\left|H\xi_{t-1}-\eta_{t}\right|^{2}\right]\\\
&=\mathrm{E}\left[\eta_{t}\xi_{t-1}^{\intercal}\right]\left(\mathrm{E}\left[\xi_{t-1}\xi_{t-1}^{\intercal}\right]\right)^{-1}\end{array}$
(4)

where $\mathrm{E}[\cdot]$ is the expectation operator. Assuming that the
system is persistently excited and that $u_{t}$ is not collinear with $x_{t}$,
there is a unique solution, $H^{\mathrm{o}}$.

Figure 2: Block diagram of the learned lifting linearization algorithm. The
loss, $J$, is used to tune the weights of neural network, $g$, and linear
dynamic model matrices $A$ and $H$ via error backpropagation (EBP). Note that
both instances of $g$ are equivalent to each other for all time $t$.

The original nonlinear dynamics $f$ can now be modeled using the dual-faceted
linear dynamics:

$x_{t+1}=A\xi_{t};\quad\eta_{t+1}\approx H\xi_{t};\quad\eta_{0}=g(x_{0});$ (5)

Practical benefits of using DFL as a tool to model systems include:

  * •

Augmented state feedback can be used to better inform controllers [27].

  * •

Linear observer design is enabled for augmented state feedback.

  * •

Model-predictive control is convex [17].

  * •

Because DFL is based in physical modelling theory, augmented state systems may
be measurable, and dynamics may have physical intuition [28].

For these reasons, DFL has proven to be a valuable tool in modeling nonlinear
systems. However, DFL requires knowledge of the structure of the physical
system. Furthermore, many systems contain no obvious, measurable observables
with which to augment the state, and there is no guarantee that, when they do
exist, physically meaningful observables make the best choices for augmenting
the system state.

###  II-C Machine Learning for Linear Latent Spaces

Like Koopman operator theory and DFL, learned latent-space dynamic modeling
techniques also involve constructing a nonlinear representation of the
original state, then using a model to evolve the new “observables” through
time. In Koopman operator theory and DFL, the dynamic model is linear.
Recently, much work has been done to explore applying deep learning techniques
to Koopman operators.

Abraham and Murphey [8] presented an active learning strategy for robotic
systems that extended observables, $g$, from Koopman operator theory to
include the control input, $u_{t}$. Their algorithm trains a neural network to
approximate an optimal $g$. In each epoch, they recompute a finite-dimensional
matrix approximation of the Koopman operator $\mathcal{K}$, but because $g$ is
a function of both $x_{t}$ and $u_{t}$, regressing such a matrix requires
knowing $u_{t+1}$ in addition to $x_{t+1}$. To solve this causality problem,
they propose replacing $u_{t+1}$ with $u_{t}$ in the Koopman operator update.

Han et al. [29] also proposed using a neural network that approximates an
optimal lifting function $g$. At the end of each epoch, after feeding
$x_{t-1}$ forward through $g$, they solve a least squares optimization problem
to regress a modified $H^{\mathrm{o}}$ where $H_{x}=0$, then use the learned
model to backpropagate the error
$\mathrm{E}\left[\left|\left|g(x_{t})-H^{\mathrm{o}}\xi_{t-1}\right|\right|_{F}\right]$
plus a penalty on the norms of the components of $H^{\mathrm{o}}$. Because
they do not track the evolution of the state, $x$, directly, they must
simultaneously learn an additional matrix to approximate $g^{-1}$.

Lusch et al. [30] and Mastia and Bemporad [31] replaced the $g^{-1}$ matrix
approximation of Han et al. with a neural network decoder, while in this work
we handle this by feeding the state forward directly into the linear dynamic
model. Their models learn to minimize cost along three axes: neural network
$g^{-1}$ must invert neural network $g$, i.e. $x=g\left(g^{-1}(x)\right)$; the
nonlinear model must be able to predict, i.e.
$x_{t+1}=g^{-1}\left(Hg(x_{t})\right)$; and the lifted state must propagate
linearly, i.e. $g(x_{t+1})=Hg(x_{t})$.

Work by Yeung et al. [32] applied deep learning to dynamic mode decomposition
(DMD) by training a neural network to approximate the “snapshot” function
mapping state to physical observable. Instead of optimizing a complicated cost
function like in [29, 30, 31], their work simply trains a model to map states
to other measurements of the system, then regresses a linear dynamic model to
propagate the measurements forward in time.

Other work in learned latent spaces for lifting linearizations leverages
neural networks to approximate functions similar to Koopman’s observables in
reinforcement learning [16], sampling-based motion planning [14], Kalman
filtering [15], and partially observable Markov decision process [13]
paradigms. However, these works all learn linearizations $A$ and $H$ that are
time- or state-dependent, further reducing the long-term robustness of the
linear model and giving up the benefits of provably convex optimal control.

###  II-D Anticausal Observables

Most lifting linearization techniques, including Koopman operator theory and
DFL, require that the lifting observables be control input-independent, i.e.
$\eta=g(x)$. If auxiliary variables depend on the control input, $u$, i.e.
$\eta=g(x,u)$, then propagating $\eta$ forward through time requires knowledge
of future values of control input:

$\eta_{t+1}\approx\eta_{t}+\frac{\partial\eta}{\partial
x}\left(x_{t+1}-x_{t}\right)+\frac{\partial\eta}{\partial
u}\left(u_{t+1}-u_{t}\right)$ (6)

In order to avoid problems with causality, most methods explicitly avoid
augmenting the system state with control input-dependent variables. There are
two common techniques.

The first solution to the causality problem is to include a state-feedback
control law in the model [8, 29]. By constraining the control input to be a
known function of state, the system becomes effectively autonomous, and the
original formulation of Koopman operator theory applies. This works for
regulators, including for controllers regulating a system to follow a
predetermined state trajectory, but no exogenous input is allowed.

The second solution to the causality problem is to assume that the auxiliary
variables are linear in $u$ [20, 19, 33]:

$\eta(x,u)=\eta^{*}(x)+Du$ (7)

where $\eta^{*}$ is exclusively state-dependent and $D$ is a fixed matrix
coefficient of $u$.

Although predicting future values of $\eta(x,u)$ remains impossible without a
control law, this formulation allows for the modeling of the evolution of
$\eta^{*}(x)$. The auxiliary state equation can be rewritten:

$\eta^{*}_{t+1}=H_{x}^{*}x_{t}+H_{\eta}^{*}\eta_{t}+H_{u}^{*}u_{t}+r_{\eta^{*},t}$
(8)

where $r_{\eta^{*},t}$ is a residual.

Substituting (7) into (8) yields:

$\eta^{*}_{t+1}=H_{x}^{*}x_{t}+H_{\eta}^{*}\eta^{*}_{t}+\left(H_{u}^{*}+H_{\eta}^{*}D\right)u_{t}+r_{\eta^{*},t}$
(9)

which is a causal, augmented state dynamic equation. The question of causality
is therefore solved by preprocessing the auxiliary state data to filter out
their dependence on $u$.

##  III Modeling Algorithm

Fig. 2 shows the overview of the modeling algorithm. The learning system
consists of three major components. The first is the Linear Dynamic Model
predicting the transition of the system. The second is a multi-layer neural
network for generating auxiliary variables, and the third is the process of
evaluating the prediction error.

We assume all state variables are accessible and sufficient data for training
are attainable. All state variables are fed into the left neural network, $g$,
to produce a set of synthetic observables, $\eta_{t}$, to be learned. The
state $x_{t}$, observables $\eta_{t}$, and input $u_{t}$ are fed into the
linear dynamic model parameterized by matrices $A$ and $H$. The linear model
produces predicted state $\hat{x}_{t+1}$ and predicted observables
$\hat{\eta}_{t+1}$. The predicted state and observables are compared to their
ground truth values, $x_{t+1}$ and $\eta_{t+1}=g(x_{t+1};\theta)$,
respectively. The right neural network in the figure is a twin copy of the
left neural network, sharing parameters $\theta$.

The squared error of the predicted state and observables, $J$, is used for
updating the linear dynamic model with respect to parameters $A$ and $H$ and
the neural network weights. The update of these parameters is computed via
error backpropagation (EBP).

The causality analysis involved in the DFL modeling allows us to examine
whether observables, called auxiliary variables, are functions of state alone
or include inputs. If some auxiliary variables are causal, having no
dependence on control input, they can be added to the state for lifting the
dynamics. Let $\zeta^{*}_{t}$ represent causal auxiliary variables. As shown
in Fig. 2, the causal observables can be fed into the neural network, so that
the synthetic auxiliary variables $\eta_{t}$ can be produced from richer data.
Note that the ground truth $\eta_{t+1}$, too, can be produced in response to
not only state $x_{t+1}$ but also $\zeta_{t+1}^{*}$. The tunable parameter
space of the neural network is expanded with the use of the causal auxiliary
variables.

In case physically measurable auxiliary variables are functions of both state
and input, such input-dependent auxiliary variables cannot be used in their
original form for lifting the dynamics. It is necessary to filter out the
input components from the observables. Fig. 2 also shows a simple filter to
eliminate the effect of input from those variables, $\zeta_{t}$, as discussed
in § III-B.

###  III-A Discrete-Time Learned Lifting Linearization

Consider the discrete-time dynamic system from (1). Let $g$ be a neural
network, illustrated in Fig. 3, defined by randomly initialized parameters
$\theta$ to generate synthetic observables:

$\eta_{t}^{\theta}=g(x_{t};\theta)$ (10)

Figure 3: Diagram of the neural network and linear dynamic model to compute
$x_{t+1}$ and $\eta_{t+1}$ given $x_{t}$ and $u_{t}$. With abuse of notation,
we hereafter include $\zeta^{*}$ in $x$.

###  III-B Extension to Anticausal Observables

Define a datum vector
$\xi_{t}^{\theta}\triangleq\left(x_{t}^{\intercal},\eta_{t}^{\theta\intercal},u_{t}^{\intercal}\right)^{\intercal}$.
Let $A\in\mathbb{R}^{l\times p}$ and $H\in\mathbb{R}^{m\times p}$ be matrix
coefficients modeling the state and augmented state transition dynamics,
respectively. We override (5) to include residuals in the original and
augmented state dynamic equations:

$\begin{cases}x_{t+1}=A\xi_{t}^{\theta}+r_{x_{t+1}}^{\theta}\\\
\eta_{t+1}^{\theta}=H\xi_{t}^{\theta}+r_{\eta_{t+1}}^{\theta}\\\
\eta_{0}^{\theta}=g(x_{0};\theta)\end{cases}$ (11)

Given observation data of $x_{t}$, $x_{t-1}$, and $u_{t-1}$, we synthesize
observations of the augmented state, $\eta_{t}^{\theta}=g(x_{t};\theta)$ and
$\eta_{t-1}^{\theta}=g(x_{t-1};\theta)$, and assemble datum vectors
$\xi_{t-1}^{\theta}$.

By applying the discrete-time shift operator $T_{-1}$ to (11) and rearranging,
we can compute a residual for each observation:
$r_{x_{t}}^{\theta}=x_{t}-A\xi_{t-1}^{\theta}$ and
$r_{\eta_{t}}^{\theta}=g(x_{t};\theta)-H\xi_{t-1}^{\theta}$.

We define a quadratic loss function, $J_{t}(\theta,A,H)$, used to train the
model:

$J_{t}(\theta,A,H)\triangleq r_{t}^{\theta\intercal}Qr_{t}^{\theta}$ (12)

where $Q$ is a symmetric matrix coefficient and $r_{t}^{\theta}$ is a total
residual given by
$r_{t}^{\theta}\triangleq\left(r_{x_{t}}^{\theta\intercal},r_{\eta_{t}}^{\theta\intercal}\right)^{\intercal}$.

Model parameter matrices $A$ and $H$, as well as the parameters of the neural
network, $\theta$, are computed by solving the following optimization problem
via error backpropagation:

$\theta^{o},A^{o},H^{o}=\operatorname*{arg\,min}_{\theta,A,H}\mathrm{E}\left[J_{t}(\theta,A,H)\right]$
(13)

As discussed in [19], augmenting the state with physical observables is often
useful. Because this learned lifting linearization is data-driven, augmenting
the state, $x$, with control-independent, physical observables is trivial.
However, as reviewed in § II, if the augmented state is dependent on the
control input, the augmented system dynamics become anticausal. In state space
modeling, output equations include inputs algebraically if there is a direct
transmission term from inputs to outputs [34, 35]. Namely, observations of the
system are functions of $x_{t}$ and $u_{t}$. Consider a vector of physical
observables, $\zeta(x,u)\in\mathbb{R}^{z}$, suspected of including a
dependence on $u$. As in [20, 19], and [33], assume that this dependence is
linear:

$\zeta(x,u)=\zeta^{*}(x)+Du$ (14)

where $\zeta^{*}(x)\in\mathbb{R}^{z}$ is exclusively a function of state and
$D\in\mathbb{R}^{z\times n}$ is a fixed matrix coefficient of $u$. Assuming
mean-zero data, because $\zeta^{*}(x)$ is uncorrelated with the control input,
$\mathrm{E}\left[\zeta^{*}(x)u^{\intercal}\right]=0$. Therefore, multiplying
(14) by $u^{\intercal}$ and taking the expectation yields
$\mathrm{E}\left[\zeta(x,u)u^{\intercal}\right]=D\mathrm{E}\left[uu^{\intercal}\right]$.

Given observations of $\zeta(x,u)$ and $u$, computing $D$ becomes a least-
squares linear regression:

$\hat{D}=\mathrm{E}\left[\zeta(x,u)u^{\intercal}\right]\mathrm{E}\left[uu^{\intercal}\right]^{-1}$
(15)

assuming that the input is persistently exciting.

Before training the learned lifting linearization model in (11), we preprocess
the data to “clean” the physical observables from any linear dependence on $u$
via

$\zeta^{*}(x_{t})=\zeta(x_{t},u_{t})-\hat{D}u_{t}$ (16)

for each observation, $t$. Then, we augment the DFL model to include
$\zeta^{*}(x)$. We override (10) with
$\eta_{t}^{\theta}=g\left(\left(\zeta^{*\intercal}(x_{t}),x_{t}^{\intercal}\right)^{\intercal};\theta\right)$
and follow the same training procedure described above to tune $g$, $A$, and
$H$ using (13).

The complete learned lifted linearization algorithm is summarized in Algorithm
1.

Result: Lifting linearization of nonlinear dynamics

Randomly initialize neural network $g$ and linear dynamic model $A$, $H$;

$\hat{D}\leftarrow\mathrm{E}\left[\zeta
u^{\intercal}\right]\mathrm{E}\left[uu^{\intercal}\right]^{-1}$ ;

while  _training_ do

get batch of $x_{t}$, $\zeta_{t}$, $u_{t}$, $x_{t+1}$, $\zeta_{t+1}$,
$u_{t+1}$ from training dataset ;

$\zeta^{*}_{t}\leftarrow\zeta_{t}-\hat{D}u_{t}$ ;

$\zeta^{*}_{t+1}\leftarrow\zeta_{t+1}-\hat{D}u_{t+1}$ ;

$\eta_{t}\leftarrow g((x_{t},\zeta^{*}_{t}))$ ;

$\eta_{t+1}\leftarrow g((x_{t+1},\zeta^{*}_{t+1}))$ ;

$\xi_{t}\leftarrow\left(x_{t}^{\intercal},\zeta^{*\intercal}_{t},\eta_{t}^{\intercal},u_{t}^{\intercal}\right)^{\intercal}$
;

$r\leftarrow\left(((x_{t+1}^{\intercal},\zeta^{*\intercal}_{t+1})^{\intercal}-A\xi_{t})^{\intercal},(\eta_{t+1}-H\xi_{t})^{\intercal}\right)^{\intercal}$
;

$J\leftarrow r^{\intercal}Qr$ ;

backpropagate $J$ to update $g$, $A$, and $H$ using Adam ;

end while

$A_{u}\leftarrow A_{u}+A_{\zeta}D$ ;

$H_{u}\leftarrow H_{u}+H_{\zeta}D$ ;

Algorithm 1 Learned Lifting Linearization

##  IV Numerical Examples

###  IV-A Toy Problem

The modeling algorithm in § III is implemented in PyTorch [36] on a laptop
running Ubuntu 18.04.5 LTS. The codebase is hosted as a git repository at
[37].

Figure 4: Bond graph of a nonlinear first-order system with state variable
$q$. In the bond graph, a nonlinear capacitor, C, and a nonlinear resistor, R,
are connected to an effort source, $\mathrm{S_{e}}$, that is an exogenous
input $u(t)$. Causality analysis of the bond graph determines that effort
variable $\mathrm{e_{C}}$ is the output of the nonlinear capacitor, while the
output of the nonlinear resistor is flow variable $f$. In the electrical
circuit analogy, the effort variable $\mathrm{e_{C}}$ is the voltage across
the capacitor, and the flow variable f is the current flowing through the
resistor. They are connected with the exogenous input voltage $u(t)$ at the
“1” junction, which is equivalent to Kirchhoff’s Voltage Law. The causality
analysis also reveals that a direct transmission path exists from input $u(t)$
to flow variable f. Therefore, auxiliary variable f is not a causal variable
for lifting the system. For more details about causality analysis, see [19,
26]. Figure 5: Results of open-loop simulation predicting the state, $q$ of
the toy problem from Fig. 4 excited by a square wave input (in gray). The
solid black line indicates the true trajectory. The dashed lines indicate the
open-loop simulated trajectories of various models: Koopman-with-control
(Koopman, dim$\ =33$), extended dynamic mode decomposition with control
(eDMDc, dim$\ =6$), dual-faceted linearization (DFL, dim$\ =4$), learned
lifting linearization (L3, dim$\ =6$), and L3 without first filtering out the
control input using (16) (L3 NoF, dim$\ =6$).

We test the learned lifting linearization (L3) algorithm on the nonlinear,
massless spring-damper illustrated in Fig. 4 with
$\Phi_{\mathrm{R}}(\mathrm{e_{R}})=2/\left(1+e^{-4\mathrm{e_{R}}}\right)-1$,
and $\Phi_{\mathrm{C}}(q)=\mathrm{sgn}(q)q^{2}$. We generate 100 5s
trajectories at 20Hz with initial conditions and control inputs drawn from
uniform random distributions. The state, $x$, consists only of the linear
position, $q$, and the control input, $u$, is the scalar effort ($n=1$). We
use the system bond graph to identify observables
$\zeta=\left(\mathrm{f,e_{C}}\right)^{\intercal}$.

The neural network, $g$, approximating the optimal synthetic observables,
$\eta$, is a fully connected network of 3 linear input neurons ($l=3$); two
hidden layers, each with 256 ReLU neurons; and 2 output neurons, creating 2
synthetic observables ($m=2$). Before training, the data are randomly divided
80-20 into a training set and a validation set. The neural network, $g$, and
the linear model consisting of $A$ and $H$ are trained in batches of 32 input-
output pairs using an Adam optimizer [38] with $\alpha=10^{-5}$,
$\beta_{1}=0.9$, $\beta_{2}=0.999$, and $\epsilon=10^{-8}$. The quadratic cost
parameter $Q=I$. Before each training epoch, the learned lifting linearization
model is evaluated using the validation dataset without backpropagating the
loss. Training continues until the validation loss begins to increase.

TABLE I: Integrated squared error of the models used to simulate the toy problem over ten seconds. Koopman | eDMDc | DFL | L3 | L3 (NoF)
---|---|---|---|---
5.9 | 14 | 0.73 | 0.48 | 12

We benchmark the learned lifting linearization algorithm against Koopman-with-
control, eDMDc, DFL, and L3 without the anticausal filter. Using the same data
from the training and validation sets described above, 32 observables using
polynomial basis functions are created for the Koopman observables and two
similar observables are created for eDMDc. We train the eDMDc and DFL models
using the same measurements, $\zeta$, as L3.

After training, we simulate all models given a zero initial condition and a
square wave input trajectory. The modeled state trajectory is compared against
the ground truth in Fig. 5. The learned lifting linearization model
outperforms the Koopman model despite the significantly lower dimensionality.
The integrated squared errors of the simulated models are recorded in Table I.

Note that the fidelity of the Koopman operator model is sensitive to
hyperparameters. As discussed in [9], without L1 regularization, high-
dimensional Koopman models quickly overfit to the training data. Both DFL and
L3 outperform eDMDc due to the anticausal filter compensating for the
dependence of $\mathrm{e_{R}}$ on $u$. Without the anticausal filter, L3
performs only marginally better than eDMDc. L3 also has a slight advantage
over DFL: in addition to penalizing nonlinearities in the state transition
equation, the cost function of the learned lifting linearization in (12) also
penalizes nonlinearities in the augmented state transition equation.

###  IV-B Excavation Process Modeling

Figure 6: Diagram of the states and physical observables included in the data.
Figure 7: Results of open-loop simulation predicting the positional state and
observable trajectories of the bucket through the soil given the initial
condition and control input. The solid black line indicates the true, observed
trajectory. The dashed green line indicates the trajectory predicted by the
Koopman model. The dashed lines indicate the open-loop simulated trajectories
of various models: dynamic mode decomposition with control (DMDc), learned
lifting linearization (L3), L3 without first filtering out the control input
using (16) (L3 NoF), and L3 without any information from $\zeta$ (L3 NoZ).

We also test the learned lifting linearization algorithm on an autonomous
excavator simulation using agxTerrain, a specialized module of the AGX
Dynamics [39] physics simulator used to test algorithms for autonomous
excavation [40, 41, 42, 43] illustrated in Fig. 1. We generate a random soil
profile by summing several 2-D Gaussians of random height and variance,
yielding soil shapes like those in Fig. 1. Soil properties are set in
accordance to the AGX “gravel” profile.

We collect 100 7.5s randomized trajectories from large sections of the
workspace at 100Hz, setting one of them aside for testing. A diagram of the
collected data is illustrated in Fig. 6. The trajectories include six states,
$x$: position along the x-axis, $\mathrm{x}$; position along the $y$-axis,
$y$; bucket angle, $\phi$; velocity along the x-axis, $v_{\mathrm{x}}$;
velocity along the $y$-axis, $v_{y}$; and rotational velocity of the bucket,
$\omega$. The trajectories also include three control inputs, $u$: force along
the x-axis, $u_{\mathrm{x}}$; force along the $y$-axis, $u_{y}$; and torque
actuating the bucket, $u_{\phi}$. The trajectories also include three physical
observables, $\zeta$: soil reaction force on the bucket along the x-axis,
$F_{\mathrm{x}}$; soil reaction force on the bucket along the $y$-axis,
$F_{y}$; and mass of the soil in the bucket, $m_{\mathrm{soil}}$.

These trajectories are generated using a naïve, noisy PID controller on the
translation forces and bucket angle:

$\displaystyle u_{\mathrm{x}}$
$\displaystyle=\mathrm{PID}(\dot{\mathrm{x}}-\mathcal{U}(\dot{\mathrm{x}}_{\mathrm{min}},\dot{\mathrm{x}}_{\mathrm{max}}))+\mathcal{U}(-w_{\mathrm{x}},w_{\mathrm{x}})$
(17) $\displaystyle u_{y}$
$\displaystyle=\mathrm{PID}(\dot{y}-\mathcal{U}(\dot{y}_{\mathrm{min}},\dot{y}_{\mathrm{max}}))+\mathcal{U}(-w_{y},w_{y})$
$\displaystyle u_{\phi}$
$\displaystyle=\mathrm{PID}(\phi-\mathcal{U}(\phi_{\mathrm{min}},\phi_{\mathrm{max}}))+\mathcal{U}(-w_{\phi},w_{\phi})$

where $\mathcal{U}$ is the uniform random distribution and $w_{\mathrm{x}}$,
$w_{y}$, and $w_{\phi}$ are bounds on additional noise added to ensure
persistent excitation. The set points are drawn from a uniform random
distribution in accordance with [44].

The learned lifting linearization and Koopman models for the terramechanics
experiment are almost identical to those of the nonlinear spring-damper
experiment, with some exceptions. The neural network, $g$, has 9 input
neurons, 1 hidden layer with 256 ReLU neurons, and 4 output neurons. On
average, training the L3 model took 2.5 hours. The domain of the Koopman
dynamic model has a dimensionality of 67, compared to 16 for the learned
lifting linearization model. We also benchmark against a DMDc model with a
dimensionality of 12 trained using the observables, $\zeta$, in addition to
the state. There was not a significant performance difference between DMDc and
eDMDc. In addition to benchmarking the learned lifting linearization model
against the Koopman and DMDc models, we also compare the results with and
without filtering out the control input from the physical observables using
(16). In this model, instead of following the procedure described in § III-B,
we incorporated $\zeta(x,u)$ directly into the state without filtering.
Effectively, this forces the model to violate causality by predicting future
values of observables dependent upon control input, $u$. Finally, we also test
L3 without any information from observations, $\zeta$, to examine the effect
of removing supplementary measurements.

After training the learned lifting linearization and Koopman models, we
simulate both models using the control input trajectory and initial condition
from the testing trajectory. The modeled trajectories of the positional
states, x, $y$, and $\omega$, and the three observables, $F_{x}$, $F_{y}$, and
$m_{\mathrm{soil}}$, are compared against the ground truth trajectories in
Fig. 7. The integrated squared error across the three states for the learned
lifting linearization model is only 5% and 18% of the same errors for the
Koopman and DMDc models, respectively.

We also retrain the learned lifting linearization model without first
filtering out the control input using (16) and simulate using the same
technique, L3 (NoF). The integrated squared error across the three states is
more than eight times the same error for the standard L3 model. Input
filtering is vital to model robustness. A similar experiment performed without
any information from supplemental measurements, L3 (NoZ), results in an error
more than ten times greater than the standard L3 model. Various experiments
retraining the L3 model without access to individual states or observables
(left out of Fig. 7) result in similar reductions in accuracy.

##  V Discussion and Conclusion

In this paper, we presented a learned lifting linearization algorithm to model
nonlinear dynamic systems. This model extended Koopman operator theory and
dual-faceted linearization by training a neural network to produce nonlinear
observables to augment the state. We also presented an algorithm to “clean”
anticausal physical observables of any linear dependence on control input so
that they can be used by the neural network to generate richer synthetic
observables. We tested this algorithm on a nonlinear, massless spring-damper
model and an autonomous excavation simulation, and we compared the results
against Koopman, DMD, and DFL models. Learned lifting linearization
outperformed all benchmarks at minimizing state prediction error.

As with many data-driven techniques, training data quality is paramount to
model accuracy. The reduced performance in modeling progression along the
$y$-axis is likely due to the small domain of the training data along that
dimension. If we initialize the system in a different configuration from what
was recorded during data collection, or if the quantity of data available for
learning is reduced by more than half, models such as L3 and DMD perform
substantially worse.

In real-world excavation tasks, soil properties can vary. While this paper
addresses more homogeneous soil profiles like gravel, intelligent use of
observables to include information about the environment could enable L3 to
learn optimal auxiliary variables and linear models for dynamic soil
properties. Additionally, L3 could be primed with simulator data to reduce the
amount of hardware-in-the-loop learning required.

Excavators often require expensive hydraulics modifications to execute end
effector force control. State-of-the-art excavators have solved this problem
[7]. This work focused on the nonlinearities involved in soil dynamics. Future
work on this topic should model the nonlinearity of hydraulic systems.

## References

  * [1] S. Dadhich, U. Bodin, and U. Andersson, “Key challenges in automation of earth-moving machines,” Automation in Construction, vol. 68, pp. 212–222, 2016.
  * [2] J. Hook, “A global forecast for the construction industry to 2030,” tech. rep., Global Construction Perspectives and Oxford Economics, Broadwall House, 21 Broadwall, London SE1 9PL United Kingdom, Sept. 2015.
  * [3] F. E. Sotiropoulos and H. H. Asada, “A model-free extremum-seeking approach to autonomous excavator control based on output power maximization,” IEEE Robotics and Automation Letters, vol. 4, no. 2, pp. 1005–1012, 2019.
  * [4] R. Fukui, T. Niho, M. Nakao, and M. Uetake, “Imitation-based control of automated ore excavator: improvement of autonomous excavation database quality using clustering and association analysis processes,” Advanced Robotics, vol. 31, no. 11, pp. 595–606, 2017.
  * [5] O. Luengo, S. Singh, and H. Cannon, “Modeling and identification of soil-tool interaction in automated excavation,” in Proceedings. 1998 IEEE/RSJ International Conference on Intelligent Robots and Systems. Innovations in Theory, Practice and Applications (Cat. No. 98CH36190), vol. 3, pp. 1900–1906, IEEE, 1998.
  * [6] R. J. Sandzimier and H. H. Asada, “A data-driven approach to prediction and optimal bucket-filling control for autonomous excavators,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 2682–2689, 2020.
  * [7] D. Jud, P. Leemann, S. Kerscher, and M. Hutter, “Autonomous free-form trenching using a walking excavator,” IEEE Robotics and Automation Letters, vol. 4, no. 4, pp. 3208–3215, 2019.
  * [8] I. Abraham and T. D. Murphey, “Active learning of dynamics for data-driven control using koopman operators,” IEEE Transactions on Robotics, vol. 35, no. 5, pp. 1071–1083, 2019.
  * [9] D. Bruder, B. Gillespie, C. D. Remy, and R. Vasudevan, “Modeling and control of soft robots using the koopman operator and model predictive control,” in Proceedings of Robotics: Science and Systems, (FreiburgimBreisgau, Germany), June 2019.
  * [10] A. Broad, T. Murphey, and B. Argall, “Learning models for shared control of human-machine systems with unknown dynamics,” in Proceedings of Robotics: Science and Systems, (Cambridge, Massachusetts), July 2017.
  * [11] Y. Susuki, I. Mezic, F. Raak, and T. Hikihara, “Applied koopman operator theory for power systems technology,” Nonlinear Theory and Its Applications, IEICE, vol. 7, no. 4, pp. 430–459, 2016.
  * [12] A. Leonard, J. Rogers, and A. Gerlach, “Koopman operator approach to airdrop mission planning under uncertainty,” Journal of Guidance, Control, and Dynamics, vol. 42, no. 11, pp. 2382–2398, 2019.
  * [13] M. Zhang, S. Vikram, L. Smith, P. Abbeel, M. Johnson, and S. Levine, “Solar: Deep structured representations for model-based reinforcement learning,” in International Conference on Machine Learning, pp. 7444–7453, PMLR, 2019\.
  * [14] B. Ichter and M. Pavone, “Robot motion planning in learned latent spaces,” IEEE Robotics and Automation Letters, vol. 4, no. 3, pp. 2407–2414, 2019\.
  * [15] T. Haarnoja, A. Ajay, S. Levine, and P. Abbeel, “Backprop kf: Learning discriminative deterministic state estimators,” in Advances in Neural Information Processing Systems (D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, eds.), vol. 29, Curran Associates, Inc., 2016.
  * [16] M. Watter, J. Springenberg, J. Boedecker, and M. Riedmiller, “Embed to control: A locally linear latent dynamics model for control from raw images,” in Advances in Neural Information Processing Systems (C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, eds.), vol. 28, Curran Associates, Inc., 2015.
  * [17] A. Mauroy, Y. Susuki, and I. Mezić, The Koopman Operator in Systems and Control.  Springer, 2020.
  * [18] S. Vijayshankar, S. Nabi, A. Chakrabarty, P. Grover, and M. Benosman, “Dynamic mode decomposition and robust estimation: Case study of a 2d turbulent boussinesq flow,” in 2020 American Control Conference (ACC), pp. 2351–2356, IEEE, 2020.
  * [19] H. Harry Asada and F. E. Sotiropoulos, “Dual faceted linearization of nonlinear dynamical systems based on physical modeling theory,” Journal of Dynamic Systems, Measurement, and Control, vol. 141, no. 2, 2019.
  * [20] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” Automatica, vol. 93, pp. 149–160, 2018.
  * [21] B. O. Koopman, “Hamiltonian systems and transformation in hilbert space,” Proceedings of the national academy of sciences of the united states of america, vol. 17, no. 5, pp. 315–318, 1931.
  * [22] S. Maćešić and N. Črnjarić-Žic, “Koopman operator theory for nonautonomous and stochastic systems,” in The Koopman Operator in Systems and Control, pp. 131–160, Springer, 2020.
  * [23] I. Mezić, “Analysis of fluid flows via spectral properties of the koopman operator,” Annual Review of Fluid Mechanics, vol. 45, pp. 357–378, 2013\.
  * [24] M. Budišić, R. Mohr, and I. Mezić, “Applied koopmanism,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 22, no. 4, p. 047510, 2012.
  * [25] J. H. Tu, C. W. Rowley, S. L. B. Dirk M. Luchtenburg, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” Journal of Computational Dynamics, vol. 1, no. 2, pp. 391–421, 2014.
  * [26] D. C. Karnopp, D. L. Margolis, and R. C. Rosenberg, System dynamics: modeling, simulation, and control of mechatronic systems.  John Wiley & Sons, 2012.
  * [27] Y. Igarashi, M. Yamakita, J. Ng, and H. H. Asada, “Mpc performances for nonlinear systems using several linearization models,” in 2020 American Control Conference (ACC), pp. 2426–2431, IEEE, 2020.
  * [28] F. E. Sotiropoulos and H. H. Asada, “Causality in dual faceted linearization of nonlinear dynamical systems,” in 2018 Annual American Control Conference (ACC), pp. 1230–1237, IEEE, 2018.
  * [29] Y. Han, W. Hao, and U. Vaidya, “Deep learning of koopman representation for control,” in 2020 59th IEEE Conference on Decision and Control (CDC), pp. 1890–1895, IEEE, 2020.
  * [30] B. Lusch, J. N. Kutz, and S. L. Brunton, “Deep learning for universal linear embeddings of nonlinear dynamics,” Nature communications, vol. 9, no. 1, pp. 1–10, 2018.
  * [31] D. Masti and A. Bemporad, “Learning nonlinear state–space models using autoencoders,” Automatica, vol. 129, p. 109666, 2021.
  * [32] E. Yeung, S. Kundu, and N. Hodas, “Learning deep neural network representations for koopman operators of nonlinear dynamical systems,” in 2019 American Control Conference (ACC), pp. 4832–4839, IEEE, 2019.
  * [33] G. Mamakoukas, M. Castano, X. Tan, and T. Murphey, “Local koopman operators for data-driven control of robotic systems,” in Proceedings of Robotics: Science and Systems, (FreiburgimBreisgau, Germany), June 2019.
  * [34] T. Kailath, Linear systems, vol. 156.  Prentice-Hall Englewood Cliffs, NJ, 1980.
  * [35] P. M. DeRusso, A. A. Desrochers, R. J. Roy, and C. M. Close, State variables for engineers.  John Wiley & Sons, Inc., 1997.
  * [36] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems 32 (H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, eds.), pp. 8024–8035, Curran Associates, Inc., 2019.
  * [37] N. S. Selby, “Dfl.” https://github.com/Darbeloff/DFL, 2021.
  * [38] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (Y. Bengio and Y. LeCun, eds.), 2015.
  * [39] M. Servin, T. Berglund, and S. Nystedt, “A multiscale model of terrain dynamics for real-time earthmoving simulation,” Advanced Modeling and Simulation in Engineering Sciences, vol. 8, pp. 11–46, May 2021.
  * [40] Y. Yang, J. Pan, P. Long, X. Song, and L. Zhang, “Time variable minimum torque trajectory optimization for autonomous excavator,” arXiv preprint arXiv:2006.00811, 2020.
  * [41] Y. Yang, L. Zhang, X. Cheng, J. Pan, and R. Yang, “Compact reachability map for excavator motion planning,” in 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2308–2313, IEEE, 2019\.
  * [42] S. Ulin, “Digging deep: A data-driven approach to model reduction in a granular bulldozing scenario,” 2018.
  * [43] S. Backman, D. Lindmark, K. Bodin, M. Servin, J. Mörk, and H. Löfgren, “Continuous control of an underground loader using deep reinforcement learning,” 2021.
  * [44] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Generalizing koopman theory to allow for inputs and control,” SIAM Journal on Applied Dynamical Systems, vol. 17, no. 1, pp. 909–930, 2018.