††thanks: Present address: Samsung Electronics, Gyeonggi–do 16677, Republic of
Korea

# Characterization of a flux-driven Josephson parametric amplifier with near
quantum-limited added noise for axion search experiments

Çağlar Kutlu<EMAIL_ADDRESS>Korea Advanced Institute of Science and
Technology, Daejeon 34051, Republic of Korea  Center for Axion and Precision
Physics Research, Institute for Basic Science, Daejeon 34051, Republic of
Korea  Arjan F. van Loo  Center for Emergent Matter Science (CEMS), RIKEN,
Wako, Saitama 351–0198, Japan  Sergey V. Uchaikin  Andrei N. Matlashov  Center
for Axion and Precision Physics Research, Institute for Basic Science, Daejeon
34051, Republic of Korea  Doyu Lee  Center for Axion and Precision Physics
Research, Institute for Basic Science, Daejeon 34051, Republic of Korea
Seonjeong Oh  Jinsu Kim  Korea Advanced Institute of Science and Technology,
Daejeon 34051, Republic of Korea  Center for Axion and Precision Physics
Research, Institute for Basic Science, Daejeon 34051, Republic of Korea
Woohyun Chung  Center for Axion and Precision Physics Research, Institute for
Basic Science, Daejeon 34051, Republic of Korea  Yasunobu Nakamura  Center for
Emergent Matter Science (CEMS), RIKEN, Wako, Saitama 351–0198, Japan  Research
Center for Advanced Science and Technology (RCAST), The University of Tokyo,
Meguro–ku, Tokyo 153–8904, Japan  Yannis K. Semertzidis  Korea Advanced
Institute of Science and Technology, Daejeon 34051, Republic of Korea  Center
for Axion and Precision Physics Research, Institute for Basic Science, Daejeon
34051, Republic of Korea

###### Abstract

The axion, a hypothetical elementary pseudoscalar, is expected to solve the
strong _CP_ problem of QCD and is also a promising candidate for dark matter.
The most sensitive axion search experiments operate at millikelvin
temperatures and hence rely on instrumentation that carries signals from a
system at cryogenic temperatures to room temperature instrumentation. One of
the biggest limiting factors affecting the parameter scanning speed of these
detectors is the noise added by the components in the signal detection chain.
Since the first amplifier in the chain limits the minimum noise, low-noise
amplification is of paramount importance. This paper reports on the operation
of a flux-driven Josephson parametric amplifier (JPA) operating at around
$2.3\text{\,}\mathrm{GHz}$ with added noise approaching the quantum limit. The
JPA was employed as a first stage amplifier in an experimental setting similar
to the ones used in haloscope axion detectors. By operating the JPA at a gain
of $19\text{\,}\mathrm{dB}$ and cascading it with two cryogenic amplifiers
operating at $4\text{\,}\mathrm{K}$, noise temperatures as low as
$120\text{\,}\mathrm{mK}$ were achieved for the whole signal detection chain.

##  I Introduction

Axions are spin-0 particles that emerge as a result of the Peccei-Quinn
mechanism which was originally proposed as a solution to the strong _CP_
problem of quantum chromodynamics[1, 2]. They were also identified as viable
candidates for all or a fraction of the cold dark matter in our universe[3, 4,
5]. It is possible to detect axions upon their conversion to microwave
photons, using resonant cavities immersed in high magnetic fields[6]. Since
the axion mass is unknown, these detectors employ a mechanism to scan
different frequencies corresponding to different axion masses. The scanning
rate of such detectors scales with $1/T_{\mathrm{sys}}^{2}$, where
$T_{\mathrm{sys}}$ is the system noise background characterized in units of
temperature. It can be decomposed as
$T_{\mathrm{sys}}=T_{\mathrm{cav}}+T_{\mathrm{add}}$, where the first term
denotes the noise temperature accompanying the signal itself and the second
one denotes the noise added by the signal detection chain. Throughout this
work, noise temperature refers to the added noise unless otherwise stated. In
order to reduce $T_{\mathrm{cav}}$, the cavity is cooled to millikelvin
temperatures. If the first amplifier has sufficiently high gain ($G_{1}$), its
noise temperature ($T_{1}$) will be the dominant contribution to
$T_{\mathrm{add}}$ as given by the well-known relation[7]:
$T_{\mathrm{add}}=T_{1}+\frac{T_{\mathrm{rest}}}{G_{1}}$ where
$T_{\mathrm{rest}}$ is the noise temperature of the whole chain except the
first amplifier. Amplifiers based on Josephson junctions including microstrip
superconducting quantum interference device amplifiers (MSA) and JPA have
already been shown to be capable of gains higher than
$30\text{\,}\mathrm{dB}$, and noise temperatures approaching the quantum
limit[8, 9]. While an MSA has an internal shunt resistor used for biasing
which hinders noise performance[10, 11], by design the JPA requires no
resistive element to operate. Several experiments presently searching for dark
matter axions have already adopted the JPA as the first amplifier[12, 13, 14].
In this work, the frequency coverage, gain and noise properties of a flux-
driven JPA for use in an axion dark matter experiment operating around
$2.3\text{\,}\mathrm{GHz}$ are investigated.

The power spectral density of the noise accompanying a signal measured in an
impedance matched environment can be given as [15] :

$\displaystyle
S_{n}(f,T)=hf\left[\frac{1}{\exp{\left(\frac{hf}{k_{B}T}\right)}-1}+\frac{1}{2}\right]$
(1)

where $h$ is Planck’s constant and $k_{B}$ is Boltzmann’s constant. The first
term in the brackets is the mean number of quanta at frequency $f$ at the bath
temperature $T$ and the second term is the contribution from zero-point
fluctuations. The lower limit on noise temperature for linear phase-
insensitive amplifiers is given by[16] $T_{Q}=\lim_{T\rightarrow
0}S_{n}(f,T)/(k_{B})=hf/(2k_{B})$ which is about $55.2\text{\,}\mathrm{mK}$ at
$2.3\text{\,}\mathrm{GHz}$. Using a $2.3\text{\,}\mathrm{GHz}$ flux-driven JPA
$T_{\mathrm{add}}\approx$120\text{\,}\mathrm{mK}$$ is achieved. This
corresponds to a $T_{\mathrm{sys}}\approx$190\text{\,}\mathrm{mK}$$ for an
axion haloscope experiment running at a bath temperature of
$50\text{\,}\mathrm{mK}$. The lower bound for $T_{\mathrm{sys}}$ is given by
the standard quantum limit[17] $T_{\mathrm{SQL}}=2T_{Q}$ which is about
$110\text{\,}\mathrm{mK}$ at $2.3\text{\,}\mathrm{GHz}$.

##  II Flux-driven JPA

The equivalent circuit diagram of the tested device is shown in Figure 1. It
consists of a superconducting quantum interference device (SQUID) attached to
the end of a coplanar waveguide $\lambda/4$ resonator that is coupled via a
capacitor ($C_{c}$) to the transmission line for the signal input and output.
The SQUID acts as a variable inductor whose value depends on the magnetic flux
passing through its loop. In the setup, a superconducting coil is used to
provide the necessary DC flux ($\phi$) through the SQUID loop in order to tune
the resonance frequency ($f_{r}$). Parametric amplification is achieved by
modulating the flux through the SQUID using a pump signal. The pump tone is
provided by a separate transmission line inductively coupled to the SQUID. The
JPA is operated in the three-wave mixing mode[18] where the pump ($f_{p}$),
the signal ($f_{s}$), and the idler ($f_{i}$) frequencies satisfy the relation
$f_{p}=f_{s}+f_{i}$. The signal input and output share the same port. A
circulator is used to separate them. Since the $\lambda/4$ resonator only
allows odd harmonics, there is no measurable pump leakage to the output line.
This prevents the stronger pump tone from saturating the rest of the
amplifiers in the chain[19]. Figure 2 shows a schematic for the axion search
experimental setup.

Figure 1: Equivalent circuit diagram of the JPA sample. The JPA was fabricated
by photolithography of a Nb layer, deposited on a $0.3\text{\,}\mathrm{mm}$
thick Si substrate. The SQUID was placed on top of the Nb layer by E-beam
lithography followed by shadow evaporation[20, 21]. The sample was attached to
a printed circuit board (PCB) and the transmission lines were bonded with Al
wires. The PCB was fixed onto a gold plated copper structure and placed inside
a superconducting coil. The whole structure was covered tightly with a lead
shield and attached to the mixing-chamber (MC) plate using a gold plated
copper rod.

##  III Measurements

When there is no pump tone present, the JPA can be modeled as a resonator with
a well-defined quality factor and resonance frequency which are functions of
flux. The resonance frequency is estimated from the frequency domain phase
response using a parameter fit[22]. The phase response is obtained by doing a
transmission S-parameter measurement using a vector network analyzer (VNA) in
the configuration as shown in Figure 2. The resonance frequency was measured
as a function of the coil current (see Figure 3). It was found that the
minimum observable resonance frequency was at $2.18\text{\,}\mathrm{GHz}$ and
the maximum was $2.309\text{\,}\mathrm{GHz}$. The lower bound is due to the
frequency band of the circulators which spans from $2.15\text{\,}\mathrm{GHz}$
to $2.60\text{\,}\mathrm{GHz}$. At the lower frequencies, the JPA becomes much
more sensitive to flux noise due to a higher $\frac{\partial
f_{r}}{\partial\phi}$. This work mainly focused on operation with frequencies
above $2.2\text{\,}\mathrm{GHz}$.

