OS-net: Orbitally Stable Neural Networks
We introduce OS-net (Orbitally Stable neural NETworks), a new family of
neural network architectures specifically designed for periodic dynamical data.
OS-net is a special case of Neural Ordinary Differential Equations (NODEs) and
takes full advantage of the adjoint method based backpropagation method.
Utilizing ODE theory, we derive conditions on the network weights to ensure
stability of the resulting dynamics. We demonstrate the efficacy of our
approach by applying OS-net to discover the dynamics underlying the R\"{o}ssler
and Sprott's systems, two dynamical systems known for their period doubling
attractors and chaotic behavior.
Marieme Ngom, Carlo Graziani
# OS-net: Orbitally Stable Neural Networks

###### Abstract

We introduce OS-net (Orbitally Stable neural NETworks), a new family of neural network architectures specifically designed for periodic dynamical data. OS-net is a special case of Neural Ordinary Differential Equations (NODEs) and takes full advantage of the adjoint method based backpropagation method. Utilizing ODE theory, we derive conditions on the network weights to ensure stability of the resulting dynamics. We demonstrate the efficacy of our approach by applying OS-net to discover the dynamics underlying the Rossler and Sprott's systems, two dynamical systems known for their period doubling attractors and chaotic behavior.

## 1 Introduction

The study of periodic orbits of systems of the form

\[\dot{\mathbf{x}}=f(\mathbf{x}),\ \mathbf{x}(0)=\mathbf{x}_{0},\ \mathbf{x}\in\mathbf{U}\subset\mathbf{R}^ {n} \tag{1.1}\]

is an important area of research within the field of nonlinear dynamics with applications in both the physical (astronomy, meteorology) and the nonphysical (economics, social psychology) sciences. In particular, periodic orbits play a significant role in chaos theory. In [6], chaotic systems are defined as systems that are sensitive to initial conditions, are topologically transitive (meaning that any region of the phase space can be reached from any other region), and have dense periodic orbits. Notably, chaotic systems are constituted of infinitely many Unstable Periodic Orbits (UPOs) which essentially form a structured framework, or a "skeleton", for chaotic attractors. A periodic orbit is (orbitally) unstable if trajectories that start near the orbit do not remain close to it. Finding and stabilizing UPOs is an interesting and relevant research field with numerous applications such as the design of lasers [18], the control of seizure activities [22] or the design of control systems for satellites [30]. An important tool when studying the stability of periodic orbits of a given system is the Poincare or return map which allows one to study the dynamics of this system in a lower dimensional subspace. It is well-known that the stability of a periodic orbit containing a point \(\mathbf{x}_{0}\) is inherently connected to the stability of \(\mathbf{x}_{0}\) as a fixed point of the corresponding Poincare map. However, explicitly computing Poincare maps has been proven to be highly challenging and inefficient [29]. With the emergence of data-driven approaches, researchers in [2] proposed a data-driven computation of Poincare maps using the SINDy method [4]. Subsequently, they leveraged this technique in to develop a method for stabilizing UPOs of chaotic systems [3].

As a matter of fact, researchers have been increasingly exploring the intersection of machine learning and differential equations in recent years. For example, Partial Differential Equations (PDEs) inspired neural network architectures have been developed for image classification in [19, 27]. On the other hand, data-driven-based PDE solvers were proposed in [24] while machine learning has been effectively utilized to discover hidden dynamics from data in [17, 4, 21]. One notable example of such intersectional work is Neural Ordinary Differential Equations (NODEs), which were introduced in [5]. NODEs are equivalent to continuous
residual networks that can be viewed as discretized ODEs [9]. This innovative approach has led to several extensions that leverage well-established ODE theory and methods [7; 34; 14; 35; 33; 9] to develop more stable, computationally efficient, and generalizable architectures.

In the present work, we aim at learning dynamics obtained from chaotic systems with a shallow network featuring a single hidden layer, wherein the network's output serves as a solution to the dynamical system

\[\dot{\mathbf{x}}=\mathbf{W}_{d}^{T}\sigma(\mathbf{W}_{e}^{T}\mathbf{x}+\mathbf{b}_{e}),\ \ \mathbf{x}(0)=\mathbf{x}_{0}. \tag{1.2}\]

