\section{Introduction}\label{sec.intro}
\input{sections/1_introduction}

\section{Results}\label{sec.results}
\input{sections/2_1_model}
\input{sections/2_2_policies}

\section{Discussion}\label{sec.discussion}
\input{sections/3_discussion}

\section{Methods}\label{sec.methods}
\input{sections/4_methods}

\section{Data availability}\label{sec.data_availability}
The data shown in this paper can be found at \href{https://doi.org/10.4121/20402037.v1}{https://doi.org/10.4121/20402037.v1}.

\section{Code availability}\label{sec.code_availability}
Our code can be found in the following GitHub repository: \href{https://github.com/AlvaroGI/optimal-homogeneous-chain}{https://github.com/AlvaroGI/optimal-homogeneous-chain}.

\bibliographystyle{unsrt}

\subsection{Optimal entanglement distribution policies}
As described above, nodes always attempt entanglement generation if there are qubits available.
Cutoffs are always applied whenever an entangled state becomes too old.
Since nodes are free to attempt swaps as soon as entangled links are available or some time later, they must agree on an \emph{entanglement distribution policy}: a set of rules that indicate when to perform a swap.
The \emph{entanglement delivery time} of a policy is the time required to generate end-to-end entanglement when following said policy.
An \emph{optimal policy} minimizes the expected entanglement delivery time.
In a \emph{global-knowledge policy}, nodes have information about all the entangled links in the chain.
In a \emph{local-knowledge policy}, the nodes only need to know the state of the qubits they hold.
An example of local policy is the \emph{swap-asap policy}, in which each node performs a swap as soon as both entangled links are available.

We model the evolution of the state of the repeater chain as an MDP.
We then formulate the Bellman equations \cite{Sutton2018} and solve them using value iteration and policy iteration to find global-knowledge optimal policies.
More details and formal definitions are provided in Section \ref{sec.methods}.

Let us now describe the relation between the expected delivery time of an optimal policy, $T_\mathrm{opt}$, and the variables of the system ($n$, $p$, $p_\mathrm{s}$, and $t_\text{cut}$).
Repeater chains with a larger number of nodes $n$ yield a larger $T_\mathrm{opt}$, since more entangled links need to be generated probabilistically.
When $p$ is small, more entanglement generation attempts are required to succeed, yielding a larger $T_\mathrm{opt}$.
Decreasing $p_\text{s}$ also increases $T_\mathrm{opt}$, since more attempts at entanglement swapping are required on average.
When $\tc$ is small, all entangled states must be generated within a small time window and therefore $T_\mathrm{opt}$ is also larger.
Figure~\ref{fig.delivery-time-n5} shows the expected delivery time of an optimal policy in a five-node chain.
Interestingly, $p_\text{s}$ has a much stronger influence on $T_\mathrm{opt}$ than $p$ and $\tc$: decreasing $p_\text{s}$ from 1 to 0.5 in a five-node chain translates into an increase in $T_\mathrm{opt}$ of an order of magnitude.
Similar behavior is observed for other values of $n$, as shown in Appendix~\ref{app.optimal}.

\begin{figure}[t!]
\centering
\includegraphics[width=0.9\textwidth]{delivery-time_n5_tol1e-07.pdf}
\caption{\textbf{The expected delivery time increases with lower $p$, $p_\mathrm{s}$, and $t_\mathrm{cut}$.}
	Expected delivery time of an optimal policy, $T_\mathrm{opt}$, versus $p$ in a five-node chain, for different values of cutoff ($t_\text{cut}=2,5,10$). Solid lines correspond to deterministic swaps ($p_\text{s}=1$) and dashed lines correspond to probabilistic swaps with $p_\text{s}=0.5$.}
	\label{fig.delivery-time-n5}
\end{figure}

To evaluate the advantages of an optimal policy, we use the swap-asap policy as a baseline.
Early swaps can provide an advantage in terms of delivery time, since swapping earlier can free up qubits that can be used to generate backup entangled links, as displayed in the first transition in Figure~\ref{fig.chain_example}.
However, the age of a swapped link may reach the cutoff time earlier than one of the input links consumed in the swap, as the swapped link assumes the age of the oldest input link. Following the example in Figure~\ref{fig.chain_example} and assuming $\tc=1$, if no swaps are performed, the links between nodes two and three and between three and four will exist for one more time slot, while the link between nodes four and five will be removed immediately since it reached the cutoff time. If both swaps are performed, the swapped link between nodes two and five will be removed immediately since it reached the cutoff time.
Since we have arguments in favor of and against swapping early, it is not trivial to determine the scenarios in which the swap-asap policy is close to optimal. Next, we compare the expected delivery times of an optimal global-knowledge policy and the swap-asap policy.

\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{FiveLinkChainEx}
\caption{\textbf{Swap-asap policies free up qubits, but swapped links expire earlier.}
Evolution of an example state when following a waiting policy versus the swap-asap policy during a single time slot.
Entangled links are represented as solid black lines, with occupied qubits in black and free qubits in white.
A waiting policy decides to not perform any swap, while the swap-asap policy decides to swap all three links.
The swap frees up qubits (marked in orange) that can be used to resume entanglement generation either if the swap is successful, as in the picture, or not.
After performing swaps, a cutoff $\tc=1$ is applied and links with age 1 are removed, causing the swapped link to expire.}
	\label{fig.chain_example}
\end{figure}

Figure~\ref{fig.advantage5} shows the relative difference between the expected delivery times of an optimal global-knowledge policy, $T_\text{opt}$, and that of the swap-asap policy, $T_\text{swap}$, in a five-node chain.
Increasing values of $(T_\text{swap}-T_\text{opt})/T_\text{opt}$ mean that the optimal policy is increasingly faster on average.
Note that we restrict our analysis to the parameter regime $p\geq0.3$ and $2\leq\tc\leq6$ due to the very large computational cost of calculating the solution for smaller $p$ and larger $\tc$ (for more details, see Section \ref{sec.methods}).
Let us first focus on deterministic swaps (Figure~\ref{fig.advantage5-ps1}).
The advantage provided by an optimal policy increases for decreasing $p$.
When $p$ is small, links are more valuable since they are harder to generate. Therefore, it is convenient to avoid early swaps, as they effectively increase the ages of the links involved and make them expire earlier.
When $\tc$ is small, a similar effect happens: all entangled links must be generated within a small time window and early swaps can make them expire too soon. For larger $\tc$, increasing the age of a link does not have a strong impact on the delivery time, since the time window is larger. Therefore, an optimal policy is increasingly better than swap-asap for decreasing $\tc$.
The maximum difference between expected delivery times in the parameter region explored is 5.25\%.

Interestingly, probabilistic swaps (Figure~\ref{fig.advantage5-ps05}) yield an opposite behavior in the parameter region explored: optimal policies are increasingly better than swap-asap for increasing $p$ and $\tc$ (except when $p\leq0.4$ and $\tc\leq3$), and the relative difference in expected delivery time can be as large as 13.2\% (achieved in a five-node chain with $p=0.9$ and $\tc=6$).
One reason for this may be the action that each policy decides to perform when the repeater chain is in a \emph{full state}, which is a situation where each pair of neighboring nodes shares an entangled link (see state at the top of Figure~\ref{fig.swapasap-vs-doubling}).
When swaps are deterministic, the optimal policy chooses to swap all links in a full state, since end-to-end entanglement will always be achieved.
However, when swaps are probabilistic, an optimal policy generally chooses to perform two separate swaps (see Figure~\ref{fig.swapasap-vs-doubling}), similar to the nested purification scheme proposed in \cite{Briegel1998}.
As an example, for $n=5$, $p=0.9$, $\tc=2$, and $p_\text{s}=0.5$, the swap-asap policy yields an expected delivery time of $T=9.35$.
If, in full states, the swap at the third node is withheld, $T$ drops to 8.34.
The swap-asap policy is on average slower than this modified policy by 12.1\%.
The action chosen in full states has a stronger influence on $T$ for increasing $p$. This is because full states are more frequent for large $p$: whenever a swap fails, a full state is soon recovered, since new entangled states are generated with high probability.
As a consequence, an optimal policy is increasingly better than swap-asap for higher $p$ when swaps are probabilistic.
A similar effect happens for large $\tc$.
Note however that the effect of the action chosen in full states is practically irrelevant in four-node chains (see Appendix~\ref{app.optimal}).
Note also that the advantage of an optimal policy in terms of delivery time is not always monotonic in $p$ and $t_\text{cut}$ (see Appendix~\ref{app.optimal}).

\begin{figure}[t!]
\captionsetup[subfigure]{justification=centering}
     \centering
     \begin{subfigure}[b]{0.9\textwidth}
         \centering
         \includegraphics[width=\textwidth]{advantage__n5_ps1.000_tolerance1e-07_randomseed2.pdf}
	\vspace{-1.5\baselineskip}
         \caption{Deterministic swaps ($p_\mathrm{s}=1$).}
	\label{fig.advantage5-ps1}
         \vspace{5pt}
     \end{subfigure}

     \begin{subfigure}[b]{0.9\textwidth}
         \centering
         \includegraphics[width=\textwidth]{advantage__n5_ps0.500_tolerance1e-07_randomseed2.pdf}
	\vspace{-1.5\baselineskip}
         \caption{Probabilistic swaps ($p_\mathrm{s}=0.5$).}
	\label{fig.advantage5-ps05}
         \vspace{5pt}
     \end{subfigure}

	\caption{\textbf{In a five-node chain, an optimal policy performs increasingly better than swap-asap for lower/higher values of $p$ and $t_\mathrm{cut}$ when swaps are deterministic/probabilistic.}
	Relative difference between the expected delivery times of an optimal policy, $T_\mathrm{opt}$, and the swap-asap policy, $T_\mathrm{swap}$, in a five-node chain, for different values of $p$ and $\tc$.
	}
	\label{fig.advantage5}
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[width=\textwidth]{figs/swapasap-vs-doubling}
	\caption{\textbf{All possible transitions after performing a swap-asap action or a nested action in a full state, depending on which swaps succeed.}
	In full states, every pair of neighbors shares an entangled link (solid black lines, with occupied qubits in black and free qubits in white).
	The swap-asap policy decides to swap all links, while the nested approach consists in swapping only at nodes 2 and 4.
	When swaps are probabilistic, the nested approach is generally optimal in terms of expected delivery time.}
	\label{fig.swapasap-vs-doubling}
\end{figure}

Optimal policies are also increasingly faster than swap-asap for increasing $n$, as shown in Figure~\ref{fig.advantage-n}.
For example, for $p=0.3$, $p_\text{s}=0.5$, and $t_\text{cut}=2$, the relative difference in expected delivery time is 1.7\%, 5.9\%, and 12.3\%, for $n=4$, 5, and 6, respectively.
This is in line with the fact that, when the number of nodes grows, there are increasingly more states in which the optimal action is to perform a strict subset of all possible swaps, as shown in Appendix~\ref{app.swapasap-states}.
Note that, in three- and four-node chains, the relative difference in expected delivery time is generally below 1\%.

\begin{figure}[t]
\captionsetup[subfigure]{justification=centering}
     \centering
     \begin{subfigure}[b]{0.8\textwidth}
         \centering
         \includegraphics[width=\textwidth]{advantageVSn_p0.3-0.9_ps1_tcut2_randseed2.pdf}
	\vspace{-1.5\baselineskip}
         \caption{Deterministic swaps ($p_\mathrm{s}=1$).}
	\label{fig.advantage-n-ps1}
         \vspace{5pt}
     \end{subfigure}

     \begin{subfigure}[b]{0.8\textwidth}
         \centering
         \includegraphics[width=\textwidth]{advantageVSn_p0.3-0.9_ps0.5_tcut2_randseed2.pdf}
	\vspace{-1.5\baselineskip}
         \caption{Probabilistic swaps ($p_\mathrm{s}=0.5$).}
	\label{fig.advantage-n-ps05}
         \vspace{5pt}
     \end{subfigure}

	\caption{\textbf{An optimal policy performs increasingly better than swap-asap in longer chains.}
	Relative difference between the expected delivery times of an optimal policy, $T_\mathrm{opt}$, and the swap-asap policy, $T_\mathrm{swap}$, for $\tc=2$ and different values of $p$, as a function of the number of nodes $n$.
	Black lines correspond to $p=0.3$, and the value of $p$ increases in steps of $0.1$ with increasing line transparency up to $p=0.9$.
	}
	\label{fig.advantage-n}
\end{figure}

\subsection{Network model}\label{subsec.model}
We analyze quantum repeater chains wherein nodes can store quantum states in the form of qubits and can perform three basic operations with them:
\\

\ssec{Entanglement generation.} Two adjacent nodes can attempt the heralded generation of an entangled link (i.e., a shared bipartite entangled state), succeeding with probability $p$. Generation of entanglement is heralded, meaning that the nodes receive a message stating whether they successfully generated an entangled link or not \cite{Barrett2005,Bernien2013}.
We assume that entanglement generation is noisy. Hence, the newly generated entangled links are not maximally entangled states but Werner states \cite{Werner1989}. Werner states are maximally entangled states that have been subjected to a depolarizing process, which is a worst-case noise model \cite{Dur2005}, and they can be written as follows:
\begin{equation}
    \rho = \frac{4F-1}{3}\ketbra{\phi^+} + \frac{1-F}{3} \mathbb{I}_4,
\end{equation}
where $\ket{\phi^+}=\frac{\ket{00}+\ket{11}}{\sqrt{2}}$ is a maximally entangled state, $F$ is the fidelity of the Werner state to the state $\ket{\phi^+}$, and $\mathbb{I}_d$ is the $d$-dimensional identity.
In our notation, the fidelity of a mixed state $\rho$ to a pure state $\ket{\phi}$ is defined as
\begin{equation}
    F(\rho,\ket{\phi}) \vcentcolon= \bra{\phi}\rho\ket{\phi}.
\end{equation}
We assume that the fidelity of newly generated entangled links is $F_\text{new}\leq1$.
\\

\ssec{Entanglement swap.} Two neighboring entangled links can be fused into a longer-distance entangled link. Consider a situation where node B shares an entangled link with node A, and another link with node C (see Figure~\ref{fig.swap}). Then, B can perform an entanglement swap to produce an entangled link between A and C while consuming both initial links \cite{Zukowski1993, Duan2001, Sangouard2011}. We refer to the link generated in a swap operation as a \emph{swapped link}.
This operation is also probabilistic: a new link is produced with probability $p_\text{s}$, and no link is produced (but both input links are still consumed) with probability $1-p_\text{s}$.

The generation of an entangled link between two end nodes without intermediate repeaters is limited by the distance between the end nodes \cite{Munro2015} -- e.g., the noise affecting a photon sent over an optical fiber grows exponentially with the length of the fiber \cite{Briegel1998}.
Therefore, a repeater chain that makes use of entanglement swapping is needed to generate end-to-end entanglement over long distances.
\\

\begin{figure}[t]
	\centering
	\includegraphics[width=0.6\textwidth]{figs/swap}
	\caption{\textbf{Entanglement swap.}
	When node B performs a swap, an entangled link between nodes A and B and an entangled link between nodes B and C are consumed to produce a single entangled link between A and C. This operation is essential for the generation of long-distance entanglement.}
	\label{fig.swap}
\end{figure}

\ssec{Cutoffs.} The fidelity of a quantum state decreases over time due to couplings to the environment \cite{Dur2005,Chirolli2008}.
These decoherence processes can be captured using a white noise model in which a depolarizing channel is applied to the entangled state at every instant. As a result, the fidelity of a Werner state at time $t$, $F(t)$, is given by
\begin{equation}\label{eq.ftT}
    F(t) = \frac{1}{4} + \Big( F(t-\Delta t) - \frac{1}{4}\Big) e^{-\Delta t/T_2},
\end{equation}
where $\Delta t$ is an arbitrary interval of time and $T_2$ is the coherence time of the qubits, a parameter that depends on the physical realization of the qubit. (\ref{eq.ftT}) is derived in Appendix \ref{app.decoherence}.

In general, quantum network applications require quantum states with fidelity above some threshold value $F_\text{min}$.
A common solution is to impose a cutoff time $\tc$ on the entangled links: all entangled links used to generate the final end-to-end link must be generated within a time window of size $\tc$ \cite{Rozpedek2018}.
Imposing memory cutoffs requires keeping track of the time passed since the creation of each entangled link. We call this time the \emph{age} of the link.
A link is discarded whenever it gets older than $\tc$.
Moreover, we assume that an entangled link generated as a result of entanglement swapping assumes the age of the oldest link that was involved in the swapping operation.
Another valid approach to calculate the age of a swapped link would be to re-compute the age based on the post-swap fidelity, although this would lead to a more complicated formulation to ensure that all the links that were used to produce a swapped link were generated within the time window of size $t_\text{cut}$.
To produce end-to-end links with fidelity above $F_\text{min}$ on a repeater chain that generates new links with fidelity $F_\text{new}$, it suffices to ensure that the worst policy, i.e., the policy that produces the lowest end-to-end fidelity, satisfies this requirement.
In Appendix \ref{app.cutoff}, we show that the policy that yields the lowest fidelity for the end-to-end link waits until the end of the $t_\text{cut}$ interval to perform all entanglement swaps.
Analyzing such a worst policy leads to the following condition for the cutoff time:
\begin{equation}\label{eq.cutoffcondition}
    \tc \leq -T_2 \ln\Bigg(\frac{3}{4F_\text{new}-1} \Big( \frac{4F_\text{min}-1}{3} \Big)^{\frac{1}{n-1}} \Bigg),
\end{equation}
where $n$ is the number of nodes.
For a full derivation of the previous condition, see Appendix \ref{app.cutoff}.
\\

In this paper, we consider quantum networks that operate with a limited number of qubits. Specifically, we use the following additional assumptions:
\begin{enumerate}[label=(\roman*)]
	\item The chain is \textbf{homogeneous}, i.e., the hardware is identical in all nodes. This means that all pairs of neighbors generate links with the same success probability $p$ and fidelity $F_\text{new}$, all swaps succeed with probability $p_\text{s}$, all states decohere according to some coherence time $T_2$, and all nodes apply the same cutoff time $t_\text{cut}$. This assumption may not hold for some long-distance quantum networks where each node is implemented using a different technology, but may be directly applicable to, e.g., small metropolitan-scale networks.
	\item We assume that each node has only \textbf{two storage qubits}, each of which is used to generate entanglement with one side of the chain. Each end node has a single storage qubit.
	This assumption is in line with the expectations for early quantum networks, in which nodes are likely to have a number of storage qubits on the order of the unit (e.g., in \cite{Pompili2021} the authors realized the first three-node quantum network using NV centers, each with a single storage qubit).
    	\item Time is discretized into non-overlapping \textbf{time slots}. During one time slot, the nodes first attempt entanglement generation; second, they can attempt entanglement swaps; and lastly they discard any entangled link that existed for longer than $\tc$ time slots.
    The unit of time used in this paper is the duration of a time slot, unless otherwise specified.
    	\item We also assume that classical communication between nodes is instantaneous. This means that every node has \textbf{global knowledge} of the state of the repeater chain in real time.
	In general, this is not a realistic assumption. However, given that classical communication delays decrease the performance of the network, our results constitute a lower bound on the expected delivery time of real setups and can be used as a benchmark.
\end{enumerate}

A repeater chain under the previous assumptions is characterized by four parameters:
\begin{itemize}
	\item[·] $n$: number of nodes in the chain, including end nodes.
	\item[·] $p$: probability of successful entanglement generation.
	\item[·] $p_\text{s}$: probability of successful swap.
	\item[·] $\tc$: cutoff time.
	Note that $F_\text{new}$, $F_\text{min}$, and $T_2$ are used to determine a proper value of cutoff time (see condition (\ref{eq.cutoffcondition})), but they are not needed after that.
\end{itemize}

In an experimental setup, the value of $p$ is determined by the inter-node distance and the type of hardware used, as quantum nodes can be realized using different technologies, such as NV centers \cite{Rozpedek2019, Bernien2013, Hensen2015, Humphreys2018, Pompili2021} and trapped ions \cite{Moehring2007,Slodicka2013}.
Linear optics setups generally perform swaps with probability $p_\text{s}=0.5$ \cite{Duan2001,Calsamiglia2001}, while other setups can perform deterministic swaps ($p_\text{s}=1$) at the cost of a slower speed of operation \cite{Pompili2021}.
The cutoff time $\tc$ can be chosen by the user, as long as condition (\ref{eq.cutoffcondition}) is satisfied. Note that (\ref{eq.cutoffcondition}) depends on $T_2$ (which depends on the hardware available), $F_\text{new}$ (which depends on the hardware and the choice of entanglement generation protocol), and $F_\text{min}$ (which is specified by the final application).
\\

The state of the repeater chain at the end of each time slot can be described using the age of every entangled link.
In Figure~\ref{fig.chain_evolution} we show an example of the evolution of the state of a chain with cutoff $t_\text{cut}=3$, over four time slots:
\begin{itemize}
	\item In the first time slot ($t\in [0,1)$), all pairs of neighbors attempt entanglement generation, but it only succeeds between nodes two and three.
	No swaps can be performed, and the only link present is younger than the cutoff, so it is not discarded.
	\item In the second time slot ($t\in [1,2)$), the age of the link between nodes two and three increases by one. All pairs of neighbors (except nodes two and three) attempt entanglement generation, which succeeds between nodes four and five.
	\item In the third time slot ($t\in [2,3)$), the age of both existing links increases by one.
	All pairs of neighbors (except nodes two and three and nodes four and five) attempt entanglement generation, and only nodes five and six succeed.
	A swap can be performed at node five but they decide to wait.
	\item In the fourth time slot ($t\in [3,4)$), the age of every existing link increases by one.
	Nodes one and two and nodes three and four attempt entanglement generation but none of the pairs succeeds.
	A swap is successfully performed at node five, and a new link between nodes four and six is generated. This new link assumes the age of the oldest link involved in the swap operation.
	Lastly, the entangled link between nodes two and three is discarded, as its age reached the cutoff time.
\end{itemize}

\begin{figure}[t]
	\centering
	\includegraphics[width=0.8\textwidth]{figs/chain_evolution}
	\caption{\textbf{Example of entangled link dynamics in a repeater chain.}
	Each row represents the state of the chain at the end of time slot $t$.
	Entangled links are represented as black solid lines, with occupied qubits as black circles. The number above each entangled link is the age of the link. We assume cutoff $t_\text{cut}=3$.}
	\label{fig.chain_evolution}
\end{figure}


\subsection{Expansion of BH thermodynamics} The original first law of BH thermodynamics, that exploits the quanutm relation between temperature \emph{T} and surface gravity  $\kappa$ and that of entropy \emph{S} and area of the black hole, for a charged rotating BH has the following form :
\begin{equation} \label{eqn : 1}
    \delta M = T\delta S + \Omega\delta J + \phi \delta Q \ ,
\end{equation}
where \emph{M} is the mass identified with the internal energy \emph{U}, $\Omega$ and \emph{J} is angular velocity and angular momentum associated with rotating black hole and $\phi$ and \emph{Q} is the electric potential and charge associated with the charged BH. But equation \ref{eqn : 1} does not yet include the \emph{pressure} and \emph{volume} terms which are commonly found in thermodynamics. It was by then well established that the cosmological constant $\Lambda$ can be treated as a thermodynamical BH parameter \cite{16,17}. Consequently, the \emph{pressure} of the BH was interpreted in terms of a negative cosmological constant and the corresponding $P \delta V$ term was finally added to the original first law  \cite{18}, leading to  a more generalized first law of the following form :
\begin{equation}
    \delta M = T\delta S + V \delta P + \Omega\delta J + \phi \delta Q\ ,
\end{equation} where \begin{equation*}
 P = \frac{-\Lambda}{8\pi}
\end{equation*} and the corresponding conjugate variable volume is \cite{19,20},
\begin{equation*}
V \equiv \left( \frac{\partial M}{\partial P}\right)_{S,Q,J}
\end{equation*} This also changed the interpretation of mass \emph{M} to the enthalpy \emph{H} \cite{21} of the BH, which now involves both the internal energy \emph{U} and energy due to the \emph{PV} term. Expansion of BH thermodynamics in such a way leads to a new extended picture of the phase space that allows us to study phase transitions with a renewed perspective. Such an investigation has led to an enhanced identification between the AdS black holes and van der Waals liquid-gas system \cite{22,23} and other interesting discoveries \cite{24,25,26,27,28}.
\subsection{Holographic heat engines}
With this kind of framework at hand, Johnson put forth the idea of "holographic heat engines" \cite{29}, wherein he utilises the \emph{pdV} term to obtain mechanical work from heat energy, black holes being the working substance. For such a thermodynamic system, work can be defined in terms of a cycle in the \emph{P - V} plane with \emph{$Q_H$} and \emph{$Q_C$} being the net input and output heat respectively and the net mechanical work is written as \emph{W} =   \emph{$Q_H$} - \emph{$Q_C$}. The heat engine's efficiency is defined in the usual thermodynamical sense as $\eta = W/Q_H = 1 - Q_C / Q_H$. As we know, the efficiency of Carnot's engine corresponds with the maximum achievable efficiency given by,
\begin{equation*}
    \eta_c = 1 - \frac{Q_C}{Q_H} = 1 - \frac{T_C}{T_H}
\end{equation*} where $T_C$ and $T_H$ are lower and higher temperature of the reservoir respectively. \par In this paper, we develop a more generalised version of thermodynamic heat engine by incorporating the idea of massive gravity and study the effects of such a modification on the efficiency of the heat engine. There are several arguments one can make in support of such a generalisation : The theory of general relativity (GR) is known to describe the gravitational force with a good agreement with the observations. However the search for alternative theories has been an on-going quest since the formulation of GR. From a modern particle physics
perspective, GR is a unique theory for a massless spin-2 particle. Taking GR as starting point, alternative theories intend to modify or extend GR to different scales in a way as to account for the observational anomalies. One such modification maintains the notion of gravity being propagated by a spin-2 particle while giving a mass to this particle, that is the so called "Massive gravity". The first attempt to add a mass to the graviton was made by Fierz and Pauli \cite{30}, who in their work showed the presence of an extra degree of freedom (called ghost) in a generic Lorentz-invariant modification of gravity. Unfortunately, presence of an unwanted mode was discovered, named as the vDVZ discontinuity \cite{31,32} that represented failure of a massive gravity theory to reduce itself into general relativity (a massless spin-2 field theory with two propagating degree of freedom ) in the massless limit. This pathology naturally posed a threat to one of the most important underlying principle when it comes to deforming a theory : Continuity in physical predictions of a theory in the parameters of the theory. This issue was then resolved by Vainshtein in \cite{33} by proposing a mechanism through which nonlinearities managed to hid the extra degrees of freedom, at short distances and has been successfully applied to several models \cite{34,35,36}.
\par Thus, to summarize, the cosmological models motivated by GR are plagued by issues like the
cosmological constant problem, late time acceleration of universe, issues at the interface of gravity and quantum mechanics, etc., forcing us to look for modifications such as mentioned above. Massive gravity theories are of immense interest as an alternative to GR, strongly motivated on both theoretical and phenomenological grounds \cite{37,38,39}.
 Undoubtedly, its inclusion in the study of the efficiency of BH heat engines would result into some interesting insights. Also there are enough evidences of such impact of a massive graviton, for example the presence of van der Waals phase transitions for non-spherical geometries, in holography \cite{40}, in the resolution to the cosmological constant problem \cite{41,42}, for enhanced polarization of gravitational waves \cite{43} and many more. It would be interesting to find out the effect of such highly influential graviton mass on the efficacy of the heat engines. We consider a regular charged black hole (RCB) solution for such an investigationand and for additional work along this direction for different black hole systems see\cite{44,45,46}. In doing so, we also calculate temperature, pressure and specific heat of the RCB heat engine. \par The paper is organised as follows : In the next section, we begin by constructing a RCB heat engine in massive gravity in extended phase space and draw some important results through their plots. In section \ref{section :3} we turn to our main aim that is to study the effect of inclusion of graviton mass on the efficiency of the RCB heat engine. In the last section, we summarize our results.

\section{Introduction}
\input{Introduction}
\section{Regular charged black holes in massive gravity as heat engines}
\input{RCB}
\section{Efficiency of the RCB heat engine}
\input{effiiciency}
\newpage
\section{Conclusion}
\input{Conclusions}


\section{Introduction}
Neural Ordinary Differential Equations (ODEs) were first introduced in 2018, as a radical new neural network design that boasted better memory efficiency, and an ability to deal with irregularly sampled data \cite{hao_2020}. The idea behind this family of deep learning models is that instead of specifying a discrete sequence of hidden layers, we instead parameterize the derivative of the hidden states using a neural network \cite{chen2018neural}. The output of the network can then be computed using a differential equation solver \cite{chen2018neural}. This work has spurred a whole range of follow-up work, and since 2018 several variants have been proposed, such as augmented neural ODEs (ANODEs) and their ensuing variants \cite{norcliffe2020sonode,dupont2019neurips,amir2019anode}. These variants provide a more expressive formalism by augmenting the state space of neural ODEs to allow for the flow of state trajectories to cross. This crossing, prohibited in the original framework, allows for the learning of more complex functions that were prohibited by the original neural ODE formulation \cite{dupont2019neurips}.

Due to the potential that neural networks boast in revolutionizing the development of intelligent systems in numerous domains, the last several years have witnessed a significant amount of work towards the formal analysis of these models.
The first set of approaches that were developed considered the formal verification of neural networks (NN), using a variety of techniques including reachability methods \cite{tran2020cavtool,wang2018reluval,bak2021nnenum,ruan2018ijcai}, and SAT techniques \cite{katz2017reluplex,katz2019marabou,ehlers2017atva}. Thereafter, many researchers proposed novel formal method approaches for neural network control systems (NNCS), where the majority of methods utilized a combination of NN and hybrid system verification techniques \cite{fan2020reachNN*,Huang2019reachNN,tran2019emsoft,Ivanov2019verisig,ivanov2021verisig20,bogomolov2019Juliareach}. Building on this work, a natural outgrowth is extending these approaches to analyze and verify neural ODEs,
and some recent studies have considered the analysis of formal properties of neural ODEs. One such study aims to improve the understanding of the inner operation of these networks by analyzing and experimenting with multiple neural ODE architectures on different benchmarks \cite{massaroli2020neurips}. Some studies have considered analyses of the robustness of neural ODEs such as \cite{carrara2019wifs} and \cite{yan2020robustness}, which evaluate the robustness of image classification neural ODEs and compare the efficacy of this class of network against other more traditional image classifier architectures. The first reachability technique targeted for neural ODEs presented a theoretical regime for verifying neural ODEs using Stochastic Lagrangian Reachability (SLR) \cite{gruenbacher2021verification}. This method is an abstraction-based technique that computes confidence intervals for the calculated reachable set with probabilistic guarantees. In a follow-up work, these methods were improved and implemented in a tool called Gotube \cite{gruenbacher2021gotube}, which is able to compute reach sets for longer time horizons than most state-of-the-art reachability tools. However, these methods only provide stochastic bounds on the reach sets, so there are no formal guarantees on the derived results.

To the best of our knowledge, this paper presents the first deterministic verification framework for a general class of neural ODEs with multiple continuous-time and discrete-time layers. In this work, we present our verification framework NNVODE that makes use of deterministic reachability approaches for the analysis of neural ODEs.
Our methods are evaluated on a set of benchmarks with different architectures and conditions in the area of dynamical systems, control systems, and image classification. We also compare our results against three state-of-the-art verification tools when possible.
In summary, the contributions of this paper are:
\begin{itemize}
    \item We introduce a general class of neural ODEs that allows for the combination of multiple continuous-time and discrete-time layers.
    \item We develop NNVODE, an extension of NNV \cite{tran2020cavtool}, to formally analyze a general class of neural ODEs using sound and deterministic reachable methods.

    \item We run an extensive evaluation on a collection of benchmarks within the context of time-series analysis, control systems, and image classification. We compare the results to Flow*, GoTube, and JuliaReach for neural ODE architectures where this is possible.

\end{itemize}

\section{Background and Problem Formulation}

Neural ODEs emerged as a continuous-depth variant of neural networks becoming a special case of ordinary differential equations (ODEs), where the derivatives are defined by a neural network, and can be expressed as:
\begin{equation}
\label{eq:omegalayer_func}
    \dot{z} = g(z), \quad \quad \quad  z(t_0) = z_0,
\end{equation}
where the dynamics of the function $g:\mathbb{R}^j\rightarrow{}\mathbb{R}^p$ are represented as a neural network, and initial state $z_0 \in \mathbb{R}^{j}$, where $j$ corresponds to the dimensionality of the states $z$. A \emph{neural network} (NN) is defined to be a collection of consecutively connected $\mathscr{NN}$-Layers as described in Definition \ref{def:nnLayer}.

\begin{definition}[$\mathscr{NN}$-Layer]
\label{def:nnLayer}
    A \textbf{$\mathscr{NN}$-Layer} is a function h: $\mathbb{R}^{j} \rightarrow{} \mathbb{R}^{p}$, with input x $\in$ $\mathbb{R}^{j}$, output y $\in$ $\mathbb{R}^{p}$

    defined as follows
    \begin{gather}
        y = h(x)
    \end{gather}
\end{definition}
where the function $h$ is determined by parameters $\theta$, typically defined as a tuple $\theta$ = $\langle$$\sigma$,\textbf{W},\textbf{b}$\rangle$ for fully-connected layers, where \textbf{W} $\in$ $\mathbb{R}^{j\times p}$, \textbf{b} $\in$ $\mathbb{R}^{p}$, and activation function $\sigma$: $\mathbb{R}^{j} \rightarrow{} \mathbb{R}^{p}$, thus the fully-connected $\mathscr{NN}$-Layer is described as
\begin{gather}\label{eq:nnlayer}
    y = h(x) = \sigma(\textbf{W}(x) + \textbf{b}).
\end{gather}

However, for other layers such as convolutional-type $\mathscr{NN}$-Layers, $\theta$ may include parameters like the filter size, padding, or dilation factor and the function $h$ in ~\eqref{eq:nnlayer} may not necessarily apply. For a formal definition and description of the reachability analysis for each of these layers that are integrated within NNVODE, we refer the reader to Section 4 of \cite{tran2020cav}. For the Neural ODEs (NODEs),
we assume that $g$ is Lipschitz continuous, which guarantees that the solution of $\dot{z}=g(z)$ exists. This assumption allows us to model $g$ with $m$ layers of continuously differentiable activation functions $\sigma_k$ for $k\in\{1,2,...,m\}$ such as sigmoid, tanh, and exponential activation functions. Described in ~\eqref{eq:omegalayer_func} is the notion of a NODE as introduced in \cite{chen2018neural} and an example is illustrated in Figure \ref{fig:omegalayer}.

\begin{definition}[NODE]
\label{def:omegaLayer}
    A \textbf{NODE} is a function $\dot{z} = g(z)$ with $m$ fully-connected $\mathscr{NN}$-Layers, and it is defined as follows
    \begin{gather}\label{eq:omegalayer_layer}
        \dot{z} = g(z) = h_m(h_{m-1}(...h_1(z))),
    \end{gather}
    where $g : \mathbb{R}^j \rightarrow{} \mathbb{R}^p$, $\sigma_k : \mathbb{R}^{j_k} \rightarrow{} \mathbb{R}^{p_k}$, $\textbf{W}_k \in \mathbb{R}^{j_k \times p_k}$, and $\textbf{b}_k \in \mathbb{R}^{p_k}$. For each layer $k$=$\{1,2,\dots,m\}$, we describe the function $h_k$ as in ~\eqref{eq:nnlayer}.
\end{definition}

\begin{figure}[!ht]
  \centering
    \includegraphics[trim={300 170 300 230}, clip, width=0.55\columnwidth]{figures/Omegalayer.pdf}
   \caption{Illustration of an example of NODE, as defined in \eqref{eq:omegalayer_func} and ~\eqref{eq:omegalayer_layer}, and Definition \ref{def:omegaLayer}.}
  \label{fig:omegalayer}
\end{figure}

An example of NODE with m~=~2 hidden layers, 2 inputs and 2 outputs is depicted in Figure~\ref{fig:omegalayer}.

\subsection{General Neural ODE}
The general class of neural ODEs (GNODEs) considered in this work is more complex than previously analyzed neural ODEs as it may be comprised of two types of layer: NODEs and $\mathscr{NN}$-Layers.
We introduce a more general framework where multiple NODEs can make up part of the overall architecture along with other $\mathscr{NN}$-Layers, as described in Definition \ref{def:neuralODE}. This is the reachability problem subject of evaluation in this work.

\begin{definition}[GNODE] \label{def:neuralODE}
A \textbf{GNODE} $~\mathcal{F}$ is any sequence of consecutively connected $N$ layers $\mathcal{L}_k$ for $k \in \{1, \dots N\}$ with $N_O$ NODEs and $N_{D}$ $\mathscr{NN}$-Layers, that meets the conditions $1 \leq N_O \leq N$, $0 \leq N_{D} < N$, and $N_D$ + $N_O$ = $N$.
\end{definition}

With the above definition, we formulate a theorem for the restricted class of neural ODEs (NODE) that we use to compare our methods against existing techniques.

\begin{observation}[Special case 1: NODE]\label{obs:omegalayer}
Let $\mathcal{F}$ be a GNODE with $N$ layers. If $N = 1$, $N_O = 1$ and $N_D = 0$, then $\mathcal{F}$ is equivalent to an ODE whose continuous-time dynamics are defined as a neural network, and we refer to it as a \textbf{NODE} (as in Definition~\ref{def:omegaLayer}).
\end{observation}

In Figure \ref{fig:gnode}, an example of a GNODE is shown, which has 1 input (\emph{x}), 1 output (\emph{y}), $N_O$ = 2 NODEs, and $N_D$ = 5 $\mathscr{NN}$-Layers, with its 5 numbered segments described as: 1) The first segment has a $\mathscr{NN}$-Layer, with one hidden layer of 2 neurons and an output layer of 2 neurons, 2) a NODE with one hidden layer of 3 neurons and an output layer of 2 neurons, 3) $\mathscr{NN}$-Layer with 3 hidden layers of 4,1, and 2 neurons respectively, and an output layer of 2 neurons, 4) NODE with a hidden layer of 3 neurons and output layer of 2 neurons, and 5) $\mathscr{NN}$-Layer with a hidden layer of 4 neurons and and output layer with 1 neurons.