Figure 2: The experimental setup used in all the characterization
measurements. SG, VNA, and SA stand for the signal generator, vector network
analyzer, and spectrum analyzer, respectively. During this work, the switch
that selects between the cavity and the noise source was always kept at the
position shown in the figure. The ports IN2 and OUT were used to directly
measure the JPA characteristics, bypassing the cavity. The microwave short
element shown next to the JPA was used to bypass the JPA for calibration
measurements. U1 and U2 are HEMT amplifiers with noise temperatures of
$1.5\text{\,}\mathrm{K}$ and $5\text{\,}\mathrm{K}$, respectively.

During the experiments, the MC plate temperature was stabilized at
$50\pm{}1\text{\,}\mathrm{mK}$. With the temperature fixed, the frequency
response of the JPA is determined by three experimental variables: the coil
current ($i_{b}$), the pump frequency ($f_{p}$), and the pump power ($P_{p}$).
The measurements shown in this work had $i_{b}$ confined to the region where
the flux through the SQUID loop is given by $-0.5\phi_{0}<\phi<0$, where
$\phi_{0}$ is the magnetic flux quantum. Therefore, $f_{r}$ can be
unambiguously converted to $\phi$ or $i_{b}$. All experiments began with a
transmission measurement, with the resonance frequency tuned to
$2.18\text{\,}\mathrm{GHz}$. This becomes the baseline measurement to be used
for the duration of the experiment. When the result was compared to a separate
measurement, in which a microwave short was put in place of the JPA, it was
found that the baseline obtained via such an off-resonance measurement was at
most $0.2\text{\,}\mathrm{dB}$ lossier than an ideal mirror. The JPA gain
($G_{J}$) was estimated by dividing the transmission magnitude response with
the baseline’s magnitude response.

To investigate the gain behavior, a sweep over the parameters $i_{b}$,
$f_{p}$, $P_{p}$ was made and the maximum gain was measured at each point.
After each $i_{b}$ tuning step, the resonance frequency is estimated by
performing a phase measurement and applying a parameter fit. With the detuning
defined as $\delta=f_{p}/2-f_{r}$, the equigain contours had a minimum in
necessary pump power around $\delta=0$, as shown in Figure 4(a). It was
observed that for resonance frequencies above $2.299\text{\,}\mathrm{GHz}$ the
minimum starts to shift to lower detunings which is attributed to pump-induced
shifts in resonance frequency[22]. Figure 4(b) shows that the slice of
$\delta=0$ can be used to achieve peak gains of up to $30\text{\,}\mathrm{dB}$
along the frequency range of the device.

Figure 3: Resonance frequency versus flux obtained by sweeping the coil
current and measuring the phase response at each step. One period corresponds
to a current of $324.4\text{\,}\mathrm{\SIUnitSymbolMicro A}$. The inset shows
the fit performed to estimate the resonance frequency for each applied flux.

(a)

(b)

Figure 4: (a) Maximum gain measured as a function of detuning and pump power
for a flux bias corresponding to $f_{r}=$2.29\text{\,}\mathrm{GHz}$$. (b)
Maximum gain as a function of frequency and pump power with $f_{p}=2f_{r}$.

To investigate noise temperature, a methodology similar to the well-known
Y-factor method[23] was used. A $50\text{\,}\mathrm{\SIUnitSymbolOhm}$
cryogenic microwave terminator was used as the noise source. A bias-tee was
attached in front for improved thermalization of its inner conductor. These
two components were fixed onto a gold-plated copper plate along with a
ruthenium oxide temperature sensor and a
$100\text{\,}\mathrm{\SIUnitSymbolOhm}$ resistor functioning as a heater. This
plate was then fixed onto the MC plate so that the dominant thermalization was
through a thin copper wire attached to the MC plate. The noise source was
connected to the switch input using a superconducting coaxial cable, which
provides thermal isolation while minimizing losses. Using a PID controller,
the terminator temperature could be adjusted from $50\text{\,}\mathrm{mK}$ to
$1\text{\,}\mathrm{K}$ without affecting the MC plate temperature. The noise
power generated by the noise source was measured using a spectrum analyzer
(SA) with $1\text{\,}\mathrm{kHz}$ resolution bandwidth after being amplified
by the JPA and the rest of the signal detection chain. The power spectra were
recorded at noise source temperatures ($T_{s}$) of
$60120180\text{\,}\mathrm{mK}$. The power values were converted into power
spectral densities (PSD) by dividing them with the noise bandwidth
corresponding to the SA settings used. Before each PSD measurement, the JPA
gain and passive resonance were measured. From these measurements, it was
concluded that there were neither gain changes nor resonance shifts. From the
obtained PSD values $S(T_{s})$, a fit was done to a function of the following
form independently for each frequency bin (see Figure 5) :

$\displaystyle
S(T_{s})=(2G_{J}-1)\frac{G_{L}G_{\mathrm{tot}}}{G_{J}}(S_{n}(T_{s})+rk_{B}T_{n}+\gamma)$
(2)

where $S_{n}(T_{s})$ is the noise PSD of the source, $G_{\mathrm{tot}}$ is the
total gain seen from the reference plane, $G_{L}$ is the loss factor between
the $50\text{\,}\mathrm{\SIUnitSymbolOhm}$ terminator and the reference plane
and $T_{n}$ is the noise temperature. The reference plane is at the end of the
superconducting cable connected to the noise source (see Figure 2). Here, $r$
and $\gamma$ are factors that are explained in Appendix A.

Figure 5: The upper plot shows the set of power spectra obtained during a
noise temperature measurement performed for a tuning at
$f_{r}=$2.305\text{\,}\mathrm{GHz}$$ with $f_{p}=2f_{r}$. The offset $\nu$ is
defined as $\nu=f-f_{r}$ where $f$ is the center of the frequency bin at which
the power was measured using the spectrum analyzer. $T_{s}$ is the temperature
of the noise source. The inset shows three vertical slices which were fit with
Equation 2. The lower plot shows the estimated noise temperature of the whole
chain as a function of $\nu$.

Since the amplifier needs to be tuned along with the cavity during the axion
experiment, the noise temperature was investigated at different frequencies.
The measurements were done in $5\text{\,}\mathrm{MHz}$ steps from
$2.282.305\text{\,}\mathrm{GHz}$. At each step, the pump power and resonance
frequency were tuned such that the JPA gain was about $20$ $\mathrm{dB}$. From
these measurements (Figure 6) a minimum noise temperature of
$120\text{\,}\mathrm{mK}$ was observed at $2.28\text{\,}\mathrm{GHz}$.

Figure 6: Total gain and the noise temperature of the whole chain for 6 tuning
points with $19.3\pm{}0.5\text{\,}\mathrm{dB}$ JPA gain. Both quantities were
estimated from noise temperature measurements. The small variations along
tuning frequencies were mainly attributed to the losses due to the microwave
components before the JPA.

Another important characteristic is the saturation that occurs when a
narrowband signal is applied. For stable and predictable operation the JPA
must be operated away from the effects of saturation. A common way to quantify
the saturation of an amplifier is to determine the input power at which the
gain is reduced by $1\text{\,}\mathrm{dB}$ ($P_{1\mathrm{dB}}$). The
$P_{1\mathrm{dB}}$ was measured at $\delta=0$ for different frequencies and
different pump powers corresponding to different gains. It is evident from the
results (see Figure 7) that an axion-like signal with an expected power of
$-180\text{\,}\mathrm{\milli}$ is far from saturating the device. While
saturation from narrowband signals is avoidable to a certain extent, it was
observed that thermal noise at the input can also saturate and alter the
behavior of the device. For frequencies below $2.28\text{\,}\mathrm{GHz}$ with
gains above $23\text{\,}\mathrm{dB}$, the device started showing saturated
behavior with thermal noise when the noise source temperature was raised above
$120\text{\,}\mathrm{mK}$, which was done to measure noise temperature. While
this does not necessarily mean that the device is unusable below these
frequencies, it renders the direct measurement of the noise temperature using
a noise source unreliable for these frequency and gain regions.

Figure 7: Saturation measurements for three different $f_{r}$. Each
measurement was done by sweeping the signal powers from the VNA and observing
at which input power the gain reduces by 1 dB. The horizontal axis corresponds
to the unsaturated gain measured with the lowest signal power available from
VNA.

##  IV Conclusion

In conclusion, a flux-driven JPA, tunable in the range
$2.22.305\text{\,}\mathrm{GHz}$ was demonstrated and determined to be
operational for use in axion search experiments. The added noise temperatures
of the receiver chain were measured using a noise source at a location as
close as possible to the origin of the axion signal. With an added noise
temperature of $120\text{\,}\mathrm{mK}$ the system was shown to reach
$T_{\mathrm{sys}}\approx 1.7T_{\mathrm{SQL}}$. This is the first record of
$T_{\mathrm{sys}}$ below $2T_{\mathrm{SQL}}$ for an axion haloscope setup
operating below $10\text{\,}\mathrm{GHz}$. The saturation input power for the
JPA was observed to be more than adequate for an axion-like signal. Currently,
the tested JPA is being used as part of a KSVZ[24, 25] sensitive axion search
experiment at the Center for Axion and Precision Physics Research (CAPP). The
system is taking physics data with a scanning speed that has been improved
more than an order of magnitude. We expect that further optimization of the
JPA design could result in improved instantaneous bandwidth and tuning range.