where \(\mathbf{W}_{e}\) are the input-to-hidden layer weights, \(\mathbf{b}_{e}\) the corresponding bias term, \(\mathbf{W}_{d}\) the hidden-to-output layer weights, and \(\sigma\) the activation function of the hidden layer. The proposed network is a specific case of NODEs and fully utilizes the adjoint method-based [16] weight update strategy introduced in [5]. Our primary objective is to establish sufficient conditions on the network parameters to ensure that the resulting dynamics are orbitally stable. We base our argument on the finding that the stability of Poincare maps is equivalent to the stability of the first variational equation associated with the dynamical system under consideration [29]. We then build on the stability results of linear canonical systems presented in [11] to derive a new regularization strategy that depends on the matrix \(\mathbf{J}=\mathbf{W}_{e}^{T}\mathbf{W}_{d}^{T}\) and not on the weight matrices taken independently. We name the constructed network OS-net for Orbitally Stable neural NETworks.

Since we are dealing with periodic data, the choice of activation function is critical. Indeed, popular activation functions such as the sigmoid or the \(\tanh\) functions do not preserve periodicity outside the training region. A natural choice would be sinusoidal activations however these do not hold desired properties such as monotonicity. Furthermore, they perform poorly [15] on the training phase because the optimization can stagnate in a local minimum because of the oscillating nature of sinusoidal functions. In [13], the authors constructed a Fourier neural network (i.e a neural network that mimics the Fourier Decomposition) [23; 36] that uses a \(\sin\) activation but had to enforce the periodicity in the loss function to ensure that periodicity is conserved outside of the training region. The activation functions \(x+\frac{1}{a}\sin^{2}(ax)\) -called snake function with frequency \(a\)- and \(x+\sin x\) were proposed in [37] for periodic data and were proven to be particularly well suited to periodic data. As such, we use both these activation functions in this work.

This paper is organized as follows: in section (2) we present the OS-net's architecture and the accompanying new regularization strategy. In section (3) we showcase its performance on simulated data from the chaotic Rossler and Sprott systems and perform an ablation study to assess the contributions of the different parts of OS-net.

## 2 Building OS-net

### Background

In this chapter, we recall the main results on the stability of periodic orbits of dynamical systems we will be using to build OS-net. We refer readers to the appendices for more details about orbits of dynamical systems.

We consider the system

\[\dot{\mathbf{x}}=f(\mathbf{x}),\ \ \mathbf{x}(0)=\mathbf{x}_{0}. \tag{2.1}\]

and suppose it has a periodic solution \(\phi(t,x_{0})\) of period \(T\). We denote \(\gamma(x_{0})\) a periodic orbit corresponding to \(\phi(t,x_{0})\). Stability of periodic orbits have been widely studied in the literature. It is, in particular well-known [29; Chapter 12] that the stability of periodic orbits of Equation (2.1) is linked to the stability of its First Variational (FV) problem

\[\dot{\mathbf{y}}=\mathbf{A}(t)\mathbf{y},\ \ \mathbf{A}(t)=d\left(f(\mathbf{x})\right)_{(\mathbf{\Phi}(t, \mathbf{x}_{0}))}\ \ \text{and}\ \mathbf{A}(t+T)=\mathbf{A}(t). \tag{2.2}\]

which is obtained by taking the gradient of Equation (2.1) with respect to \(x\) at \(\phi(t,x_{0})\). As such, the first variational problem describes the dynamics of the state variable \(y=d\left(\phi(t,t_{0},\mathbf{x})\right)\) and is a linear system as the matrix \(\mathbf{A}(t)\) does not depend on \(y\).

To assess the stability of OS-net, we investigate the first variational equation associated with Equation (1.2). OS-net's FV is given by

\[\dot{\mathbf{y}}=\mathbf{W}_{d}^{T}diag\left(\sigma^{{}^{\prime}}\left(\mathbf{W}_{e}^{T} \phi(t,\mathbf{x}_{0})+\mathbf{b}_{e}\right)\right)\mathbf{W}_{e}^{T}\mathbf{y},\]