\begin{figure}[!ht]
  \centering
    \includegraphics[trim=0 0 0 0, clip, width=0.95\columnwidth]{figures/gnode.png}
    \vspace{-1em}
  \caption{Example of a GNODE. The filled circles represent weighted neurons in each layer, non-filled neurons represent the inputs of each layer-type segment, shown for visualization purposes. }
  \label{fig:gnode}
\end{figure}

Given this GNODE architecture, we are capable of encoding the originally proposed neural ODE \cite{chen2018neural}, as well as several of its model improvements including higher-order neural ODEs like the second-order neural ODE (SONODE) \cite{norcliffe2020sonode}, augmented neural ODEs (ANODEs) \cite{dupont2019neurips}, and input-layer augmented neural ODEs (ILNODEs).
This general class of neural ODEs from Definition~\ref{def:neuralODE} is applicable to many applications including image classification, time series prediction, dynamical system modeling and normalizing flows.

\begin{observation}[Special case 2 - NNCS]\label{obs:nncs}
We can consider the NNCS as a special case of the \emph{GNODE}, as the NNCS models also consist of $\mathscr{NN}$-Layers and NODE, under the assumption that the plant is described as an NODE. NNCS in the general architecture have $\mathscr{NN}$-Layers followed by an ODE or NODE, which connects back to the $\mathscr{NN}$-Layers in a feedback loop manner for $cp$ control steps.
Thus, by unrolling the NNCS $cp$ times, we consecutively connect the $\mathscr{NN}$-Layers with the NODE, creating a \emph{GNODE}.
\end{observation}