This work was supported by the Institute for Basic Science
(IBS-R017–D1–2021–a00) and JST ERATO (Grant No. JPMJER1601). A. F. van Loo is
supported by a JSPS postdoctoral fellowship.

##  Appendix A Noise Temperature Estimation

Figure 8: The simplified model used for the noise temperature estimations
conducted in this work. Bold letters denote the power gains of components. The
reference plane marks the input of the detector chain. Arrows denote the flow
of power entering the nodes shown with a small circle. $G_{L}$ is a composite
gain factor for everything between the noise source and the reference plane.
$G_{c}$ is the circulator gain factor. $G_{J}$ is the signal and $G_{I}$ is
the idler gain for the JPA. The amplifier gain $G_{R}$ and noise temperature
$T_{nR}$ contain the effects of all elements after the last circulator,
including SA noise. For simplicity, circulators are assumed to have complete
rotational symmetry with respect to their ports and to be completely identical
to each other.

The output PSD from a component with its input connected to a matched source
can be written as :

$\displaystyle S_{\mathrm{O}}=GS_{\mathrm{in}}+S_{\mathrm{added}}$ (3)

where $S_{\mathrm{in}}$ is the source PSD, $G$ is the power gain of the
component, and $S_{\mathrm{added}}$ is the noise added by it. The noise
temperature ($T_{n}$) is a measure of the added noise at the output of a
component. By convention, it is defined as if it is for noise entering the
device itself: $T_{n}=S_{\mathrm{added}}/(k_{B}G)$. The entire detection chain
(see Figure 8), from the reference plane to the spectrum analyzer, can be
described as a single composite component with $G=G_{\mathrm{tot}}$ and noise
temperature $T_{n}$.

The noise temperature can be defined for a situation similar to the
experimental one where a narrowband axion signal with power $A$ is present.
This signal enters the chain from a source connected to the reference plane.
Assuming the source is thermalized to the MC plate with the temperature
$T_{f}$, then the defining relation for $T_{n}$ can be written as :

$\displaystyle
S_{O}=G_{\mathrm{tot}}(\underbrace{A\delta(f-f_{s})+S_{n}(T_{f})}_{S_{\mathrm{in}}}+k_{B}T_{n})$
(4)

where $S_{O}$ is the PSD at the output, $G_{\mathrm{tot}}$ is the total power
gain for the signal from the reference plane, $S_{n}$ is the noise coming from
the source itself. The main idea here is that if one has a reliable estimate
of $T_{n}$, and understands the source environment well ($S_{n}(T_{f})$), it
is straightforward to estimate $A$ without the precise knowledge of
$G_{\mathrm{tot}}$. This is possible since $S_{O}$ can be easily measured at
two frequencies $f_{s}$ and $f_{s}^{\prime}$ using a spectrum analyzer.
Provided that $|f_{s}-f_{s}^{\prime}|$ is small enough so that $T_{n}$ is
approximately the same for both frequencies, $A$ can be estimated from these
two measurements. This approach forms the basis of the analysis methods
applied in axion dark matter search experiments[26, 27, 28, 29].

The detection chain consists of passive components, the JPA and the HEMT
amplifiers. Each one of these adds noise in a different way. A passive
component at physical temperature $T_{f}$ has
$S_{\mathrm{added}}=(1-G)S_{n}(T_{f})$. The HEMT amplifier noise is usually
estimated from measurements. The JPA adds noise by two main mechanisms. The
first one is by amplifying the input noise at the idler mode onto the signal
mode. The second one is via the losses or other dissipation mechanisms inside
or before the sample. Ideally, the latter can be made zero, whereas the former
will approach to the half-photon added noise in the limit of a
$0\text{\,}\mathrm{K}$ bath temperature.

Using the model shown in Figure 8, it is straightforward to write a relation
for the output PSD. For clarity, the explicit frequency dependence of the
thermal noise $S_{n}$ and of the gains will be omitted. Also, the
approximation $S_{n}(f,T)\approx S_{n}(f_{p}-f,T)$ will be denoted with the
shorthand $S_{nf}=S_{n}(f,T_{f})$. This approximation has less than 30 ppm
error given that $|2f-f_{p}|<$100\text{\,}\mathrm{kHz}$$. Note that
$100\text{\,}\mathrm{kHz}$ is the typical bandwidth for the JPA tested in this
work. Furthermore, the transmission characteristics of the microwave
components will be assumed to not vary on a scale of
$100\text{\,}\mathrm{kHz}$. Using the gain symbols for components as shown in
Figure 8, the power flow at each node in terms of their PSD is written as :

$\displaystyle\begin{split}S_{B}&=G_{L}S_{A}+(1-G_{L})S_{nf}\\\
S_{C}&=G_{c}S_{B}+(1-G_{c})S_{nf}\\\ S_{D}&=G_{c}S_{C}+(1-G_{c})S_{nf}\\\
S_{E}&=G_{J}S_{D}+G_{I}S_{D}+G_{J}S_{j}\\\
S_{F}&=G_{c}S_{E}+(1-G_{c})S_{nf}\\\
S_{O}&=G_{R}(S_{F}+k_{B}T_{nR})\end{split}$ (5)

The idler gain is denoted by $G_{I}$, and is substituted using
$G_{I}=G_{J}-1$[30] in the following derivations. As shown in Equation 5, the
idler contribution to the noise appears as $G_{I}S_{D}$. The symbol $S_{j}$
denotes an unknown noise density added at the JPA stage which does not
contribute to the quantum limit but rather contains losses or other mechanisms
of stationary noise. Note that the noise propagating back from the later
stages is also included in $S_{j}$. The output $S_{O}$ can be written for two
cases. In the first case, the noise source is operational at temperature
$T_{s}$, and in the second case, a signal source at temperature $T_{f}$ is
connected to the reference plane. The former case describes the measurement
situation, whereas the latter case is only used to define $T_{n}$ in terms of
the parameters in the model. For the first case, i.e. $S_{A}=S_{n}(T_{s})$,
the output PSD can be written as :

$\displaystyle S_{O}^{(1)}$
$\displaystyle=\underbrace{G_{R}G_{c}^{3}G_{L}(2G_{J}-1)}_{G_{\mathrm{noise}}}\left[S_{n}(T_{s})+S_{\alpha}\right]$
(6) $\displaystyle S_{\alpha}$
$\displaystyle=\lambda^{(1)}S_{nf}+\frac{G_{J}S_{j}}{(2G_{J}-1)G_{c}^{2}G_{L}}+\left.\frac{k_{B}T_{nR}}{G_{c}^{3}G_{L}(2G_{J}-1)}\right.$
(7) $\displaystyle\lambda^{(1)}$
$\displaystyle=\beta_{l}+\frac{\beta_{c}}{G_{L}}+\frac{\beta_{c}}{G_{c}G_{L}}+\frac{\beta_{c}}{G_{c}^{2}G_{L}}$
(8) $\displaystyle\beta_{\mathord{\color[rgb]{0,0,0}\bullet}}$
$\displaystyle\equiv\frac{1}{G_{\mathord{\color[rgb]{0,0,0}\bullet}}}-1$ (9)

The output for the second case, where $S_{B}=A\delta(f-f_{s})+S_{nf}$, is
written as :

$\displaystyle S_{O}^{(2)}$
$\displaystyle=\underbrace{G_{J}G_{c}^{3}G_{R}}_{G_{\mathrm{tot}}}\left[S_{B}+k_{B}T_{n}\right]$
(10) $\displaystyle k_{B}T_{n}$
$\displaystyle=\lambda^{(2)}S_{nf}+\frac{S_{j}}{G_{c}^{2}}+k_{B}\frac{T_{nR}}{G_{c}^{3}G_{J}}$
(11) $\displaystyle\lambda^{(2)}$
$\displaystyle=\left[\frac{G_{J}-1}{G_{J}}+\left(\beta_{c}+\frac{\beta_{c}}{G_{c}}\right)\frac{2G_{J}-1}{G_{J}}+\frac{\beta_{c}}{G_{c}^{2}}\right]$
(12)

Here, the unknowns are $S_{j}$, $T_{nR}$ and $G_{R}$. It is clear from
Equations 11 and 12 that $T_{n}$ approaches $T_{Q}$ as expected in the limits
of $G_{J}\gg 1$, $G_{c}\rightarrow 0$, and $G_{L}\rightarrow 0$. Using
Equations 6, 7 and 11 $,S_{O}^{(1)}$ can be rewritten as :

$\displaystyle\begin{split}S_{O}^{(1)}&=\frac{G_{\mathrm{tot}}}{r}(S_{n}(T_{s})+rk_{B}T_{n}+\gamma)\\\
r&=\frac{G_{J}}{G_{L}(2G_{J}-1)}\\\
\gamma&=\left(\lambda^{(1)}-r\lambda^{(2)}\right)S_{nf}\end{split}$ (13)