and if we make the change of variables \(\mathbf{z}=W_{e}^{T}\mathbf{y}\), this equation becomes

\[\dot{\mathbf{z}}=\mathbf{J}\mathbf{H}(t)\mathbf{z}, \tag{2.3}\]

where \(J=\mathbf{W}_{e}^{T}\mathbf{W}_{d}^{T}\) and \(\mathbf{H}(t)=diag\left(\sigma^{{}^{\prime}}\left(\mathbf{W}_{e}^{T}\phi(t,\mathbf{x}_{0} )+\mathbf{b}_{e}\right)\right)\) is periodic. This formulation can be seen as a generalization of linear canonical systems with periodic coefficients

\[\dot{\mathbf{y}}=\lambda\mathbf{J}_{m}\mathbf{H}(t)\mathbf{y} \tag{2.4}\]
where

\[\mathbf{J}_{m}=\begin{pmatrix}0&\mathbf{I}_{m}\\ -\mathbf{I}_{m}&0\end{pmatrix},\mathbf{I}_{m}\text{ is the identity matrix of size }m,\]

\(\mathbf{H}\) is a periodic matrix-valued function and \(\lambda\in\mathrm{R}\). Stability of such systems was extensively studied in [12] and in particular in [11]. We recall the main definitions and results from [11] and build upon these to derive stability conditions for OS-net. In particular, we give the definition of stability zones for Equation (2.4) and provide the main stability results we will base our study on.

**Definition 1**.: _A point \(\lambda=\lambda_{0}\)\((-\infty<\lambda_{0}<\infty)\) is called a \(\lambda\)-point of stability of Equation (2.4) if for \(\lambda=\lambda_{0}\) all solutions of Equation (2.4) are bounded on the entire t-axis. If, in addition, for \(\lambda=\lambda_{0}\) all solutions for any equation_

\[\dot{\mathbf{y}}=\lambda\mathbf{J}_{m}\mathbf{H}_{1}(t)\mathbf{y}\]

_with a periodic symmetric. matrix valued function \(H_{1}(t)=H_{1}(t+T)\) sufficiently close to \(H(t)\) are bounded, then \(\lambda_{0}\) is a \(\lambda\)-point of strong stability of Equation 2.4._

The set of \(\lambda\)-point of strong stability of Equation 2.4 is an open-set that decomposes into a system of disjoint open intervals called \(\lambda\)-zones of stability of Equation2.4. If a zone of stability contains the point \(\lambda=0\) then it is called a central zone of stability.

**Definition 2**.: _We say that Equation (2.4) is of positive type if_

\[\mathbf{H}\in\mathrm{P}_{n}(T)=\{\mathbf{A}(t)\text{ symmetric }s.t\text{ }\mathbf{A}(t) \geq 0\text{ }(0\leq t\leq T)\text{ and }\int_{0}^{T}\mathbf{A}(t)dt>0\}.\]

\(\mathbf{A}(t)\geq 0\) _means \(\forall x\in\mathbf{R}^{n},\langle\mathbf{A}(t)\mathbf{x},\text{ }\mathbf{x}\rangle\geq 0\) and \(\int_{0}^{T}\mathbf{A}(t)dt>0\) means \(\int_{0}^{T}\langle\mathbf{H}(t)\mathbf{x},\text{ }\mathbf{x}\rangle dt>0\)_

**Definition 3**.: _Let \(\mathbf{A}\) be a square matrix with non-negative elements. We denote by \(\mathcal{M}(\mathbf{A})\) the least positive eigenvalue among its eigenvalues of largest modulus. Note that Perron's theorem (1907) guarantees the existence of \(\mathcal{M}(\mathbf{A})\)[10]._

We can now state the main result we will derive our regularization from:

**Theorem 1** ([11] section 7, criterion \(I_{n}\)).: _A real \(\lambda\) belongs to the central zone of stability of an Equation (2.4) of positive type, if_

\[|\lambda|<2\mathcal{M}^{-1}(\mathbf{C})\]