\subsection{NODE: Applications} \label{sec:omegalayerapps}
There are two application modes of a NODE, model reduction and dynamical systems. On one hand, we can use it as a substitute for multiple similar layers to reduce the depth and overall size of a model \cite{chen2018neural}. In this context, we treat the NODE as an input-output mapping function, and set the integration time $t_f$ to a fixed number during training, which will be then used during simulation or reachability. Typically, this value is set to 1. On the other hand, we can make use of NODEs to capture the behavior of time-series data or dynamical systems. In this sense, $t_f$ is not fixed, and it will be determined by the user/application. In summary, we can use NODEs for: 1) time-series or dynamical system modeling, and 2) model reduction. For time-series, $t_f$ is variable, user or application dependent. For model reduction, $t_f$ is a parameter of the model fixed before training. Only the output value at time = $t_f$ is used, while for the time-series models, we are interested in the interval [0,$t_f$].

\subsection{Reachability Analysis} \label{sec:reachAnalysis}

The main focus of this manuscript is to introduce a framework for the reachability analysis of GNODEs. This framework
combines set representations and methods from Neural Network (NN), Convolutional Neural Network (CNN) and hybrid systems reachability analysis. Similar to neural network reachability in NNV \cite{tran2020cavtool}, we consider the set propagation through each layer as well as the set conversion between layers for all reachability methods for the supported layers. Hence, the reachability problem of GNODE is defined as follows:

\begin{definition}[Reachable Set of a GNODE]
    \label{def:NODEreach}
    Let $\mathcal{F}$: $\mathbb{R}^j \rightarrow{} \mathbb{R}^p$ be a \textbf{GNODE} with $N$ layers. The output reachable set $\mathcal{R}_{N}$ of a \textbf{GNODE} with input set $\mathcal{R}_0$ is defined as:
    \begin{equation}
    \begin{split}
        \mathcal{R}_{1} &\triangleq \{y_1~|~y_1 = f_1(y_0), ~ y_0 \in \mathcal{R}_0\}, \\
      \mathcal{R}_{2} &\triangleq \{y_2~|~y_2 = f_2(y_1), ~ y_1 \in \mathcal{R}_{1}\}, \\
     &\vdots \\
    \mathcal{R}_{N} &\triangleq \{y_N ~|~ y_N = f_N(y_{N-1}), ~y_{N-1} \in \mathcal{R}_{N-1} \}, \\
    \end{split}
    \end{equation}
where $f_i$: $\mathbb{R}^{j_i} \rightarrow{} \mathbb{R}^{p_i}$ is the function $f$ of layer $i$, where $f$ is either a NODE (g) or a $\mathscr{NN}$-Layer (h).
\end{definition}

The proposed solution to this problem is sound and incomplete, in other words, an over-approximation of the reachable set. This means that, given a set of inputs, we compute an output set of which, given any point in the input set, we simulate the GNODE and the output point is always contained within the reachable output set (sound). However, it is incomplete because the opposite is not true; there may be points in the output reachable set that cannot be traced back to be within the bounds of the input set.