This relation is used to perform a fit with $G_{\mathrm{tot}}$ and $T_{n}$ as
the fit parameters. For the estimations, the parameters $G_{L}$ and $G_{c}$
were taken as $-0.05\text{\,}\mathrm{dB}$, and $-0.4\text{\,}\mathrm{dB}$
respectively. Some typical values for $r$ and $\gamma/k_{B}$ can be found in
Table 1.

Table 1: Calculated $r$ and $\gamma/k_{B}$ for typical component losses and $T_{f}=$50\text{\,}\mathrm{mK}$$. $G_{J}($\mathrm{dB}$)$ | $G_{L}($\mathrm{dB}$)$ | $G_{c}($\mathrm{dB}$)$ | $r$ | $\gamma/k_{B}($\mathrm{mK}$)$
---|---|---|---|---
15 | 0 | -0.3 | 0.508 | -31.1
15 | -0.2 | -0.5 | 0.532 | -26.8
20 | 0 | -0.3 | 0.503 | -31.4
20 | -0.2 | -0.5 | 0.526 | -27.1
19 | -0.05 | -0.4 | 0.509 | -29.8

## References

  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Preskill, Wise, and Wilczek [1983] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Sikivie [1983] P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * Friis [1944] H. Friis, Proceedings of the IRE 32, 419 (1944).
  * Castellanos-Beltran and Lehnert [2007] M. A. Castellanos-Beltran and K. W. Lehnert, Appl. Phys. Lett. 91 (2007), 10.1063/1.2773988.
  * Kinion and Clarke [2011] D. Kinion and J. Clarke, Appl. Phys. Lett. 98, 10 (2011).
  * Uchaikin  _et al._ [2019] S. Uchaikin, A. Matlashov, D. Lee, W. Chung, S. J. Oh, Y. Semertzidis, V. Zakosarenko, Ç Kutlu, A. van Loo, Y. Urade, S. Kono, M. Schmelz, R. Stolz, and Y. Nakamura, in  _2019 IEEE International Superconductive Electronics Conference (ISEC)_ (2019) pp. 1–3.
  * André  _et al._ [1999] M. O. André, M. Mück, J. Clarke, J. Gail, and C. Heiden, Appl. Phys. Lett. 75, 698 (1999).
  * Braine  _et al._ [2020] T. Braine, R. Cervantes, N. Crisosto, N. Du, S. Kimes, L. J. Rosenberg, G. Rybka, J. Yang, D. Bowring, A. S. Chou, R. Khatiwada, A. Sonnenschein, W. Wester, G. Carosi, N. Woollett, L. D. Duffy, R. Bradley, C. Boutan, M. Jones, B. H. Laroque, N. S. Oblath, M. S. Taubman, J. Clarke, A. Dove, A. Eddins, S. R. O’kelley, S. Nawaz, I. Siddiqi, N. Stevenson, A. Agrawal, A. V. Dixit, J. R. Gleason, S. Jois, P. Sikivie, J. A. Solomon, N. S. Sullivan, D. B. Tanner, E. Lentz, E. J. Daw, J. H. Buckley, P. M. Harrington, E. A. Henriksen, and K. W. Murch, Phys. Rev. Lett. 124, 101303 (2020).
  * Brubaker  _et al._ [2017a] B. M. Brubaker, L. Zhong, Y. V. Gurevich, S. B. Cahn, S. K. Lamoreaux, M. Simanovskaia, J. R. Root, S. M. Lewis, S. Al Kenany, K. M. Backes, I. Urdinaran, N. M. Rapidis, T. M. Shokair, K. A. van Bibber, D. A. Palken, M. Malnou, W. F. Kindel, M. A. Anil, K. W. Lehnert, and G. Carosi, Phys. Rev. Lett. 118, 061302 (2017a).
  * Crescini  _et al._ [2020] N. Crescini, D. Alesini, C. Braggio, G. Carugno, D. D’Agostino, D. Di Gioacchino, P. Falferi, U. Gambardella, C. Gatti, G. Iannone, C. Ligi, A. Lombardi, A. Ortolan, R. Pengo, G. Ruoso, and L. Taffarello, Phys. Rev. Lett. 124, 171801 (2020).
  * Callen and Welton [1951] H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).
  * Caves [1982] C. M. Caves, Phys. Rev. D 26, 1817 (1982).
  * Roy and Devoret [2016] A. Roy and M. Devoret, Comptes Rendus Physique 17, 740 (2016), quantum microwaves / Micro-ondes quantiques.
  * Yamamoto  _et al._ [2008] T. Yamamoto, K. Inomata, M. Watanabe, K. Matsuba, T. Miyazaki, W. D. Oliver, Y. Nakamura, and J. S. Tsai, Appl. Phys. Lett. 93, 042510 (2008).
  * Dolan [1977] G. J. Dolan, Appl. Phys. Lett. 31, 337 (1977).
  * Zhong  _et al._ [2013] L. Zhong, E. P. Menzel, R. Di Candia, P. Eder, M. Ihmig, A. Baust, M. Haeberlein, E. Hoffmann, K. Inomata, T. Yamamoto, Y. Nakamura, E. Solano, F. Deppe, A. Marx, and R. Gross, New Journal of Physics 15 (2013), 10.1088/1367-2630/15/12/125013.
  * Krantz  _et al._ [2013] P. Krantz, Y. Reshitnyk, W. Wustmann, J. Bylander, S. Gustavsson, W. D. Oliver, T. Duty, V. Shumeiko, and P. Delsing, New Journal of Physics 15, 105002 (2013).
  * Engen [1970] G. F. Engen, IEEE Transactions on Instrumentation and Measurement 19, 344 (1970).
  * Kim [1979] J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman, Vainshtein, and Zakharov [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nuclear Physics, Section B 166, 493 (1980).
  * Lee  _et al._ [2020] S. Lee, S. Ahn, J. Choi, B. R. Ko, and Y. K. Semertzidis, Phys. Rev. Lett. 124, 101802 (2020).
  * Jeong  _et al._ [2020] J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J. E. Kim, and Y. K. Semertzidis, Phys. Rev. Lett. 125, 221302 (2020).
  * Asztalos  _et al._ [2001] S. Asztalos, E. Daw, H. Peng, L. J. Rosenberg, C. Hagmann, D. Kinion, W. Stoeffl, K. van Bibber, P. Sikivie, N. S. Sullivan, D. B. Tanner, F. Nezrick, M. S. Turner, D. M. Moltz, J. Powell, M.-O. André, J. Clarke, M. Mück, and R. F. Bradley, Phys. Rev. D 64, 092003 (2001).
  * Brubaker  _et al._ [2017b] B. M. Brubaker, L. Zhong, S. K. Lamoreaux, K. W. Lehnert, and K. A. van Bibber, Phys. Rev. D 96, 123008 (2017b).
  * Yurke  _et al._ [1989] B. Yurke, L. R. Corruccini, P. G. Kaminsky, L. W. Rupp, A. D. Smith, A. H. Silver, R. W. Simon, and E. A. Whittaker, Phys. Rev. A 39, 2519 (1989).


# Neural Differential Algebraic Equations

James Koch*, Madelyn Shapiro, Himanshu Sharma, Draguna Vrabie, Ján Drgoňa
*Contact<EMAIL_ADDRESS>Pacific Northwest National Laboratory, Richland, WA

###### Abstract

Differential-Algebraic Equations (DAEs) describe the temporal evolution of
systems that obey both differential and algebraic constraints. Of particular
interest are systems that contain implicit relationships between their
components, such as conservation relationships. Here, we present Neural
Differential-Algebraic Equations (NDAEs) suitable for data-driven modeling of
DAEs. This methodology is built upon the concept of the Universal Differential
Equation; that is, a model constructed as a system of Neural Ordinary
Differential Equations informed by theory from particular science domains. In
this work, we show that the proposed NDAEs abstraction is suitable for
relevant system-theoretic data-driven modeling tasks. Presented examples
include (i) the inverse problem of tank-manifold dynamics and (ii) discrepancy
modeling of a network of pumps, tanks, and pipes. Our experiments demonstrate
the proposed method’s robustness to noise and extrapolation ability to (i)
learn the behaviors of the system components and their interaction physics and
(ii) disambiguate between data trends and mechanistic relationships contained
in the system.

##  I INTRODUCTION

Modeling and simulation of systems via Differential-Algebraic Equations (DAEs)
can be difficult, both in terms of formulation as well as numerical
implementation. Common issues include (i) over- or under-constraining states
and (ii) incorrect model specification, either of which can lead to a
prohibitively challenging modeling task [1]. While many such algebraic
relationships may be known (e.g., mass is conserved through a pipe), the
relationships between system states and these algebraic constraints are often
difficult to reconcile in a modeling framework (see, for example, index of a
DAE and choosing appropriate initial conditions [1]). Incorporating these
algebraic relationships into well-established data-driven modeling paradigms
is also non-trivial. Here, measurement noise, partial observations, and
missing or incorrect physics all confound the ability to effectively build
models that obey algebraic constraints.

Furthermore, it is generally rare to construct a DAE-based model of an
engineering system from first principles. More common is to adopt simplified
models that can capture improperly specified or missing physics at a given
fidelity. A common approach is to employ a closure model [2], i.e., a
surrogate for the full-fidelity physics inserted into a lower-fidelity model.
Closure models have enabled major advancements in scientific computing but
still traditionally rely upon rigorous reduction from first principles physics
to derive required simplified relationships [3].

To address this issue, we present a framework for the data-driven modeling of
DAEs with a novel neural timestepper scheme. Our framework leverages themes
from operator splitting methods [4] and Physics-Informed Neural Network
(PINN)-style [5] optimization problems to accommodate algebraic relationships.
We demonstrate the utility of these methods on two closure modeling tasks: (i)
a tank-manifold problem, where one needs to ‘invert’ the model to provide an
estimate of a property of a model component, and (ii) a tank-pipe-pump network
for learning unknown nonlinear interactions. Our specific contributions are:

  * •

Taking inspiration from the well-established field of operator splitting, we
extend neural timesteppers to appropriately integrate DAEs.

  * •

We accomplish this via an algebra surrogate layer that modifies algebraic
states of a system, while an ODE Solver integrates differential states.

  * •

We demonstrate the utility of the method in applications to canonical closure
and inverse modeling of constrained networked system dynamics.

The paper is organized as follows. In Section II, we give a brief overview of
recent methods proposed for data-driven modeling of agent-based systems.
Section III details our proposed methodology and its implementation. Example
results are provided in Section IV, followed by a discussion in Section V.

##  II BACKGROUND

###  II-A Differential-Algebraic Equations

Differential-Algebraic Equations (DAEs) are systems comprising differential
and algebraic relationships that describe the evolution of a system’s states.
Specifically, we are concerned with semi-explicit DAEs written as:

$\begin{split}\frac{dx}{dt}&=f\left(x,y,u\right)\\\
0&=g\left(x,y,u\right),\end{split}$ (1)

where $x\in\mathbb{R}^{n_{x}}$ are the differential states of a system that
evolve according to the vector field defined by
$f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{y}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}^{n_{x}}$,
$y\in\mathbb{R}^{n_{y}}$ are the algebraic states, and
$u\in\mathbb{R}^{n_{u}}$ are exogenous inputs. The algebraic states $y$ evolve
such that the algebraic relationships, defined by $g(x,y,u)=0$, are satisfied
for all time $t$.