_where \(\mathbf{C}=\mathbf{J}_{m_{n}}\int_{0}^{T}\mathbf{H}_{a}(t)\). If \(K\) is a matrix, \(K_{a}\) is the matrix obtained by replacing the elements of \(K\) by their absolute values._

The proof of this theorem is recalled in Appendix (B).

### Architecture and stability of OS-net

To base the stability of OS-net on stability theory for systems of type Equation (2.4), we need the matrix-valued function \(\mathbf{H}(t)\) and the matrix \(\mathbf{J}\) in Equation (2.3) to be respectively of positive type and skew-symmetric. To ensure \(\mathbf{H}(t)\) is of positive type, it is sufficient to use activation functions that are increasing since they have positive derivatives and diagonal matrices with positive elements are of positive type. Fortunately, many common activation functions (\(\tan\) or \(sigmoid\)) have that property. In this paper, we use the strictly monotonic activation functions \(x+\sin(x)\) and (the snake function) \(x+\frac{1}{a}\sin^{2}(ax),\ a\in\mathbf{R}\) displayed in Figure (2.1). These activation functions were proved to be able to learn and extrapolate periodic functions in [37].

Let us now pay attention to the matrix \(\mathbf{J}=\mathbf{W}_{e}^{T}\mathbf{W}_{d}^{T}\). To ensure \(J\) is skew-symmetric, we introduce the matrices \(\mathbf{W}\in\mathrm{M}_{n,2m}(\mathbf{R})\) and \(\mathbf{K}\in\mathrm{M}_{2m}(\mathbf{R})\) where \(n\) is the input size (i.e the size of \(x\)) and \(2m\) the size of the hidden layer. We then set \(\mathbf{W}_{e}^{T}=\mathbf{W}\), \(\mathbf{W}_{d}^{T}=\mathbf{\Omega}\mathbf{W}^{T}\), and \(\mathbf{\Omega}=\mathbf{K}-\mathbf{K}^{T}\). Note that the size of the hidden layer which is the size of \(\Omega\) needs to be even. Otherwise, \(\mathbf{\Omega}\) would be a singular matrix. The elements of the matrices \(\mathbf{W}\) and \(\mathbf{K}\) are the hyperparameters of the network that will be optimized during training. Now, knowing that any real skew-symmetric matrix \(\mathbf{J}\) is congruent to \(\mathbf{J}_{m}\)[32], there exists a real invertible matrix \(\mathbf{S}\) such that

\[\mathbf{J}_{m}=\mathbf{S}^{T}\mathbf{J}\mathbf{S}\]

and Equation (2.3) is equivalent to Equation (2.4). In fact, let \(\mathbf{u}=\mathbf{S}^{T}\mathbf{z}\) in Equation (2.3), we obtain

\[\dot{\mathbf{u}}=\lambda\mathbf{J}_{m}\mathbf{S}^{-1}\mathbf{H}(t)(\mathbf{S}^{-1})^{T}\mathbf{u}= \lambda\mathbf{J}_{m}\tilde{\mathbf{H}}(t)\mathbf{u}\]
where \(\tilde{\mathbf{H}}(t)=\mathbf{S}^{-1}\mathbf{H}(t)(\mathbf{S}^{-1})^{T}\in\mathrm{P}_{n}(T)\). We can now apply Theorem (1) to OS-net and state that OS-net is stable if

\[1<2\mathcal{M}^{-1}\left(\mathbf{J}_{a}\int_{0}^{T}\mathbf{H}(t)dt\right). \tag{2.5}\]

Note that since \(\mathbf{H}(t)\) is a diagonal matrix with positive elements, \(\mathbf{H}_{a}(t)=\mathbf{H}(t)\). We can now prove the following result that will justify our regularization strategy:

**Corollary 1.1**.: _Suppose the activation function \(\sigma\) is strictly increasing with a uniformly bounded derivative. Then, OS-net is stable if_

\[||\mathbf{J}_{a}||_{2}<\frac{2}{LT} \tag{2.6}\]

_where \(L\) is the superior bound of the derivative of the activation function._