\begin{definition}[Soundness] \label{def:sound}
Let  $\mathcal{F}$: $\mathbb{R}^j$ $\rightarrow{\mathbb{R}^p}$ be a GNODE with an input set $\mathcal{R}_0$ and output reachable set $\mathcal{R}_f$. The computed $\mathcal{R}_f$ given $\mathcal{F}$ and $\mathcal{R}_0$ is \textbf{sound} iff $\forall$$x$ $\in$ $\mathcal{R}_0$, $|$ $y$ = $\mathcal{F}$($x$), $y$ $\in$ $\mathcal{R}_f$.
\end{definition}

\begin{definition}[Completeness] \label{def:complete}
Let  $\mathcal{F}$: $\mathbb{R}^j$ $\rightarrow{\mathbb{R}^p}$ be a GNODE with an input set $\mathcal{R}_0$ and output reachable set $\mathcal{R}_f$. The computed $\mathcal{R}_f$ given $\mathcal{F}$ and $\mathcal{R}_0$ is \textbf{complete} iff $\forall$$x$ $\in$ $\mathcal{R}_0$, $\exists$$y$ = $\mathcal{F}$($x$) $|$ $y$ $\in$ $\mathcal{R}_f$ and  $\forall$$y$ $\in$ $\mathcal{R}_f$, $\exists$$x$ $\in$ $\mathcal{R}_0$ $|$ $y$ = $\mathcal{F}$($x$).
\end{definition}

\begin{definition}[Reachable set of a $\mathscr{NN}$-Layer] \label{def:reachlayer}
Let h: $\mathbb{R}^{j} \rightarrow{} \mathbb{R}^{p}$ be a $\mathscr{NN}-Layer$ as described in Definition \ref{def:nnLayer}. The reachable set $\mathcal{R}_{h}$, with input $\mathcal{X} \subset \mathbb{R}^{n}$ is defined as
    \begin{equation*}
    \centering
        \mathcal{R}_{h} = \{y~|~ y = h(x),  x \in \mathcal{X}\}.
    \end{equation*}
\end{definition}

Definition \ref{def:reachlayer}  applies to any discrete-time layer that is part of a GNODE, including, but not limited to the following supported $\mathscr{NN}-Layers$ supported in NNVODE: fully-connected layers with ReLU, tanh, sigmoid and leaky-ReLU activation functions, and convolutional-type layers such as batch normalization, 2-D convolutional, and max-pooling.

The reachability analysis of NODEs is akin to the general reachability problem for any continuous-time system modeled by an ODE.
If we represent the NODE as a single layer $i$ of a GNODE with continuous dynamics described by ~\eqref{eq:omegalayer_func}, and assume that for a given initial state $z_0 \in \mathbb{R}^{j_i}$, the system admits a unique trajectory defined on $\mathbb{R}^+_0$, described by $\zeta(., z(t_0))$, then the reachable set of the given NODE can be characterized by Definition \ref{def:ODEreach}.

\begin{definition}[Reachable set of a NODE] \label{def:ODEreach} Let $g$ be a \textbf{NODE} with solution $\zeta$(t; z$_0$) to ~\eqref{eq:omegalayer_func} for initial state z$_0$. The reachable set, $\mathcal{R}_{g}$ at $t = t_F$,  $\mathcal{R}_{g}(t_F)$, with initial set $\mathcal{R}_{0} \subset \mathbb{R}^{n}$ at time $t=t_0$ is defined as
    \begin{equation*}
    \centering
        \mathcal{R}_{g}(t_F) = \{ \zeta(t; z_0) \in \mathbb{R}^{n} ~|~ z_0 \in \mathcal{R}_{0}, t \in [0,t_f] \}.
    \end{equation*}
We also describe the reachable set for a time interval [$t_0$,$t_f$] as follows
    \begin{equation*}
        \mathcal{R}_{g}([t_0,t_f]) := \bigcup_{t\in[t_0,t_f]}\mathcal{R}_{g}(t)
    \end{equation*}
\end{definition}

Now that we have outlined the general reachability problem of a NODE, whether we compute a single reachable set for the NODE at $t$ = $t_F$, or over an interval $t$ $\in$ [$t_0$,$t_F$] (computed by \emph{get$\_$time()}), depends on how the neural ODE was trained and the specific application of its use. However, the core computation remains the same. This is outlined in Algorithm \ref{alg:reach}.