A standard DAE solution approach is to transform the system into a set of ODEs
to be integrated with standard numerical solvers. This process is termed index
reduction [1]. The index of a DAE is a distance measure between a DAE and its
associated ODE: an index-1 DAE requires one differentiation of algebraic
relationships to yield a consistent set of ODEs, index-2 requires two
differentiations, etc. After differentiation and substitution into other
system relationships, the governing equations can be written as a set of
first-order ODEs.

###  II-B Neural Differential Equations

A Neural Ordinary Differential Equation (NODE) approximates the flow map of a
dynamical system with a neural network architecture. Consider a dynamical
system with states $x\in\mathbb{R}^{n_{x}}$. A NODE approximates the temporal
dynamics of $x$ according to the ODE:

$\frac{dx}{dt}=f\left(x;\theta\right),$ (2)

where the parametric flow $f:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{x}}$
maps the states to their time derivatives with tunable parameters $\theta$.
The solution to Eq. (2) is the initial value problem:

$x^{(t+\Delta t)}=x^{(t)}+\int_{t}^{t+\Delta t}f\left(x;\theta\right)dt,$ (3)

where the superscript contained in parentheses denotes indexing by time, not
exponentiation. Practically, this integration is performed with a standard
numerical ODE solver, which we denote as:

$x^{(t+\Delta t)}=\operatorname{ODESolve}\left(f,x^{(t)};\theta\right),$ (4)

where the limits of temporal integration are understood to be $[t,t+\Delta
t]$. In general, the NODE concept extends $f$ to any differentiable mapping by
implementing Eq. (4) in a differentiable manner [6]. This is achieved through
either (i) backpropagating residuals through the elementary functions of the
explicit time-stepping integrator (reverse-mode automatic differentiation) [7,
8], (ii) an adjoint sensitivity definition and solution [9], or (iii) forward-
mode automatic differentiation [10].

A fundamental assumption of a NODE model is that a direct relationship exists
between each system state and its temporal derivative. This is a reasonable
assumption for many cases, but poses challenges for partially-observed systems
[11] and application to DAEs. Index reduction is a reliable technique for
recasting a DAE into a set of ODEs for standard numerical integration [12],
but this is incompatible with the present data science problem. When governing
equations are unknown or incomplete, the equivalent set of ODEs and required
state augmentations are also unknown.

###  II-C Fractional Step Methods

Our methodology is inspired by fractional step methods for evolving complex
multiphysics problems. In fractional step methods (also termed ‘Operator
Splitting’), the evolution of a system is split into simpler, independent sub-
tasks. The underlying ODE or PDE is then solved explicitly as the sequential
solution of these sub-tasks, albeit at a decreased level of accuracy (i.e.,
splitting error [13]). There are several well-established mathematical
techniques for decomposing the solution of an ODE into sub-tasks. The simplest
scheme, Lie-Trotter operator splitting [4], assumes a differential operator
can be decomposed into two sub-tasks. For the ODE of Eq. 2, assume $f$ can be
decomposed into operators $A$ and $B$ such that:

$\frac{dx}{dt}=A(x)+B(x)$ (5)

for the interval $t\in[t,t+\Delta t]$. The Lie-Trotter splitting scheme
defines the sub-tasks as:

$\displaystyle\hat{x}^{(t+\Delta t)}$
$\displaystyle=\operatorname{ODESolve}\left(A(x),x^{(t)}\right)$ (6)
$\displaystyle x^{(t+\Delta t)}$
$\displaystyle=\operatorname{ODESolve}\left(B(x),\hat{x}^{(t+\Delta
t)}\right)$ (7)

The first sub-task solves the ODE defined by the operator $A(x)$ on the time
interval $[t,t+\Delta t]$. The second sub-task solves the ODE defined by the
operator $B(x)$ on the same time interval, but with the initial conditions
defined by the solution of the first sub-task. Many variations of this scheme
exist for different applications and orders of accuracy. For example, Strang
splitting [14] is a second-order accurate scheme that uses half-steps (i.e.,
steps of $\frac{\Delta t}{2}$) in the above scheme.

###  II-D Related Work

####  II-D1 Machine Learning for Dynamical Systems

Methods for learning dynamics from data exist on a continuum from ‘black-box’
models to ’white-box’ models [15]. The ‘black-box’ modeling task is performed
without regard to the underlying system or domain knowledge and includes
Recurrent Neural Networks (RNNs) and their variants (e.g., LSTMs) [16], neural
state space models [17], neural timesteppsers [8], sparse symbolic regression
[18], and Neural ODEs [9]. Often, such models extract salient features in the
training data to infer future states. However, they are at risk of poor
generalization (overfitting) and lack of interpretability, except for symbolic
regression techniques [18]. ‘White-box’ models use domain knowledge to
constrain models to match expected dynamics by construction and often lead to
parameter-tuning problems [19]. In ‘gray-box’ models, recently referred to as
physics-informed machine learning [20], prior knowledge is incorporated in a
manner to constrain a ‘black-box model’. Neural ODEs constructed with domain
knowledge, termed Universal Differential Equations [6, 21], are an attractive
modeling platform for these problems because of their flexibility in encoding
physics and other structural priors.

####  II-D2 Data-driven Modeling of DAEs

If the DAE is well-characterized, the parameter estimation, inverse design,
and control tasks are possible through an adjoint sensitivity analysis [22,
23, 24]. This strategy is generally successful if (i) the index of the
underlying DAE is low and (ii) all constituent physics are known.

Recent works have begun to address some common issues regarding data-driven
DAE-based models. Moya and Lin [25] use explicit Runge-Kutta time-stepping to
combine semi-explicit DAEs into machine learning frameworks, which they
demonstrate by modeling power networks. Xiao et al. [26] construct neural ODE
and DAE modules to forecast a regional power network. Huang et al. [27] adopt
a simplified approach to the construction of a surrogate system for sequence-
to-sequence mappings for DAEs. These studies share features of expressing the
learning task in terms of a semi-explicit DAE with varying degrees of
parameterization and integration of domain knowledge. In addition to these
studies, Zhong et al. [28] addressed contact modeling in multi-body dynamics,
the authors in the work handled a constraint update architecture with similar
explicit time-stepping as in Moya [25] and Xiao [26].

While the specific treatment of each of these domain-focused problems varies,
each leverages a form of a multi-step integration technique that separates the
evolution of differential and algebraic states. Here, we recognize that the
commonality of these works holds the potential for a more powerful
generalization of standard neural integration. Specifically, we recast these
ideas under the umbrella of fractional-step methods, where integration is
split into simpler, sequential sub-tasks. In effect, we propose a new class of
neural integrators that are well-suited for a wide variety of modeling tasks,
including DAEs, boundary value problems, and other multi-physics systems.

##  III METHODS

Consider the DAE system of Eq. 1. Similar to the approach of a Lie-Trotter
fractional step scheme, we split the evolution of a system’s differential and
algebraic states into sequential sub-tasks. The first sub-task updates the
algebraic states using a trainable deep neural network (DNN). The second sub-
task takes as input the updated algebraic states and performs a single-step
integration. These two sequential sub-tasks define one complete forward pass;
all system states are updated exactly once. More formally, consider a
differentiable ODE layer, denoted as $\operatorname{ODESolve}$. Using a
fractional step method as our solution technique for Eq. 1, we formulate a
Neural DAE (NDAE) as:

$\begin{split}y^{(t+\Delta
t)}&=h\left(x^{(t)},y^{(t)},u^{(t)};\theta_{h}\right)\\\ x^{(t+\Delta
t)}&=\operatorname{ODESolve}\left(f,\\{x^{(t)},y^{(t+\Delta
t)},u^{(t)}\\};\theta_{f}\right),\end{split}$ (8)

where
$h:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{y}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{y}}$
pairs with $g$ as an explicit algebra solver or trainable neural surrogate
thereof and the temporal integration bounds are assumed to be $[t,t+\Delta
t]$. We denote this complete set of operations as $\operatorname{DAESolve}$,
with one complete time step captured as:

$\begin{split}\\{x^{(t+\Delta t)},y^{(t+\Delta t)}\\}&=\\\
\operatorname{DAESolve}&\left(f,h,\\{x^{(t)},y^{(t)},u^{(t)}\\};\Theta\right),\end{split}$
(9)

where $\Theta=\\{\theta_{h},\theta_{f}\\}$ is the set of all model parameters.

Note that in this formulation, the updates to the algebraic and differential
states are sequential; that is, the $\operatorname{ODESolve}$ step is
dependent upon the update to the algebraic states. The sub-tasks do not
commute. This is an essential assumption that enables explicit updates for
both algebraic states and differential states. Figure 1 depicts the handling
of differential and algebraic states through the Neural DAE. The module
generates an internal state $y^{(t+\Delta t)}$ that is processed by the
subsequent $\operatorname{ODESolve}$ as an exogenous input, similar to the
notion of Neural Controlled Differential Equations [29].

\begin{overpic}[width=433.62pt]{fig_01.png} \end{overpic} Figure 1: Schematic
of the proposed DAE Timestepper.

Training the algebraic update sub-task is accomplished by formulating a PINN-
style penalty loss [5] in which known constraints are used to ‘penalize’ the
discrepancy between constraint satisfaction and predicted states. These
penalties are assessed in a multi-objective optimization loss function:

$\mathcal{L}(\Theta)=\lambda_{1}\mathcal{L}_{\text{residual}}(X,\hat{X};\Theta)+\lambda_{2}\mathcal{L}_{\text{constraints}}(X,\hat{X};\Theta),$
(10)

where $\mathcal{L}_{\text{residual}}$ is a measure of the residual between the
training dataset and the model’s output, $\mathcal{L}_{\text{constraints}}$ is
a regularization term penalizing the constraint violations. The $\lambda$’s
represent weightings for balancing the multi-objective loss function. We
optimize model parameters $\Theta$ with an off-the-shelf stochastic gradient-
based optimizer [30] to solve
$\Theta^{*}=\underset{\Theta}{\text{argmin}}~{}\mathcal{L}(\Theta)$.

For this work, we choose minimal architectures for ease of implementation and
exposition, as outlined in Alg. 1. We use a first-order explicit timestepper
for the ODE integration and use a DNN as a neural surrogate of an algebra
solver.

Data:
$\hat{X}=\\{(\hat{x}^{(0)},\hat{y}^{(0)}),...,(\hat{x}^{(N)},\hat{y}^{(N)})\\}$:
Time series dataset with spacing $\Delta t$

Parameters : $\theta_{f},\theta_{h}$: Tunable parameters for mappings $f$ and
$h$

/* Pick appropriate integrator */

1 Function  _ODESolve(_$f,\\{x,y,u\\}$_)_:

2 $x^{*}\longleftarrow x+\Delta t\times f\left(x,y,u;\theta_{f}\right)$

3 return $x^{*}$

/* Sequential sub-tasks */

4 Function  _DAESolve(_$f,h,\\{x,y,u\\}$_)_:

5 $y^{*}\longleftarrow h\left(x,y,u;\theta_{h}\right)$

6
$x^{*}\longleftarrow\textnormal{{ODESolve}}\left(f,\\{x,y^{*},u\\};\theta_{f}\right)$

7 return $(x^{*},y^{*})$

/* Training loop */

8 while  _not converged_ do

9

/* Forward pass */

10 $X\longleftarrow\\{(\hat{x}^{(0)},\hat{y}^{(0)})\\}$

11 for  _k $\in 1,\dots,N$_ do

12 $X\longleftarrow$Concat ($X$,

13 DAESolve(_$f,h,\\{x^{(k-1)},y^{(k-1)},u^{(k-1)}\\}$_) )

14

15

/* Evaluate loss */

16 $\mathcal{L}\longleftarrow$ Criterion ( $X,\hat{X}$ )

17

18 if  _converged_ then

19 break;

20

21 else

22 obtain gradients;

23 update model weights $\\{\theta_{f},\theta_{h}\\}$;

24

25

Algorithm 1 Neural DAE training

##  IV NUMERICAL CASE STUDIES

All experiments were performed using the NeuroMANCER Scientific Machine
Learning (SciML) package [31] built on top of PyTorch [32]111Repository
containing all training and visualization scripts will be open-sourced and
publicly available by the time of publication.. Each presented model was
trained using an Intel I9-9980 CPU laptop with a wall-clock training time on
the order of 5 minutes.

###  IV-A Tank-Manifold Property Inference

####  IV-A1 Problem Description

A manifold is a physical realization of a conservation relationship that
manifests as an algebraic constraint in a DAE. Figure 2 depicts a pair of
tanks fed by a manifold. A pseudo-infinite reservoir supplies the pump. At the
manifold, the inflow is exactly equal to the sum of the outflows; that is,
$y_{\text{in}}=y_{1}+y_{2}$, where $y$ represents volumetric flow rates and
the subscripts $1$ and $2$ denote the two exit streams of the manifold. In
this example, the two tanks have different area-height profiles but have a
common datum such that the column heights are identical. As a fluid is pumped
into the two tanks, the flow rates in the two streams varies to ensure (i) the
conservation relationship of the manifold and (ii) that each of the tanks have
identical pressure head, i.e $x_{1}=x_{2}$, where $x$ is the height of the
fluid in the respective tank. The full DAE that describes the evolution of
this system is given as:

$\begin{split}\frac{dx_{1}}{dt}&=\frac{y_{1}}{\phi_{1}(x_{1})}\\\
\frac{dx_{2}}{dt}&=\frac{y_{2}}{\phi_{2}(x_{2})}\end{split}$ (11)
$\displaystyle 0=y_{\text{in}}-y_{1}-y_{2}\quad$ $\displaystyle
0=x_{1}-x_{2},$

where $\phi_{1}$ and $\phi_{2}$ denote the area-height profiles for the two
tanks. The solution of Eq. 11 can be obtained with many off-the-shelf standard
scientific computing packages. For this study, the area-height profiles are
set to $\phi_{1}(x)=3.0$ and $\phi_{2}(x)=\sqrt{x}+0.1$ respectively. The goal
of the learning task is to recover $\phi_{2}(x)$ given observations of tank
heights $x$ and pump volumetric flow rates $y$. Note that if the task were to
recover both $\phi_{1}(x)$ and $\phi_{2}(x)$ simultaneously, the problem would
be ill-posed in the sense that there does not exist a unique set of tank
profiles that recovers the behavior seen in observations.

\begin{overpic}[width=303.53267pt]{fig_02.png} \end{overpic} Figure 2: Tank-
manifold-pump schematic. A reservoir supplies fluid to a pump-manifold system.
The tanks have different area-height profiles but a common datum. Conservation
is enforced at the manifold.

####  IV-A2 Model Construction

In the system described by Eq. 11, the differential states are column heights
$x_{1}$ and $x_{2}$ and the algebraic states are the volumetric flows $y_{1}$
and $y_{2}$. The input flow is treated as an exogenous control input, i.e.,
$u=y_{\text{in}}$.

We aim to approximate the area-height profile of Tank #1 as well as construct
a surrogate for the behavior of the manifold such that (i) fluid is conserved
and (ii) the fluid column height constraint is satisfied. There are two neural
networks that we seek to train:
$NN_{1}:\mathbb{R}^{1}\rightarrow\mathbb{R}^{1}$, an approximation of the
sought area-height profile $\phi_{1}$, and
$NN_{2}:\mathbb{R}^{4}\rightarrow(0,1)$, the surrogate of the behavior of the
manifold. Each are implemented as a feed-forward neural network with one
hidden layer of 5 nodes, each with sigmoid activation functions. The output
layer of $NN_{1}$ is linear while the output layer of $NN_{2}$ is sigmoidal,
constraining the output to $(0,1)$.

####  IV-A3 Optimization Problem

A time series dataset is constructed by simulation of the system given in Eq.
11 regularly sampled on the interval $t\in[0,500]$ with $\Delta t=1.0$ and
with a constant inflow of $y_{\text{in}}=u=0.5$:

$\hat{X}=\\{(\hat{x}^{(0)},\hat{y}^{(0)}),...,(\hat{x}^{(N)},\hat{y}^{(N)})\\},$
(12)

where $N$ is the time series prediction horizon and the $\widehat{\cdot}$
denotes membership to the dataset. The optimization problem seeks the minimum
of the differential state trajectory reconstruction loss and algebraic state
reconstruction loss when predicted by our $\operatorname{DAESolve}$:

$\begin{split}\underset{\theta_{1},\theta_{2}}{\text{min}}&~{}~{}\lambda_{1}\sum_{k=1}^{N}||x^{(k)}-\hat{x}^{(k)}||_{2}^{2}+\lambda_{2}\sum_{k=1}^{N}||y^{(k)}-\hat{y}^{(k)}||_{2}^{2}\\\
\text{s.t.~{}~{}}&\\{x^{(k+1)},y^{(k+1)}\\}\\\
~{}&=\operatorname{DAESolve}\left(f,h,\\{x^{(k)},y^{(k)},u^{(k)}\\};\\{\theta_{1},\theta_{2}\\}\right),\\\
\end{split}$ (13)

The functions $f$ and $h$ are defined as:

$\begin{split}f&:\begin{cases}\frac{y_{1}}{3}\\\
\frac{y_{2}}{\operatorname{NN}_{1}(x_{2};\theta_{1})},\\\ \end{cases}\\\
h&:\begin{cases}u\operatorname{NN}_{2}(x_{1},x_{2},y_{1},y_{2};\theta_{2})\\\
u\left(1-\operatorname{NN}_{2}(x_{1},x_{2},y_{1},y_{2};\theta_{2})\right),\\\
\end{cases}\\\ \end{split}$ (14)

Note that the mapping $h$ defines a convex combination at the manifold: the
outflows $y_{1}$ and $y_{2}$ will sum to $u$ by construction.

For all examples contained in this work, the optimizer is Adam [30] with fixed
learning rate of 0.001. Training is deemed ‘converged’ after 20,000 epochs or
until an early stop is triggered after 20 epochs of increasing loss.

####  IV-A4 Results

The time series data for the tank heights and volumetric flows are given in
Figs. 3(a) and 3(b) respectively. The trained model is able to accurately
reproduce the time series data seen in training. Additionally, it is able to
accurately learn the model of the unknown area-height profile of Tank #2.
Figure 4(a) depicts the ground-truth area-height profile and the model’s
reconstruction. Furthermore, to demonstrate the generalization, the trained
model is tested on an unseen inlet mass flow rate;
$u(t)=\frac{1}{2}+\frac{1}{4}\sin{\left(\frac{t}{100}\right)}$. The resulting
time series of volumetric flows from the trained model rollout is compared
with the ground-truth conditions in Fig. 4(b). Table-I gives the Mean Squared
Error (MSE) of the trained neural DAE model.

Figure 3: (a) Tank heights. (b) Volumetric flow. Figure 4: Trained model (a)
Area-height relationships (b) Volumetric flow time series.

Extrapolation to a new flow condition demonstrates the ability to (i)
endogenize the behaviors of the components in the network and their
interaction physics and (ii) disambiguate between data trends and mechanistic
relationships contained in the network.

TABLE I: Mean Square error for the differential state and algebraic states for tanks-manifold-pump experiment. Manifold |  | Mean Square Error
---
(MSE)
Tanks Height | 9e-03
Volumetric Flow | 2e-02
Area-Height | 6e-03
Unseen inlet mass-flow | 1e-01

###  IV-B Tank Network Modeling

####  IV-B1 Problem Description

Figure 5 expands upon the tank-manifold problem of Fig. 2 by adding additional
tanks, feedback mechanisms, and ‘closing the loop’ in the network flow. In
this example, the pump flow rate depends on the fluid heights of Tank #1 and
the reservoir. Tanks #1 and #2 have a height constraint (equivalent pressure
head) and are fed by a common manifold. The outlets of tanks #1 and #3 depend
on the square root of the fluid column height in the respective tanks. The
full system can be described by the DAE:

$\begin{split}\frac{dx_{1}}{dt}&=\frac{1}{\phi_{1}(x_{1})}\left(y_{1}-y_{3}\right)\\\
\frac{dx_{2}}{dt}&=\frac{1}{\phi_{2}(x_{2})}\left(y_{2}\right)\\\
\frac{dx_{3}}{dt}&=\frac{1}{\phi_{3}(x_{3})}\left(y_{3}-y_{4}\right)\\\
\frac{dx_{4}}{dt}&=\frac{1}{\phi_{4}(x_{4})}\left(y_{4}-y_{0}\right)\\\
\end{split}$ (15) $\displaystyle 0=y_{0}-p(x_{1},x_{4})\quad$ $\displaystyle
0=y_{0}-y_{1}-y_{2}$ $\displaystyle 0=y_{3}-\alpha_{1}\sqrt{x_{1}}\quad$
$\displaystyle 0=y_{4}-\alpha_{2}\sqrt{x_{3}}$ $\displaystyle
0=x_{1}-x_{2}\quad$

where $x$ denotes the fluid column height for tanks #1 through #4 (the
reservoir), $p$ is a pump response function, $y$ represents volumetric flows,
and $\alpha$ describe the discharge coefficients of the tank outlets. The
differential states are $x=\\{x_{1},x_{2},x_{3},x_{4}\\}\in\mathbb{R}^{4}$ and
the algebraic states are
$y=\\{y_{0},y_{1},y_{2},y_{3},y_{4}\\}\in\mathbb{R}^{5}$.

In this example, several characteristics make the dynamics more complex.
First, the system is closed, as there is no material transfer into or out of
the system. Second, there is the potential for backflow through the manifold
from Tank #2. Lastly, the feedback mechanism from the tank level transducers
to the pump controller is such that oscillations can occur.

\begin{overpic}[width=303.53267pt]{fig_05.png} \end{overpic} Figure 5:
Networked system defined by tanks, pipes, a manifold, and a pump. The system
exhibits oscillatory dynamics for some parameter choices and sensing schemes.

####  IV-B2 Optimization Problem

The system given in Eq. (15) was simulated with the following specifications:

$\displaystyle\phi_{1}(x_{1})=2.0\quad$ $\displaystyle\phi_{2}(x_{2})=1.0$
$\displaystyle\phi_{3}(x_{3})=1.0\quad$ $\displaystyle\phi_{4}(x_{4})=10.0$
$\displaystyle p(x_{1},x_{4})=0.1x_{1}x_{4}\quad$
$\displaystyle\alpha_{1},\alpha_{2}=0.1$ (16)

A single-trajectory dataset $\hat{X}$ was constructed through sampling this
simulation at $\Delta t=0.1$ for $t\in[0,20]$. The optimization problem
definition is identical to that of Eq. 13. For this problem, it is reasonable
to assume that the basic physics for a gravity-fed pipe network are known,
which provides a starting point for constructing a model. However, discharge
coefficients and quantifying feedback mechanisms are not as straightforward.
Here, we attempt to approximate (i) the Level Controller (LC) dynamics for
pump control with a neural network
$\operatorname{NN}_{1}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$, (ii) manifold
dynamics with with a second neural network
$\operatorname{NN}_{2}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$, and (iii) the
discharge coefficients $\alpha_{1}$ and $\alpha_{2}$. The mappings $f$ and $h$
are defined as:

$\begin{split}f&:\begin{cases}\frac{y_{1}-y_{3}}{\phi_{1}(x_{1})}\\\
\frac{y_{2}}{\phi_{2}(x_{2})}\\\ \frac{y_{3}-y_{4}}{\phi_{3}(x_{3})}\\\
\frac{y_{4}-y_{0}}{\phi_{4}(x_{4})},\\\ \end{cases}\\\
h&:\begin{cases}\operatorname{NN}_{1}(x_{1},x_{4};\theta_{1})\\\
y_{0}\operatorname{NN}_{2}(x_{2},x_{3};\theta_{2})\\\
y_{0}(1-\operatorname{NN}_{2}(x_{2},x_{3};\theta_{2}))\\\
\alpha_{1}\sqrt{x_{1}}\\\ \alpha_{2}\sqrt{x_{3}}.\\\ \end{cases}\end{split}$
(17)

For this example, both $\operatorname{NN}_{1}$ and $\operatorname{NN}_{2}$ are
simple multi-layer perceptrons each with a two hidden layers of size 30 and
ReLU activations.

####  IV-B3 Results

Figure 6(a) details the time history of the four tank heights for both the
ground truth data as well as the open-loop model rollout. At these process
conditions, the pump delivers an oscillatory flow to the manifold, which our
model is able to capture. Additionally, the tuned model is able to infer the
values for the discharge coefficients: for the case of zero noise the
recovered values are $\alpha_{1}=0.1027$ and $\alpha_{2}=0.1024$.

To demonstrate generalizability, the trained model is deployed on an unseen
area-height tank profile; for this, we redefine $\phi_{1}(x_{1})=1.0$. In
doing so, the system’s limit cycle is destroyed, and the true dynamics
approach a steady state. Figure 6(b) demonstrates this extrapolation for the
same time horizon. The reconstructed time series data matches the behavior of
the unseen data: there exist damped oscillations that decay to a steady state
value. Table-II shows the MSE for the neural DAE models trained for the tank
network problem.

Noise tolerance is an important consideration in the evaluation of system
identification tasks, especially in settings with nonlinearities and feedback
mechanisms. For this example, we investigate noise tolerance by introducing
Gaussian white noise. Figures 7(a) and 7(b) show similar training and
extrapolation tasks for the Signal-to-Noise Ratio (SNR) of $20$ dB. Despite
the significant noise, the learned dynamics are faithful to the underlying
physics, as demonstrated in the extrapolation task, and result in an MSE of
$7e-02$.