Proof.: Let \(\mu\) be any eigenvalue of \(\mathbf{J}_{a}\int_{0}^{T}\mathbf{H}(t)dt\) then \(|\mu|\leq\left\|\mathbf{J}_{a}\int_{0}^{T}\mathbf{H}(t)dt\right\|_{2}\). Knowing that any norm in \(\mathbf{R}^{n,n}\) can be rescaled to be submultiplicative (i.e. \(||AB||_{2}\leq||A||_{2}||B||_{2}\)), we obtain

\[|\mu|\leq\|\mathbf{J}_{a}\|_{2}\left\|\int_{0}^{T}\mathbf{H}(t)dt\right\|_{2}\]

which leads to

\[\mathcal{M}\left(\mathbf{J}_{a}\int_{0}^{T}\mathbf{H}(t)dt\right)\leq||\mathbf{J}_{a}||_{2 }\left\|\int_{0}^{T}\mathbf{H}(t)dt\right\|_{2}\]

If \(L\) is the superior bound of the derivative of the activation function, then, since \(\mathbf{H}(t)\) is a diagonal matrix, we have \(\mathcal{M}\left(\mathbf{J}_{a}\int_{0}^{T}\mathbf{H}(t)dt\right)\leq LT||\mathbf{J}_{a}|| _{2}\). Therefore, OS-net is stable if \(1<\frac{2}{LT}||\mathbf{J}_{a}||_{2}^{-1}\) i.e \(||\mathbf{J}_{a}||_{2}<\frac{2}{LT}\). 

All in all our minimization problem becomes

\[L(\mathbf{x}_{o},g(f(\mathbf{x}_{0})))=||g(f(\mathbf{x}_{0}))-\mathbf{x}_{0}||_{2}^{2}\ \ \text{s.t.}\ \ ||\mathbf{J}_{a}||_{2}<\frac{2}{LT} \tag{2.7}\]

and this formulation is equivalent [1] to

\[L(\mathbf{x}_{0},g(f(\mathbf{x}_{0})))=||g(f(\mathbf{x}_{0}))-\mathbf{x}_{0}||_{2}^{2}+\alpha|| \mathbf{J}_{a}||_{2}^{2} \tag{2.8}\]

where \(\alpha\in\mathbf{R}\) can be fine-tuned using cross-validation. We thus have derived a new regularization strategy that stabilizes the network. By controlling the norm of \(J_{a}=|\mathbf{W}_{e}^{T}\mathbf{W}_{d}^{T}|\), we ensure solutions of Equation (2.3) and consequently periodic orbits of Equation (1.2) are stable. We validate these claims in the next section by running a battery of tests on simulated data from dynamical systems known for their chaotic behavior.

Figure 2.1: Left: Activation functions \(x+\sin(x)\) and \(x+\frac{1}{0.2}\sin^{2}(0.2x)\). Right: Derivatives of the activation functions
Numerical results

In this section, we showcace the learning capabilities and stability of OS-net on different regimes of the Rossler [20] and of the Sprott systems [25]. In all of the following experiments, the data was generated using Matlab's ode45 solver. We take snapshots at different time intervals to obtain the data used to train OS-net..

We used the \(LBFGS\) optimizer with a learning rate \(lr=1\). and the strong Wolfe [31] line search algorithm for all the experiments. Our code uses Pytorch and all the tests were performed on a single GPU1 using Argonne Leadership Computing Facility (ALCF)'s Theta/ThetaGPU [8].

Footnote 1: We base our code on the Neural ode implementation in [28]

### _The Rossler system_

As in [3], we consider the Rossler system

\[\dot{x} =-y-z \tag{3.1}\] \[\dot{y} =x+0.1y\] (3.2) \[\dot{z} =0.1+z(x-c) \tag{3.3}\]

where \(c\in\mathbf{R}\). Rossler introduced this system as an example of simple chaotic system with a single nonlinear term (\(zx\)). As \(c\) increases, this system displays period doubling bifurcations leading to chaotic behavior. Here, we consider the values \(c=6\) and \(c=18\).

#### 3.1.1 c = 6, period-2 attractor