\begin{algorithm}[h]
	\caption{Reachability analysis of a GNODE.}\label{alg:reach}
	\begin{algorithmic}[1]
        \small
		\Input{$\mathcal{F}, \mathcal{R}_0$}  ~{\scriptsize \textcolor{blue}{ $//$ GNODE, input set}}
		\Output{$\mathcal{R}_f$} ~{\scriptsize \textcolor{blue}{$//$ output reachable set}}
		\Procedure{\texttt{$\mathcal{R}_f$ = reach}}{$\mathcal{F}, \mathcal{R}_0$}
		\State N = $\mathcal{F}$.layers ~{\scriptsize \textcolor{blue}{ $//$ number of layers}}
		\For{\textit{i = 1 : N}} ~{\scriptsize \textcolor{blue}{ $//$ loop through every layer}}
		\State $\mathscr{L}_i$ = $\mathcal{F}$.layer(i)
		\If{$\mathscr{L}_i$ $is$ NODE }~{\scriptsize \textcolor{blue}{ $//$ check layer type}}
		\State $\textbf{t}_i$ = get$\_$time($\mathscr{L}_i$) ~{\scriptsize \textcolor{blue}{ $//$ get integration time bounds of layer $i$}}
		\State $\mathcal{R}_i$ =  reach$_{NODE}$($\mathscr{L}_i$,$\mathcal{R}_{i-1}$, $\textbf{t}_i$) ~{\scriptsize \textcolor{blue}{ $//$ reach set of layer $i$}, Definition \ref{def:ODEreach}}
		\Else \State $\mathcal{R}_i$ =  reach$_{\mathscr{NN}}$($\mathscr{L}_i$, $\mathcal{R}_{i-1}$) ~{\scriptsize \textcolor{blue}{ $//$ reach set of layer $i$}, Definition \ref{def:reachlayer}}
		\EndIf
        \EndFor
        \State $\mathcal{R}_f$ = $\mathcal{R}_N$ ~{\scriptsize \textcolor{blue}{ $//$ output reachable set}}
		\EndProcedure
              \end{algorithmic}
\end{algorithm}

\subsection{Reachability Methods}

In the previous sections, we defined the reachable set of a GNODE using a layer-by-layer approach. However, the computation of these reachable sets is defined by the reachability methods and set representations utilized in their construction. For instance, the same operations are not utilized to compute the reachable set for a fully-connected layer with a hyperbolic tangent activation function as with a fully convolutional layer. In the following section, we describe the set of methods in the NNVODE tool available to the community to compute the reachable set of each specific layer.

We begin with the NODE approaches, where we make a distinction based on the underlying dynamics. If the NODE is nonlinear, we make use of zonotope and polynomial-zonotope based methods that are implemented and available in CORA \cite{Althoff2015arch,althoff2013hscc}. If the NODE is purely linear, then we utilize the star set-based methods introduced in \cite{bak2017cav}, which are more scalable than other zonotope-based methods and possess soundness guarantees as well.

For the $\mathscr{NN}$-Layers, there are several methods available, including zonotope-based and star-set based methods. However, we limited our implementation only using star sets to handle these layers, as this representation was demonstrated to be more computationally efficient and to enable tighter over-approximations of the derived reachable sets than zonotope methods \cite{tran2019fm}. Additionally, in using star set methods, we allow for the use of both approximate methods (\emph{approx-star}) and exact methods (\emph{exact-star}).
A summary of these methods and supported layers is depicted in Table \ref{tab:LayerSupport}, and for a complete description of the reachability methods utilized in our work, we refer the reader to \cite{tran2020cavtool} and the manual\footnote{NNV manual is available at: \url{https://github.com/verivital/nnv/blob/master/docs/manual.pdf}}\footnote{CORA manual (Release 2021) is available at: \url{https://tumcps.github.io/CORA/data/Cora2021Manual.pdf}}.

\begin{table}[h!]
    \caption{Layers supported in NNVODE and reachability sets and methods available.}
    \label{tab:LayerSupport}
    \centering
    \begin{adjustbox}{width=0.9\columnwidth, center}
    \begin{tabular}{c|l}

    \centering
     & \textbf{Layer Type} -- \textbf{Set Rep. (method name)} \\
    \hline
    NODE & Linear -- Star-set ("direct") \cite{bak2017cav} \\
    NODE & Nonlinear -- Zonotope, Polynomial Zonotope * \cite{Althoff2015arch,althoff2013hscc} \\

    $\mathscr{NN}$-Layer & FC: linear, ReLU -- Star-set ("approx-star", "exact-star") \cite{tran2019fm} \\
    $\mathscr{NN}$-Layer & FC: leakyReLU, tanh, sigmoid, satlin -- Star-set ("approx-star") \cite{tran2020cavtool} \\
    $\mathscr{NN}$-Layer & Conv2D -- ImageStar ("approx-star", "exact-star") \cite{tran2020cavtool} \\
    $\mathscr{NN}$-Layer & BatchNorm -- ImageStar ("approx-star", "exact-star") \cite{tran2020cavtool} \\
    $\mathscr{NN}$-Layer & MaxPooling2D -- ImageStar ("approx-star", "exact-star") \cite{tran2020cavtool} \\
    $\mathscr{NN}$-Layer & AvgPooling2D -- ImageStar ("approx-star", "exact-star") \cite{tran2020cavtool} \\
    \hline
    \end{tabular}
    \end{adjustbox}
    \begin{adjustwidth}{0.2cm}{}\scriptsize{* We support several methods available using Zonotope and PolyZonotopes, which includes user-defined fixed reachability parameters ("ZonoF", "PolyF") as well as adaptive reachability methods ("ZonoA" and "PolyA"), which require no prior knowledge on reach methods or systems to verify to produce relevant results.}
    \end{adjustwidth}
\end{table}

\noindent\textbf{Implementation.}
One of the key aspects of the verification of GNODEs is the proper encoding of the NODEs within reachability schemes. Depending on the software that is used, this process may vary and require distinct steps. As an example, some tools, like NNVODE, are simpler and allow for matrix multiplications within the definition of equations. However, for tools like Flow*, the set of steps required to properly encode this problem are more complex as it requires a definition for each individual equation of the state derivative. Thus, a more general conversion is needed, which is illustrated in the Appendix.

\section{Evaluation}\label{sec:eval}

Having described the details of our reachability definitions, algorithm, and implementation, we now present the experimental evaluation of our proposed work.
We begin by presenting, a method and tool comparison analysis against
GoTube\footnote{GoTube can be found at \url{https://github.com/DatenVorsprung/GoTube}} \cite{gruenbacher2021gotube}, Flow*\footnote{Flowstar version 2.1.0 is available at \url{https://flowstar.org/}} \cite{Chen2013Flow*} and JuliaReach\footnote{JuliaReach can be found at \url{https://juliareach.github.io/}} \cite{bogomolov2019Juliareach}.
Then, we present a case study of an Adaptive Cruise Control system, and conclude with an evaluation of the scalability of our techniques using a random set of architectures for dynamical system applications as well as a set of classification models for MNIST. The GNODE architectures for each benchmark can be found in the Appendix. To facilitate the reproducibility of our experiments, we set a timeout of 2 hours (7200 seconds) for each reach set computation\footnote{Code to reproduce all results can be found here: \url{https://github.com/verivital/nnv/tree/master/code/nnv/examples/Submission/FORMATS2022}}. All our experiments were conducted on a desktop with the following configuration: Intel Core i7-7700 CPU @ 3.6GHz 8 core Processor, 64 GB Memory, and 64-bit Ubuntu 16.04.3 LTS OS.

\subsection{Method and Tool Comparison}
We have implemented several methods within NNVODE, and the first evaluation consists of comparing the available methods for nonlinear NODEs, fixed-step zonotope and polynomial zonotopes (zono-F,poly-F) and adaptive zonotope and polynomial zonotope (zono-A,poly-A) based methods \cite{Althoff2015arch}. For all other $\mathscr{NN}$-Layers in the GNODEs, we use the star-set over-approximate methods.
We considered multiple models, all inspired by the ILNODE representation that was introduced by Massaroli. They consist of a set of models with a varying number of augmented dimensions.
All these models are instances of the GNODE class presented in Definition \ref{def:neuralODE}, which present an architecture of the form $\mathscr{NN}$-Layers + NODE + $\mathscr{NN}$-Layers. In this context, we were concerned with how the methods scale with respect to the number of dimensions of each model.

\vspace{-1em}
\begin{table}[h!]
    \caption{Computation time of the reachability analysis of the Damped Oscillator benchmark. Results are shown in seconds with up to one decimal place. }
    \label{tab:doRes}

    \begin{adjustbox}{width=0.5\columnwidth,center}
    \begin{tabular}{cccccc}
    \toprule

    Aug. Dims & Zono-\emph{F} & Zono-\emph{A} & Poly-\emph{F} & Poly-\emph{A} \\
    \hline
    0 & \textbf{34.0} & 574.5 & 201.7 & 654.0  \\
    1 & \textbf{146.4} & 4205.0 & 1573.9 & 3440.9 \\
    2 & \textbf{441.0} & -- & -- & -- \\
    \hline
    \end{tabular}
    \end{adjustbox}%
\end{table}
\vspace{-4em}
\begin{figure}[!ht]
  \centering
    \subfigure[Aug. Dims = 0]{

        \includegraphics[trim=95 230 115 230, clip, width=0.23\columnwidth]{figures/dampedOsc_compare0.pdf}\label{fig:osc_0}}
    \subfigure[Aug. Dims = 0]{

        \includegraphics[trim=98 230 120 230, clip, width=0.23\columnwidth]{figures/dampedOsc_compare0_zoomed.pdf}
        \label{fig:osc_0z}}
  \subfigure[Aug. Dims = 1]{

        \includegraphics[trim=98 230 115 240, clip, width=0.23\columnwidth]{figures/dampedOsc_compare1.pdf}\label{fig:osc_1}}
    \subfigure[Aug. Dims = 1]{

        \includegraphics[trim=98 230 100 230, clip, width=0.23\columnwidth]{figures/dampedOsc_compare1_zoomed.pdf}
        \label{fig:osc_1z}}
    \vspace{-1em}
    \caption{Reach sets comparison of the Damped Oscillator benchmark of the model with 0 and 1 augmented dimensions. Plots (a) and (c) show the reachable set at $t$ = 1s, and plots b) and d) show the zoomed-in reach sets to observe the minor size differences between the four methods: \textbf{\textcolor{red}{ZonoA}}, \textbf{ZonoF}, \textbf{\textcolor{blue}{PolyF}} and \textbf{\textcolor{Green}{PolyA}}.}
  \label{fig:damOsc_compare}
\end{figure}
\vspace{-1em}

In Table \ref{tab:doRes}, we observe that the Zono-$F$ method is the fastest across all models, while the adaptive methods are the slowest. Moreover, only Zono-$F$ is able to complete the reach set computation for all three models, while the other methods time out. In terms of the size of the computed reach sets, we compare the last reach set obtained in Figure \ref{fig:damOsc_compare}. In the subplots \ref{fig:osc_0z} and \ref{fig:osc_1z} we see the zoomed in reach sets, and observe that the Poly-$A$ method computes the smallest over-approximate reach set across both experiments.
Based on these results, in all subsequent experiments, we use the \emph{zono-F} method for nonlinear NODEs, the \emph{direct} method for linear NODEs, and the over-approximate star-set methods for all the $\mathscr{NN}$-Layers.

The next part of our evaluation consists of  comparing NNVODE's methods for NODEs to those of Flow* \cite{Chen2013Flow*}, Gotube \cite{gruenbacher2021gotube} and JuliaReach \cite{bogomolov2019Juliareach} across a collection of benchmarks that include a linear and nonlinear 2-dimensional spiral \cite{chen2018neural}, a Fixed-Point Attractor (FPA) \cite{musau2018arch} and a controlled cartpole \cite{gruenbacher2020lagrangian}. The computationn reachability results are displayed in Table \ref{tab:res1} with the intention to characterize major differences between tools, i.e., to show some tools are 10$\times$ to 20$\times$ faster than others for some benchmarks. It is worth noting, that we are not experts on every tool that we considered. Thus it may be possible to optimize the reachable set computation for each benchmark with depending on the tool. However, we did not do so. Instead, we attempted around 3 to 5 different parameter combinations for each benchmark, and used the best results we could obtain when comparing against the other tools. The details can be found in the Appendix.
The first three rows correspond to the linear spiral 2D model, and the subsequent three rows to the nonlinear one. The first set of observations that can be made in this context is that Flow* times out on all problems involving nonlinear neural ODEs, and that GoTube cannot obtain a solution to the reachability problem for the linear model. In terms of computation time, the results vary. Flow* is the fastest for the linear Spiral 2D model, regardless of the size of the initial set. In general, JuliaReach and NNVODE are much faster than GoTube, with JuliaReach being the fastest tool across the board. Additionally, there is not a significant difference between NNVODE and JuliaReach in the treatment of the nonlinear spiral model and FPA models. However, JuliaReach is an order of magnitude faster than NNVODE and two orders of magnitude faster than GoTube on the cartpole benchmark.
In Figure \ref{fig:res1compare} we display a subset of the reachability results from Table \ref{tab:res1} where we observe that
JuliaReach is able to compute smaller over-approximations of the reachable set on all the benchmarks except for the linear spiral model. Notably, in most cases,
GoTube computes the largest over-approximation, and this effect grows far more significantly than the other tools when the complexity of the model increases. This can be  observed in Figures \ref{fig:cartpole_a} and \ref{fig:cartpole_b}.

\begin{table}[h!]
    \caption{Results of the reachability analysis of the NODE benchmarks. All results are shown in seconds. TH stands for Time Horizon and $\delta_{\mu}$ to the input uncertainty.}
    \label{tab:res1}
    \centering
    \addtolength{\tabcolsep}{5pt}
    \begin{adjustbox}{width=0.75\columnwidth,center}
    \begin{tabular}{c|cc|cccc}

    \centering
    Name & TH ($s$) & $\delta_{\mu}$ & Flow* & GoTube & JuliaReach & NNVODE (ours) \\
    \hline
    Spiral$_{L_1}$ & 10 & 0.01 & \textbf{4.0} & -- & 10.2 & 8.0 \\
    Spiral$_{L_2}$ & 10 & 0.05 & \textbf{3.7} & -- & 7.2 & 7.4 \\
    Spiral$_{L_3}$ & 10 & 0.1 & \textbf{3.7} & -- & 7.2 & 7.3 \\ \hline
    Spiral$_{NL_1}$ & 10 & 0.01 & -- & 106.4 & 58.9 & \textbf{47.4} \\
    Spiral$_{NL_2}$ & 10 & 0.05 & -- & 106.7 & \textbf{41.7} & 46.2 \\
    Spiral$_{NL_3}$ & 10 & 0.1 & -- & 106.5 & \textbf{41.9} & 46.2 \\ \hline
    FPA$_1$ & 0.5 & 0.01 & - & 13.6 & 1.9 & \textbf{1.3}  \\
    FPA$_2$ & 2.5 & 0.01 & - & 42.4 & \textbf{5.6} & 6.6 \\
    FPA$_3$ & 10.0 & 0.01 & - & 140.2 & 28.4 & \textbf{7.5} \\ \hline
    Cartpole$_1$ & 0.1 & 1$e$-4 & - & 183.0 & \textbf{3.2} & 67.9  \\
    Cartpole$_2$ & 1.0 & 1$e$-4 & - & 1590.9 & \textbf{13.8} & 404.3  \\
    Cartpole$_3$ & 2.0 & 1$e$-4 & - & 3065.7 & \textbf{35.5} & 834.7  \\
    \hline
    \end{tabular}
    \end{adjustbox}
\end{table}
\begin{figure}[!ht]
  \centering
    \subfigure[Spiral$_{L_3}$]{

        \includegraphics[trim=98 240 120 250, clip, width=0.23\columnwidth]{figures/spiralL_compare_0.1.pdf}\label{fig:spiral_la}}
    \subfigure[Spiral$_{NL_3}$]{

        \includegraphics[trim=100 230 120 240, clip, width=0.23\columnwidth]{figures/spiral_compare_0.1_nl.pdf}
        \label{fig:spiral_nla}}
     \subfigure[FPA$_3$]{

        \includegraphics[trim=100 240 130 250, clip, width=0.23\columnwidth]{figures/dfp_compare.pdf}\label{fig:fpa_a}}
    \subfigure[Cartpole$_3$]{

        \includegraphics[trim=95 230 115 240, clip, width=0.23\columnwidth]{figures/cartpole_compare.pdf}\label{fig:cartpole_a}}
    \subfigure[Spiral$_{L_3}$]{

        \includegraphics[trim=98 230 120 240, clip, width=0.23\columnwidth]{figures/spiralL_compare_0.1_zoom.pdf}\label{fig:spiral_lc}}
    \subfigure[Spiral$_{NL_3}$]{

        \includegraphics[trim=100 230 120 240, clip, width=0.22\columnwidth]{figures/spiral_compare_0.1_nl_last.pdf}
        \label{fig:spiral_nlb}}
    \subfigure[FPA$_3$]{

        \includegraphics[trim=95 230 120 240, clip, width=0.23\columnwidth]{figures/dfp_compare_zoom.pdf}
        \label{fig:fpa_b}}
    \subfigure[Cartpole$_3$]{

        \includegraphics[trim=95 230 115 240, clip, width=0.23\columnwidth]{figures/cartpole_compare_zoom.pdf}
        \label{fig:cartpole_b}}
    \vspace{-1em}
  \caption{NODE reach set comparisons. The top 4 figures ([$a$,$d$]) show the complete reachable sets of each benchmark, while the bottom 4 correspond to the zoomed-in reach sets ($e$, $g$) and to the zoomed-in figure of the last reach set ($f$,$h$). The figures show the computed reach sets of \textbf{\textcolor{red}{GoTube}}, \textbf{NNVODE}, \textbf{\textcolor{blue}{JuliaReach}} and \textbf{\textcolor{OliveGreen}{Flow*}}.}
  \label{fig:res1compare}
\end{figure}

\subsection{Case Study: Adaptive Cruise Control (ACC)}

This case study was selected to evaluate an original NNCS benchmark used in all the AINNCS ARCH-Competitions \cite{ARCH2021ainncs,ARCH20AINNCS,manzanas2019arch_ainncs,tran2020cavtool} against NNCS with NODEs as dynamical plants learned from simulation data using 3$^{rd}$ order NODEs \cite{norcliffe2020sonode}. We demonstrate the verification of the ACC with different GNODEs learned as the plant model of the ACC and compare against the original benchmark, while using the same NN controller across all three models.
The details of the original ACC NNCS benchmark can be found in \cite{tran2020cavtool}, and the architectures of the third order neural ODEs can be found in the Appendix.

\begin{figure}[!ht]
  \centering

    \subfigure[Original]{

        \includegraphics[trim = 0 0 0 0,clip, width=0.31\columnwidth]{figures/acc/reach_orig.pdf}\label{fig:acc_orig}}
    \subfigure[Linear NODE]{

        \includegraphics[trim = 0 0 0 0,clip, width=0.31\columnwidth]{figures/acc/reach_plant.pdf}\label{fig:acc_lin}}
    \subfigure[Nonlinear NODE]{

        \includegraphics[trim = 0 0 0 0,clip, width=0.31\columnwidth]{figures/acc/reach_tanhplant.pdf}\label{fig:acc_nl}}
    \vspace{-1em}
  \caption{Adaptive Cruise Control comparison. In \textcolor{red}{red} we display the relative distance of the original plant, in \textcolor{blue}{blue} the linear NODE, in \textbf{black} the nonlinear NODE, and in \textcolor{magenta}{magenta}, the safe distance.}
  \label{fig:acc}
\end{figure}

We considered all three models using the same initial conditions and present the results in Figure \ref{fig:acc}. The reachable sets obtained from all three plants are largely similar, and we can guarantee that all the models are safe since the intersection between the safe distance and the relative distance is empty. When we consider the size of the reachable sets, one can see that the original model returns the smallest reach set, whereas the linear model boasts the largest one. In all, the biggest difference between the plant dynamics is the computation time. Here, the linear 3$^{rd}$ order NODE boasts the fastest computation time of 0.86 seconds, whereas the original plant takes 14.9 seconds, and the nonlinear 3$^{rd}$ order NODE takes 998.7 seconds.

\subsection{Classification NODEs}

Our second set of experiments considered performing a robustness analysis for a set of MNIST classification models with only fully-connected $\mathscr{NN}$-Layers and other 3 models with convolutional layers as well. There is one linear NODE in each model, and we vary the number of parameters and states across all of them to study the scalability of our methods. We evaluate the robustness of these models under an L$_{\infty}$ adversarial perturbation of $\epsilon$ = $\{$0.05, 1, 2$\}$ over all the pixels, and $\epsilon$ = $\{$2.55, 12.75, 25.5$\}$ over a subset of pixels (80). \footnote{Adversarial perturbations are applied before normalization, pixel values z$_p$ $\in$ [0,255]. } The complete evaluation of this benchmark consists of a robustness analysis using 50 randomly sampled images for each attack. We compare the number of images the neural ODEs are robust to, as well as the total computation time to run the reachability analysis for each model in Table \ref{tab:mnistRes}.

\begin{definition}[Robustness]\label{def:robustness}
    Given a classification-based GNODE $\mathcal{F}(z)$, input $z \in \mathbb{R}^j$, perturbation parameter $\epsilon \in \mathbb{R}$ and an input set $Z_p$ containing $z_p$ such that $Z_p = \{z: ||z - z_p|| \leq \epsilon\}$ that represents the set of all possible perturbations of $z$. The neural ODE is locally \textbf{robust} at $z$ if it classifies all the perturbed inputs $z_p$ to the same label as $z$, i.e., the system is \textbf{robust} if $\mathcal{F}(z_p) = \mathcal{F}(z)$ for all $z_p \in Z_p$.
\end{definition}
\begin{table}[h!]
    \caption{Robustness analysis of MNIST classification GNODEs under L$_{\infty}$ adversarial perturbations. The accuracy and robustness results are described as percentage values between 0 and 1, and the time computation corresponds to the average time to compute the reachable set per image. Columns 3-8 corresponds to $\epsilon$ = \{0.5,1,2\} over all pixels in the image, columns 9-14 corresponds to $\epsilon$ = \{2.55,12.75,25.5\} attack over a subset of pixels (80) n each image.}
    \label{tab:mnistRes}
    \centering
    \begin{adjustbox}{width=\columnwidth,center}
    \addtolength{\tabcolsep}{3pt}
    \begin{tabular}{lc|cc|cc|cc|cc|cc|cc|}

    \centering

     & & \multicolumn{6}{c|}{$\ell$$_{\infty}$} & \multicolumn{6}{c|}{$\ell$$_{\infty}$ $_{(80)}$} \\

    \cline{3-14}
     & & \multicolumn{2}{c|}{\textbf{0.5}} & \multicolumn{2}{c|}{\textbf{1}} & \multicolumn{2}{c|}{\textbf{2}} & \multicolumn{2}{c|}{\textbf{2.55}} & \multicolumn{2}{c|}{\textbf{12.75}} & \multicolumn{2}{c|}{\textbf{25.5}}\\
    \cline{3-14}
     \textbf{Name} & \textbf{Acc.} & \textbf{Rob.} & \textbf{T(s)} & \textbf{Rob.} & \textbf{T(s)} & \textbf{Rob.} & \textbf{T(s)} & \textbf{Rob.} & \textbf{T(s)} & \textbf{Rob.} & \textbf{T(s)} & \textbf{Rob.} & \textbf{T(s)} \\
    \hline
    FNODE$_{S}$ & 0.9695 & 1 & 0.0836 & 0.98 & 0.1062 & 0.98 & 0.1186 & 1 & 0.0192 & 0.98 & 0.0193 & 0.98 & 0.0295\\
    FNODE$_{M}$ & 0.9772 & 1 & 0.0948 & 1 & 0.1177 & 0.98 & 0.2111 & 1 & 0.0189 & 1 & 0.0225 & 0.98 & 0.0396 \\
    FNODE$_{L}$ & 0.9757 & 1 & 0.1078 & 1 & 0.1674 & 0.96 & 0.5049 & 1 & 0.0187 & 1 & 0.0314 & 0.96 & 0.0709 \\
    CNODE$_{S}$ & 0.9706 & 0.98 & 13.95 & 0.96 & 15.82 & 0.86 & 17.58 & 0.98 & 1.943 & 0.94 & 3.121 & 0.78 & 4.459 \\
    CNODE$_{M}$ & 0.9811 & 1 & 176.5 & 1 & 193.8 & 1 & 293.8 & 1 & 21.45 & 1 & 83.28 & 1 & 182.0 \\
    CNODE$_{L}$ & 0.9602 & 1 & 1064 & 1 & 1089 & 1 & 1736 & 1 & 234.6 & 1 & 522.7 & 1 & 779.3 \\
    \hline
    \end{tabular}
    \addtolength{\tabcolsep}{-3pt}
    \end{adjustbox}
\end{table}
\vspace{-1em}

Finally, we performed a scalability study using a set of random GNODE architectures with multiple NODEs. Here, we focus on the nonlinear methods for both the $\mathscr{NN}$-Layers and the NODEs and evaluate how our methods scale with the number of neurons, inputs, outputs and dimensions in the NODE. The main challenge of these benchmarks is the presence of multiple NODEs within the GNODE. There are a total of 6 GNODEs (XS,S,M,L,XL, and XXL), all with the same number of layers. However, we increase the number of inputs, outputs, and parameters across all the layers of the GNODEs. Here, XS corresponds to the smallest model and XXL to the largest. A description of these architectures can be found in the Appendix.

Several trends can be observed from Table \ref{tab:randomRed}. The first is that in general, smaller models have smaller reach set computation times, with two notable exceptions: the first run with XS and the last experiment with XL. Furthermore, one can observe that the largest difference in the reach set computation times comes from increasing the number of states in the NODEs, which are $\{2,3,4,4,5,5\}$ for both NODEs in every model in $\{$XS,S,M,L,XL, and XXL$\}$ respectively, while increasing the input and output dimensions ($\{1,2,2,3,3,4\}$ respectively) does not affect the reachability computation as much. This is because of the complexity of nonlinear ODE reachability as state dimensions increase, while for $\mathscr{NN}-Layers$, increasing the size of the inputs or neurons by 1 or a few units does not affect the reachability computation as much.

\vspace{-1em}
\begin{table}[h!]
    \caption{Computation time of the reachability analysis of the randomly generated GNODEs. Results are shown in seconds.}
    \label{tab:randomRed}
    \centering
    \begin{adjustbox}{width=0.7\columnwidth,center}
    \addtolength{\tabcolsep}{8pt}
    \begin{tabular}{l|cccccc}

    \centering
     & \textbf{XS} & \textbf{S} & \textbf{M} & \textbf{L} & \textbf{XL} & \textbf{XXL} \\
    \hline
    $\delta_{\mu}$ = 0.01 & 57.0 & 16.2 & 59.3 & 61.8 & 168.3 & 223.9 \\
    $\delta_{\mu}$ = 0.02 & 3.4 & 11.4 & 42.2 & 41.0 & 262.4 & 115.4 \\
    $\delta_{\mu}$ = 0.04 & 3.1 & 10.3 & 37.6 & 72.9 & 1226.3 & 243.6 \\
    \hline
    \end{tabular}
    \addtolength{\tabcolsep}{-8pt}
    \end{adjustbox}
\end{table}
\vspace{-1em}
\vspace{-1em}

\section{Related Work}

\textbf{Analysis of Neural ODEs.} To the best of our knowledge, this is the first empirical study of the formal verification of neural ODEs as presented in the general neural ODE (GNODE) class. Some other works have analyzed neural ODEs, but are limited to a more restricted class of neural ODEs with only purely continuous-time models. We refer to these models in this paper as NODEs. The most comparable work is a theoretical inquiry of the neural ODE verification problem using Stochastic Lagrangian Reachability (SLR) \cite{gruenbacher2021verification}, which was later extended and implemented in a stochastic reachability analysis tool called GoTube \cite{gruenbacher2021gotube}. The SLR method is an abstraction-based technique that is able to compute a tight over-approximation of the set of reachable states, and provide stochastic guarantees in the form of confidence intervals for the derived reachable set.
Beyond reachability analysis, there have also been several works investigating the robustness of neural ODEs. In \cite{carrara2019wifs}, the robustness of neural ODE image classifiers is empirically evaluated against residual networks, which are a standard deep learning model for image classification tasks. Their analysis demonstrates that neural ODEs are more robust to input perturbations than residual networks. In a similar work, a robustness comparison of neural ODEs and standard CNNs was performed \cite{yan2020robustness}. Their work considers two standard adversarial attacks, Gaussian noise and FGSM \cite{goodfellow2015iclr}, and their analysis illustrate that neural ODEs are more robust to adversarial perturbations.
In terms of the analysis of GNODEs, to the best of our knowledge, no comparable work as been done, although
for specific models where all the $\mathscr{NN}$-Layers of a GNODE have fully-connected layers with continuous differentiable activation functions like sigmoid or tanh, it may be possible to compare our methods to other tools (with minor modifications) like Verisig \cite{Ivanov2019verisig,ivanov2021verisig20} or JuliaReach \cite{bogomolov2019Juliareach}. However, that would restrict the more general class of neural ODEs (GNODEs) that we evaluate in this manuscript.

\vspace{-1em}
\begin{table}[h!]
    \caption{Summary of related verification tools. A \checkmark means that the tool  supports verification of this class, a $\bigcirc$ means that it may be supported it, but some minor changes may be needed, a \textbf{--} means it does not support it, and a $\odot$ means that some small changes have been made to the tool for comparison and the tools has been used to verify at least one example on this class.}
    \label{tab:comparison}

    \begin{adjustbox}{width=0.86\columnwidth,center}
    \addtolength{\tabcolsep}{8pt}
    \begin{tabular}{lccccc}
    \toprule
    \centering
    \textbf{Tool} & \textbf{ODE}\footnote{} & \textbf{NODE}\footnote{} & \textbf{NN}\footnote{} & \textbf{NNCS} & \textbf{GNODE} \\
    \hline

    CORA \cite{Althoff2015arch} & \checkmark & $\odot$ & -- & -- & -- \\

    ERAN \cite{singh2018neurips} & -- & -- & \checkmark & -- & -- \\
    Flow* \cite{Chen2013Flow*} & \checkmark & $\odot$ & -- & -- & -- \\
    GoTube \cite{gruenbacher2021gotube} & \checkmark & \checkmark & -- & -- & -- \\

    JuliaReach \cite{bogomolov2019Juliareach} & \checkmark & $\odot$ & \checkmark & \checkmark & -- \\

    Marabou \cite{katz2019marabou} & -- & -- & \checkmark & -- & -- \\
    nnenum \cite{bak2021nnenum} & -- & -- & \checkmark & -- & -- \\

    ReachNN \cite{Huang2019reachNN,fan2020reachNN*} & \checkmark & $\bigcirc$ & \checkmark & \checkmark & -- \\
    Reluplex \cite{katz2017reluplex} & -- & -- & \checkmark & -- & -- \\
    ReluVal \cite{wang2018reluval} & -- & -- & \checkmark & -- & -- \\
    Sherlock \cite{Dutta2019hscc} & \checkmark & $\bigcirc$ & \checkmark & \checkmark & -- \\

    SpaceEx \cite{Frehse2011spaceex} & \checkmark & $\bigcirc$ & -- & -- & -- \\
    Verisig \cite{Ivanov2019verisig,ivanov2021verisig20} & \checkmark & $\bigcirc$ & \checkmark & \checkmark & -- \\
    NNV \cite{tran2020cavtool} & \checkmark & $\odot$ & \checkmark & \checkmark & -- \\
    \textbf{NNVODE} (ours) & \checkmark & \checkmark & \checkmark & \checkmark & \checkmark \\
    \hline
    \end{tabular}
    \addtolength{\tabcolsep}{-8pt}
    \end{adjustbox}
    \scriptsize{$^8$ODE verification is considered to be supported for any tool than can verify at least one of linear or nonlinear continuous-time ODEs. $^{9}$NODE is a specific type of ODE, so in theory any tool that supports ODE verification may be able to support NODE. However, these tools are optimized for NODE and in practice may not be able to verify most of these models as seen on our comparison with Flow*. $^{10}$We include in this category any tool that supports one or more verification methods for fully-connected, convolutional or pooling layers.}
\end{table}
\vspace{-1em}

\textbf{Verification of Neural Networks and Dynamical systems.}
The considered GNODEs are a combination of dynamical systems equations (ODEs) modeled using neural networks, and neural networks. When these subjects are treated in isolation, one finds that there are numerous studies that consider the verification of dynamical systems, and correspondingly there are numerous works that deal with the neural network verification problem. With respect to the former, the hybrid systems community has considered the verification of dynamical systems for decades and developed tools such as SpaceEx \cite{Frehse2011spaceex}, Flow* \cite{Chen2013Flow*}, CORA \cite{Althoff2015arch}, and JuliaReach \cite{bogomolov2019Juliareach} that deal with the reachability problem for discrete-time, continuous-time or even hybrid dynamics. A more comprehensive list of tools can be found in the following paper \cite{Doyen2018}. Within the realm of neural networks, the last several years have witnessed numerous promising verification methods proposed towards reasoning about the correctness of their behavior. Some representative tools include Reluplex \cite{katz2017reluplex}, Marabou \cite{katz2019marabou}, ReluVal \cite{wang2018reluval}, NNV \cite{tran2020cavtool} and ERAN \cite{singh2018neurips}. The tools have drawn inspiration from a wide range of techniques including optimization, search, and reachability analysis \cite{Liu2021}. A discussion of these approaches, along with pedagogical implementations of existing methods, can be found in the following paper \cite{Liu2021}.
Building on the advancements of these fields, as a natural progression, frameworks that consider the reachability problem involving neural networks and dynamical systems have also emerged. Problems within the space are typically referred to as Neural Network Control Systems (NNCS), and some representative tools include NNV \cite{tran2020cavtool}, Verisig \cite{Ivanov2019verisig,ivanov2021verisig20}, and ReachNN* \cite{Huang2019reachNN,fan2020reachNN*}. These tools have demonstrated their capabilities and efficiency in several works including the verification of an adaptive cruise control (ACC) \cite{tran2020cavtool}, an automated emergency breaking system \cite{tran2019emsoft}, and autonomous racing cars \cite{ivanov2020car,ivanov2021emsoft} among others, as well as participated in the yearly challenges of NNCS verification \cite{manzanas2019arch_ainncs,ARCH20AINNCS,ARCH2021ainncs}. These studies and their respective frameworks are very closely related to the work contained herein. In a sense, they deal with a restricted GNODE architecture, since their analysis combines the basic operations of neural network and dynamical system reachability in a feedback-loop manner. However, this restricted architecture would only consist of a fully-connected neural network followed by a NODE.
In Table \ref{tab:comparison}, we present a \textbf{comparison to verification tools} that can support one or more of the following verification problems: the analysis of continuous-time models with linear, nonlinear or hybrid dynamics (ODE), neural networks (NN), neural network control systems (NNCS), neural ODEs (NODE) and GNODEs as presented in Figure \ref{fig:gnode}.

\section{Conclusion \texorpdfstring{$\&$}{} Future Work.}
We have presented a verification framework to compute the reachable sets of a general class of neural ODEs (GNODE) that no other existing methods are able to solve. We have demonstrated through a comprehensive set of experiments the capabilities of our methods on dynamical systems, control systems and classification tasks, as well as the comparison to state-of-the-art reachability analysis tools for continuous-time dynamical systems (NODE).
One of the main challenges we faced was the scalability of the nonlinear ODE reachability analysis as the dimension complexity of the models increased, as observed in the cartpole and damped oscillator examples. Possible improvements include integrating other methods into our framework such as \cite{li2020formats}, which improves the current nonlinear reachability analysis via an improved hybridization technique that reduces the sizes of the linearization domains, and therefore reduces overapproximation error.
Another approach would be to make use of the Koopman Operator linearization prior to analysis in order to compute the reach sets of the linear system, which are easier and faster to compute as observed from the experiments conducted using the two-dimensional spiral benchmark \cite{Bak2021adhs}.
In terms of models and architectures, latent neural ODEs \cite{rubanova2019latent} and some of its proposed variations such as controlled neural DEs \cite{kidger2020neuralcde} and neural rough DEs \cite{morrill2021RDEs} have demonstrated great success in the area of time-series prediction, improving the performance of NODEs, ANODEs and other deep learning models such as RNNs or LSTMs in time-series tasks. The main idea behind these models is to learn a \emph{latent} space from the input of the neural ODE from which to sample and predict future values. In the future, we will analyze these models in detail and explore the addition of verification techniques that can formally analyze their behavior.

\subsubsection*{Acknowledgements}
The material presented in this paper is based upon work supported by the National Science Foundation (NSF) through grant numbers 1910017 and 2028001, the Defense Advanced Research Projects Agency (DARPA) under contract number FA8750-18-C-0089, and the Air Force Office of Scientific Research (AFOSR) under contract number FA9550-22-1-0019. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of AFOSR, DARPA, or NSF.

\small
\bibliographystyle{splncs04}