Figure 6: Tank Network problem : (a) Ground truth data vs Neural-DAE model
prediction. (b) Extrapolation case: inferring behavior for unseen component
parameters using Neural-DAE. Figure 7: Tank Network problem : (a) noisy data
vs Neural-DAE model prediction. (b) Extrapolation case: inferring behavior for
unseen component parameters using Neural-DAE. TABLE II: Mean Square error for
the differential state and algebraic states for tanks network experiment.

Trained Model Type | Tank-Network Experiment |  | Mean Square Error
---
(MSE)
No Noise | Tanks Height | 1e-02
Unseen initial Data : Tanks Height | 6e-02
20 dB Noise | Tanks Height | 1.03
Unseen initial Data : Tanks Height | 7e-2

##  V CONCLUSIONS AND FUTURE WORKS

This work presents a framework for the data-driven modeling of Differential
Algebraic Equations (DAEs). Specifically, we introduce a novel method called
Neural Differential Algebraic Equations (NDAEs) that extends the notion of
Universal Differential Equations to endogenize the evolution of algebraic
relationships. Inspired by Lie-Trotter operator splitting, we achieve this by
splitting the integration of system states into two sequential sub-tasks:
first, updating algebraic states via a differentiable surrogate of an algebra
solver, followed by a conventional temporal integration step for a system of
Ordinary Differential Equations.

The effectiveness of these methods was demonstrated by (i) modeling networked
dynamical systems with underlying constraints and conservation relationships
and (ii) showing generalization to unseen process conditions or component
variation. We demonstrate that the proposed NDAE method is robust to noise and
possesses the extrapolation ability to (i) learn the behaviors of the DAE
system components and their interaction physics and (ii) disambiguate between
data trends and mechanistic relationships contained in the DAE system.

While the methods and example application problems were successful in this
work, future research and development are needed in several areas. First, a
central assumption of our method is that the underlying dynamics of the
systems of interest are sufficiently smooth, i.e., not containing jump
discontinuities or shocks, such that the underlying dynamics are non-stiff. As
presented, our method leverages an explicit fractional-step approach to
decoupling the evolution of the differential and algebraic states. To increase
the ability of the method to handle stiff dynamics or discontinuities, moving
to an implicit integration scheme may be warranted. Second, further analysis
is needed to investigate the interplay between the modeling task and the index
of the underlying DAE (for example, higher-order DAEs that are characteristic
of design or parameter optimization problems). Lastly, we have made no claims
about computational scalability nor partial observability, both of which are
critical for modeling real-world complex systems with hundreds of degrees of
freedom. We reserve these topics for future work.

## ACKNOWLEDGMENTS

This research was supported by the U.S. Department of Energy through the
Building Technologies Office under the Advancing Market-Ready Building Energy
Management by Cost-Effective Differentiable Predictive Control projects. PNNL
is a multi-program national laboratory operated for the U.S. Department of
Energy (DOE) by Battelle Memorial Institute under Contract No. DE-
AC05-76RL0-1830.

## References

  * [1] P. Kunkel, Differential-algebraic equations: analysis and numerical solution, vol. 2.  European Mathematical Society, 2006.
  * [2] G. L. Mellor and T. Yamada, “Development of a turbulence closure model for geophysical fluid problems,” Reviews of Geophysics, vol. 20, no. 4, pp. 851–875, 1982.
  * [3] R. H. Bush, T. S. Chyczewski, K. Duraisamy, B. Eisfeld, C. L. Rumsey, and B. R. Smith, “Recommendations for future efforts in rans modeling and simulation,” in AIAA scitech 2019 forum, p. 0317, 2019.
  * [4] H. Holden, Splitting methods for partial differential equations with rough solutions: Analysis and MATLAB programs, vol. 11.  European Mathematical Society, 2010.
  * [5] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nature Reviews Physics, vol. 3, no. 6, pp. 422–440, 2021.
  * [6] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal differential equations for scientific machine learning,” arXiv preprint arXiv:2001.04385, 2020.
  * [7] S. H. Rudy, J. N. Kutz, and S. L. Brunton, “Deep learning of dynamics and signal-noise decomposition with time-stepping constraints,” Journal of Computational Physics, vol. 396, pp. 483–506, 2019.
  * [8] Y. Liu, J. N. Kutz, and S. L. Brunton, “Hierarchical deep learning of multiscale differential equation time-steppers,” Philosophical Transactions of the Royal Society A, vol. 380, no. 2229, p. 20210200, 2022.
  * [9] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” Advances in neural information processing systems, vol. 31, 2018.
  * [10] J. Revels, M. Lubin, and T. Papamarkou, “Forward-mode automatic differentiation in julia,” arXiv preprint arXiv:1607.07892, 2016.
  * [11] E. Dupont, A. Doucet, and Y. W. Teh, “Augmented neural odes,” Advances in neural information processing systems, vol. 32, 2019.
  * [12] S. E. Mattsson, H. Elmqvist, and M. Otter, “Physical system modeling with modelica,” Control engineering practice, vol. 6, no. 4, pp. 501–510, 1998.
  * [13] J. B. Perot, “An analysis of the fractional step method,” Journal of Computational Physics, vol. 108, no. 1, pp. 51–58, 1993.
  * [14] G. Strang, “On the construction and comparison of difference schemes,” SIAM journal on numerical analysis, vol. 5, no. 3, pp. 506–517, 1968.
  * [15] C. Legaard, T. Schranz, G. Schweiger, J. Drgoňa, B. Falay, C. Gomes, A. Iosifidis, M. Abkar, and P. Larsen, “Constructing neural network based models for simulating dynamical systems,” ACM Comput. Surv., vol. 55, feb 2023.
  * [16] K.-i. Funahashi and Y. Nakamura, “Approximation of dynamical systems by continuous time recurrent neural networks,” Neural networks, vol. 6, no. 6, pp. 801–806, 1993.
  * [17] E. Skomski, S. Vasisht, C. Wight, A. Tuor, J. Drgoňa, and D. Vrabie, “Constrained block nonlinear neural dynamical models,” in 2021 American Control Conference (ACC), pp. 3993–4000, 2021.
  * [18] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the national academy of sciences, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [19] J. O. Ramsay, G. Hooker, D. Campbell, and J. Cao, “Parameter estimation for differential equations: a generalized smoothing approach,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 69, no. 5, pp. 741–796, 2007.
  * [20] T. X. Nghiem, J. Drgoňa, C. Jones, Z. Nagy, R. Schwan, B. Dey, A. Chakrabarty, S. Di Cairano, J. A. Paulson, A. Carron, M. N. Zeilinger, W. Shaw Cortez, and D. L. Vrabie, “Physics-informed machine learning for modeling and control of dynamical systems,” in American Control Conference (ACC), pp. 3735–3750, 2023.
  * [21] J. Koch, Z. Chen, A. Tuor, J. Drgoňa, and D. Vrabie, “Structural inference of networked dynamical systems with universal differential equations,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 33, no. 2, p. 023103, 2023.
  * [22] L. Petzold, R. Serban, S. Li, S. Raha, and Y. Cao, “Sensitivity analysis and design optimization of differential-algebraic equation systems,” in NATO Advanced Research Workshop on Computational Aspects of Nonlinear Structural Systems with Large Rigid Body Motion, Pultusk, Poland, 2000.
  * [23] Y. Cao, S. Li, L. Petzold, and R. Serban, “Adjoint sensitivity analysis for differential-algebraic equations: The adjoint dae system and its numerical solution,” SIAM journal on scientific computing, vol. 24, no. 3, pp. 1076–1089, 2003.
  * [24] J. B. Jorgensen, “Adjoint sensitivity results for predictive control, state- and parameter-estimation with nonlinear models,” in 2007 European Control Conference (ECC), pp. 3649–3656, 2007.
  * [25] C. Moya and G. Lin, “Dae-pinn: a physics-informed neural network model for simulating differential algebraic equations with application to power networks,” Neural Computing and Applications, vol. 35, no. 5, pp. 3789–3804, 2023.
  * [26] T. Xiao, Y. Chen, S. Huang, T. He, and H. Guan, “Feasibility study of neural ode and dae modules for power system dynamic component modeling,” IEEE Transactions on Power Systems, 2022.
  * [27] Y. Huang, C. Zou, Y. Li, and T. Wik, “Minn: Learning the dynamics of differential-algebraic equations and application to battery modeling,” arXiv preprint arXiv:2304.14422, 2023.
  * [28] Y. D. Zhong, B. Dey, and A. Chakraborty, “Extending lagrangian and hamiltonian neural networks with differentiable contact models,” Advances in Neural Information Processing Systems, vol. 34, pp. 21910–21922, 2021.
  * [29] P. Kidger, J. Morrill, J. Foster, and T. Lyons, “Neural controlled differential equations for irregular time series,” Advances in Neural Information Processing Systems, vol. 33, pp. 6696–6707, 2020.
  * [30] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [31] J. Drgoňa, A. Tuor, J. Koch, R. Birmiwal, M. Shapiro, and D. Vrabie, “NeuroMANCER: Neural Modules with Adaptive Nonlinear Constraints and Efficient Regularizations,” 2023.
  * [32] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., “Pytorch: An imperative style, high-performance deep learning library,” Advances in neural information processing systems, vol. 32, 2019.