First, we set \(c=6\) and initialize the trajectory at \([x_{0},y_{0},z_{0}]=[0,-9.1238,0]\). In this regime, the Rossler system possesses a period-2 attractor [3].

We generate the training data by solving the Rossler system using Matlab's ode45 solver with a time step of \(0.001\) from \(t=0\) to \(t=10\). We then take snapshots of the simulated data every \(50\) step and feed it to OS-net. We build OS-net using the Runge-Kutta 4 (RK4) algorithm with a time step of \(0.005\). We chose the snake activation function \(x+\frac{1}{0.2}\sin^{2}(0.2x)\) and set the number of nodes in the hidden layer to be \(2\times 16\). We set \(\alpha=0.07\) in Equation (2.8) and use \(10\) epochs.

Figure (3.1) (left) shows the training output for the \(y\) component. OS-net was able to learn the dynamics accurately by the end of training. The norm of \(J_{a}\) is approximately \(0.99\) after training as recorded in Table (3.1). In this case, Inequality (1.1) is not strictly enforced but the norm of the matrix \(J_{a}\) is controlled enough so

Figure 3.1: Left: Training data in black along with the learned dynamics in red in the training time interval \([0,10]\). Right: Target dynamics in gray along with the data generated by OS-net on the time interval \([0,100]\).
that OS-net renders stable orbits. The elements of the matrix \(\Omega\) are concentrated in \([-0.7,\ 0.7]\) as shown in Figure 3.3.

We validate OS-net by propagating a trajectory initialized \([x_{0},y_{0},z_{0}]=[0,-9.1238+0.01,0]\) using the learned dynamics. Figure (3.1) (right) shows prediction using OS-net up to \(t=100\) and displays the accuracy of this prediction when compared to the correct dynamics. We assess the stability of OS-net by propagating the trajectory to \(t=10000\). OS-net converges to a stable period-1 attractor while the Rossler system converges to a period-2 one as show

\begin{table}
\begin{tabular}{l l l} \hline \multicolumn{3}{c}{Part} \\ \hline System & Attractor type & \(||J_{a}||\) \\ \hline Rössler, \(c=6\) & Period-2 & 0.9937 \\ Rössler, \(c=18\) & Chaotic & 0.6318 \\ Sprott, \(\mu=2.1\) & Period-2 & 0.0085 \\ \hline \end{tabular}
\end{table}
Table 3.1: Norm of \(J_{a}\)

Figure 3.2: Left: Rössler’s period-2 attractor. Right: stable OS-net period-1 attractor (right)
Ablation studyWe compare OS-net with a network obtained by keeping the same architecture and settings as in Section (3.1.1) but with the regularization in Equation (2.8) switched off. The left side of Figure (3.4) shows that training was successful while the right side shows the dynamics learned by the unregularized network diverge from the true dynamics in the time interval \([0,10]\). This shows the role of the regularization term in stabilizing the dynamics learned by OS-net.

#### 3.1.2 c = 18, chaotic behavior

We now set \(c=18\) and initialize the trajectory at \([x_{0},y_{0},z_{0}]=[0,-22.9049,0]\). The Rossler system displays a chaotic behavior in this regime. We generate the training data as before but take snapshots every \(10\) steps. For OS-net, we use RK4 with a step size of \(0.005\) and \(x+\sin(x)\) as an activation function. The hidden layer size is \(2\times 32\) and the penalty coefficient \(\alpha=2\).

Figure (3.5) shows the training output and confirms the ability of OS-net to learn the target dynamics. We then use the learned dynamics to generate a trajectory starting at \([x_{0},y_{0},z_{0}]=[0,-22.9049+0.01,0]\). Since we are dealing with a chaotic system, the learned dynamics should not be expected to reproduce the training data [2]. Figure (3.5) shows that OS-net was able to track the chaotic system up to \(t\approx 30\). The norm of the matrix \(J_{a}\) was approximately \(0.6318\) at the end of training as recorded in Table (3.1). Furthermore, the elements of the matrix \(\Omega\) were concatenated between \(-0.5\) and \(+0.5\). Figure (3.6) displays the chaotic Rossler system and the stable attractor obtained by propagating OS-net's learned dynamics from \(t=0\) to \(t=10000\).

### Simplest quadratic (Sprott's) chaotic flow

We consider the following system

\[\dot{x} =y \tag{3.4}\] \[\dot{y} =z\] \[\dot{z} =-\nu z-x+y^{2}\]

where \(\nu\in\mathbf{R}\). This system was introduced in [25] and also has period doubling bifurcations as \(\nu\) varies. Here we set \(\nu=2.1\) which yields a period-2 attractor for Equation (3.4).

We initialize the trajectory at \([x_{0},y_{0},z_{0}]=[5.7043,0.0,-2.12778]\) and solve the system using ode45 on the time interval \([0,15]\) with a step size of \(0.001\). We then take snapshots every \(10\) step and use the data for

Figure 3.4: Left: Training data in black along with the learned dynamics in red in the training time interval \([0,10]\). Right: Target dynamics in gray along with the data generated by OS-net on the time interval \([0,100]\).
training. OS-net is solved using RK4 with a step size of \(0.01\) and \(x+\frac{1}{0.3}\sin^{2}(0.3x)\) as an activation function. The hidden layer has \(2\times 16\) nodes and the penalty coefficient \(\alpha=1\).

We show in Figure (3.7) (left) the dynamics learned by OS-net for the \(y\) component after \(20\) epochs. Figure (3.7) (right) also shows how well OS-net tracks the original system in the interval \(t=0\) to \(t=100\). In this case, the norm of the matrix \(J_{a}\) was approximately \(8e-3\) and the elements of the matrix \(\Omega\) are in the interval \([-0.7,\ 0.7]\). Inequality (1.1) is strictly enforced here. We then assess the stability of the learned dynamics by generating a trajectory starting at \([x_{0},y_{0},z_{0}]=[5.7043+0.01,0.0,-2.12778]\) and evolving it from \(t=0\) to \(t=10000\). Figure (3.8) shows the period-2 attractor of the original system and the stable period-1 OS-net orbit.

NoteThe current implementation of OS-net uses the adjoint method presented in [5] which accumulates numerical errors when integrating backward. We circumvent that by using RK4 with a small step size. This

Figure 3.5: Left: Training data in black along with the learned dynamics in red in the training time interval \([0,10]\). Right: Target dynamics in gray along with the data generated by OS-net on the time interval \([0,100]\).

Figure 3.6: Left: Chaotic Rössler attractor. Right: Stable period-1 OS-net attractor
results in a computationally expensive implementation that can be improved using the methods proposed in [14, 35, 34] that we plan on incorporating into OS-net in the future.

## 4 Conclusion

We have presented a new family of stable neural network architectures for periodic dynamical data. The proposed architecture is a particular case of NODES with dynamics represented by a shallow neural network. We leveraged well-grounded ode theory to propose a new regularization scheme that controls the norm of the product of the weight matrices of the network. We have validated our theory by learning the Rossler and Sprott's systems in different regimes including a chaotic one. In all the regimes considered, OS-net was able to track the exact dynamics and converge to a stable period-1 attractor. That indicates that OS-net is a promising network architecture that can handle highly complex dynamical systems. In the future, we aim at controlling the parameters of the systems of interest by incorporating them into the state vectors that OS-net

Figure 3.7: Left: Training data in black along with the learned dynamics in red in the training time interval \([0,15]\). Right: Target dynamics in gray along with the data generated by OS-net on the time interval \([0,100]\).

Figure 3.8: Left: Sprott’s period-2 attractor. Right: stable period-1 OS-net attractor
aims at learning. Additionally, we plan on using OS-net to learn and monitor the orbits of celestial objects that have short orbital periods such as certain exoplanets or three-body systems like Mars-Phobos. This extension of OS-net's applications holds great potential in providing a broader range of stable periodic orbits for the design of spatial missions.

## Acknowledgments

This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program and the FASTMath Institute under Contract No. DE-AC02-06CH11357 at Argonne National Laboratory. Government License. The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ("Argonne"). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. [http://energy.gov/downloads/doe-public-access-plan](http://energy.gov/downloads/doe-public-access-plan).
