Prediction and control of spatiotemporal chaos by learning conjugate tubular neighborhoods

I present a data-driven predictive modeling tool that is applicable to high-dimensional chaotic systems with unstable periodic orbits. The basic idea is using deep neural networks to learn coordinate transformations between the trajectories in the periodic orbits' neighborhoods and those of low-dimensional linear systems in a latent space. I argue that the resulting models are partially interpretable since their latent-space dynamics is fully understood. To illustrate the method, I apply it to the numerical solutions of the Kuramoto--Sivashinsky partial differential equation in one dimension. Besides the forward-time predictions, I also show that these models can be leveraged for control.


Introduction
The trade-off between interpretability and accuracy is encountered in many facets of data science [1,2].Typically, the more complicated a model is, the higher is its prediction accuracy; and conversely, human-interpretable models perform poorly in complex tasks.In the data-driven studies of dynamical systems, this conundrum presents itself as unexplainable (black-box) models of chaos, especially in systems with many degrees of freedom, e.g.spatiotemporal phenomena [3,4].While such black-box models could be immensely useful for making predictions, especially in systems without known models, their contribution to the basic understanding of a system is little.
While several interpretable data-driven dynamical modeling tools exist in the literature, their applications have been predominantly limited to low-dimensional chaos or highdimensional non-chaotic systems.Arguably the most popular interpretable dynamical system modeling method is the dynamic mode decomposition (DMD) [5][6][7], with which one seeks a best-fit linear system to a given dataset and perform a modal expansion for predictions.Although a rigorous connection between DMD modes and the eigenfunctions of the Koopman operator exists [8], in practice the conditions for this correspondence are rarely satisfied except in the close neighborhoods of fixed points [9] and periodic solutions [8].Alternative approaches [10,11] to the Koopman formalism utilized deep learning for identifying the Koopman eigenfunctions, which evolve linearly.While the applications of Ref. [10] included non-chaotic nonlinear systems such as nonlinear pendulum and wake flow behind a cylindrical obstacle, Ref. [11] demonstrated that such models can also predict the evolution of the spatiotemporally chaotic Kuramoto-Sivashinsky system for short times.Another recent study [12] showed that one can develop reduced-order models of non-linearizable dynamics in terms of spectral submanifolds, which are smooth extensions of the eigenspaces of invariant solutions, e.g.equilibria and (quasi-)periodic orbits.Once again, the applications that were considered in [12] were systems that did not exhibit chaos.For low-dimensional chaotic systems, sparse identification methods [13] were shown to be able to identify the governing equations as well as the intrinsic coordinates on which the dynamics is sparse [14].These methods, however, rely on a set of candidate terms which increase in number very quickly with system dimension, thus, are not appropriate for systems with high-dimensional attractors.
One of the defining features of chaotic systems is the presence of a dense set of unstable periodic orbits [15].A central motivation for the present work stems from the numerical discoveries [16][17][18][19][20][21][22][23] of unstable periodic orbits in two-and three-dimensional Navier-Stokes simulations.These papers illustrated various similarities between the weakly turbulent flows and the periodic orbits and suggested a Markovian description of turbulence as transitions among the neighborhoods of these solutions in the state space.In addition, several [17,19,20,23] of these papers also reported that the periodic orbits had varying number of unstable dimensions, which can be taken as a numerical evidence for the lack of hyperbolicity [24].An important consequence of the absence of hyperbolicity is the breakdown of shadowing [25], which implies that the long trajectories arising in Navier-Stokes simulations might not correspond to "true" trajectories of the system as they explore the state space regions with different number of unstable directions.In light of the above-cited numerical results and their theoretical implications, I believe that reduced-order modeling of high-dimensional chaotic systems such as turbulent flows should not aim to model the system as a whole but rather attempt to produce an ensemble of models each applying to distinct regions of the state space.
Reduced-order data driven models of high-dimensional systems that are local in the state space can be obtained by performing clustering prior to modelling, see [26,27] for examples.Both of these studies k-means clustering [28], which groups states of the system according to their distances (measured in some norm in the state space, typically L 2 ) to one another.This approach to the state space partitioning is, however, blind to the dynamics since two nearby trajectories of a nonlinear system can evolve towards very different future states if, for example, they lie at two sides of a basin boundary.This can be readily seen in Ref. [26]'s application to the Lorenz system where the trajectories that evolve towards different "ears" of the attractor are grouped into the same cluster.In contrast, the standard binary symbolic partition of the Lorenz attractor from which one can accurately estimate fractal properties of the attractor using the periodic orbit theory [29] encodes the trajectories according to the order in which they visit each ear, see e.g.[30,31].
The purpose of the present paper is to demonstrate through an application to the onedimensional Kuramoto-Sivashinsky equation that one can produce simple data-driven models that are capable of predicting spatiotemporally chaotic dynamics by anchoring the models at the periodic orbits.Similar to [10] and [11], I utilize deep autoencoders to find transforma-tions between the system's states and a low-dimensional latent space in which the dynamics is linear.The latent space dynamics is determined to be conjugate to that of the periodic orbit's leading Floquet eigenspace.In this sense, the method is conceptually similar to the spectral subspaces of Ref. [12], although here I only aim to capture linearizable dynamics.The paper is organized as follows.Section 2 presents a brief description of the Kuramoto-Sivashinsky equation and the computational tools and introduces the periodic orbits.The modeling method and its training is described in Section 3 followed by the main results in Section 4. Results are discussed in Section 5 and the paper concludes in Section 6.In one space dimension, the Kuramoto-Sivashinksy [32,33] equation reads where the subscripts t and x denote the partial derivatives with respect to time and space, respectively; and the real-valued scalar field u(x, t) satisfies the periodic boundary condition u(x, t) = u(x + L, t).The domain length L is the sole control parameter of the system and when it is sufficiently large, the Kuramoto-Sivashinsky system exhibits chaos [34,35].(1) can be simulated by computing the truncated discrete Fourier expansion u(x) where q k = 2πk/L, and integrating the set of ordinary differential equations (ODEs) which is obtained by substituting u(x) with its Fourier series in (1).Noting that k = 0 mode is decoupled from the rest in (2) (Galilean invariance), one can set it to 0 without loss of generality to obtain a 2K-dimensional dynamical system defined by the nonlinear ODEs (2).Since v k = v * −k due to real-valuedness of u(x, t), the number of independent degrees of freedom of this system is also d = 2K.For the numerical integration of (2), I utilized the general purpose integrator odeint of scipy.integrate[36] evaluating the nonlinear term pseudospectrally [37] , where F denotes the discrete Fourier transformation.I chose L = 22.0 and truncated the Fourier expansion at K = 15 following Cvitanović et al. [38], who reported chaotic dynamics and (relative) periodic orbits at these parameters.In all results presented in the following, the simulated trajectories are sampled at time steps of δt = 0.01.As an example, Figure 1 shows a spacetime plot of u(x, t) simulated starting from a random initial condition on the system's attractor.
When u(x, t) describes a velocity field in one dimension, the L 2 inner product is interpreted as the kinetic energy density and its rate of change is given by [38] KE = P − ε , where One can interpret the observables P and ε as the instantaneous rates of power input and dissipation, respectively.

Symmetry reduction
The Kuramoto-Sivashsinky equation (1) in the periodic domain preserve its form (equivariance) under the translations u(x, t) → u(x − δx, t) and the reflection u(x, t) → −u(−x, t).Consequently, its chaotic solutions come in families of symmetry copies that can be generated by these transformations.Additionally, the presence of continuous symmetries give rise to relative periodic orbits, satisfying where δx p ∈ [0, L), n p ∈ {0, 1}, T δx and M are the translation and reflection (mirror) operators, respectively, and Φ t denotes the finite-time flow implied by the simulation of (1), i.e. u(t) = Φ t (u(0)).In the state space, the relative periodic orbits with translations correspond to two-tori (quasiperidoicity) parametrized by spatial and temporal shifts.If δx p = 0 and n p = 1, then the orbit is said to be pre-periodic because u p repeats itself after two periods since M 2 = I, where I is the identity operator.Note that the orbits with δx p ̸ = 0 and n p = 1 can be transformed to preperiodic orbits by a translation by −δx p /2 since T δx M = MT −δx .For the model reduction method presented in the following, these symmetry multiplicities are undesirable or in some cases prohibitive.Thus, one has to eliminate them by a symmetry-reduction prior to the dynamical modeling.Symmetry reduction of the Kuramoto-Sivashinsky equation by means of first Fourier mode slice and invariant polynomials was formulated in [39,40].Already in [39] it was shown that eliminating the translation degree of freedom by fixing the phase of the first Fourier mode to a set value results in fast oscillations.Although this can be remedied by a time-rescaling transformation, it may not be possible in all cases, especially if the data is collected experimentally.Therefore, here I opt for a two-step transformation starting with one that fixes the phase of the second Fourier mode to a constant value, followed by construction of complex polynomials that eliminate the remaining discrete symmetries.Leaving the details to the Appendix A, hereafter, I use ξ and S(u) to denote the symmetry-reduced states and the symmetry-reducing transformation satisfying where γ ∈ Γ = {T δx , T δx M}, such that an inverse transformation S −1 (ξ) = u ′ = γ ′ u, where ξ ∈ Γ can be found.As a result, no information other than the symmetry multiplicity is lost by the symmetry reduction since it is revertible up to a choice of a symmetry copy.

The attractor and the periodic orbits
Figure 2A illustrates the chaotic attractor of the Kuramoto-Sivashinsky system as the projection of a symmetry-reduced trajectory integrated for a time interval t ∈ [0, 10 4 ] onto the first three modes obtained from the principal component analysis (PCA) [41] of the trajectory.
In addition to the attractor, Figure 2A also shows two periodic orbits PO 1 and PO 2 as blue and orange closed curves, respectively.Figure 2B and C show the spacetime visualizations of PO 1 and PO 2 , respectively, where each orbit is plotted for t ∈ [0, 50] so that adjacent panels have the same timescale and the periods T 1 = 32.80 and T 2 = 43.61 of the orbits are indicated by the horizontal dashed lines.While PO 1 shifts by δx 1 = 10.96 after one period, PO 2 is pre-periodic, namely its final state is the reflection of its initial one.These orbits are the two that I found most-frequently by initiating Newton searches from close recurrences of the symmetry reduced chaotic trajectories, and they form the bases of the reduced-order models described in the following.
The linear stability of a relative periodic orbit is determined by the Jacobian of the ( 5)'s right-hand side given by  whose eigenvalues Λ i and eigenvectors V i are called the Floquet multipliers and Floquet vectors, respectively [42].Corresponding exponents λ i that satisfy Λ i = e λ i Tp are known as the Floquet exponents and those with negative and positive real parts correspond to the stable and unstable directions, respectively.Each relative periodic orbit of the Kuramoto-Sivashinsky system has at least two marginal eigenvalues with Re λ i = 0 corresponding to directions corresponding to temporal and spatial translations.Table 1 lists the leading (ordered in descending |Λ i |) two non-marginal Floquet multipliers and exponents of the periodic orbits.These eigenvalues were obtained by approximating the Jacobian dΦ t (u)/du along the periodic orbit in the Fourier space by integrating the gradient of (2) starting from the identity matrix as the initial condition.
3 Conjugate tubular neighborhoods 3.1 Basic assumptions and the model architecture I would like to begin with an overview of the modeling approach.As Figure 2A illustrates, the unstable periodic orbits of the Kuramoto-Sivashinsky system appear to be embedded in its chaotic attractor.Remembering also that these orbits are identified by Newton searches initiated from the states on the attractor, it is sensible to assume that these orbits can form the basis of models for approximating the chaotic trajectories in their vicinity.Guided by these observations, I construct two models each aiming to predict the dynamics in the neighborhood of one of the periodic orbits.The series of operations for obtaining predictions from each model is summarized in Figure 3, where an initial state u(t) is first symmetry reduced and then encoded into a three-dimensional latent state η.The latent state is evolved linearly to η ′ (t + τ ), and decoded back to a symmetry-reduced state ξ ′ (t + τ ) which is finally transformed back to u ′ (t + τ ) by the inverse symmetry reduction and state alignment.Here, ′ signifies that the quantity is a prediction.The transformations to and from the latent space are achieved by an autoencoder [43] which is trained indivdiually for each periodic orbit.The linear dynamics in the latent space is prescribed such that the periodic orbit is mapped to the unit circle and the transverse dynamics is conjugate [15] to that of the leading Floquet eigenspace of the periodic orbit associated with the exponents with the largest real parts.In doing so, the model aims to capture the dynamics associated with the exponentially dominant part of the periodic orbits' tangent space.
The symmetry reduction S and its inverse S −1 is described in the Section 2 and Appendix A. The symmetry-aligning transformation A(ξ ′ (t+τ ); u(t)) shown in Figure 3 is an abstraction for the following series of operations.Given a symmetry-reduced prediction ξ ′ (t + τ ) for τ ∈ {0, δt, 2δt, . ..}, first the inverse symmetry reduction S −1 (ξ ′ ) = ũ′ is obtained.Next, for each time step with τ > 0, one out of four discrete symmetry copies that can be reached by σ ∈ Σ = {I, M, T (L/2), MT } is selected, such that the resulting trajectory is the smoothest one.This is implemented by approximating the partial derivative ũ′ t via finite differences (second-order central differences for intermediate points and first-order finite-differences at the beginning and the end of the interval) and transforming ũ′ (τ The resulting trajectory is one where each state has a fixed second Fourier mode phase, which can be interpreted as a slice [39] and, thus, the remaining translations can be determined by the reconstruction equation [44].Finally, γ * = arg min γ∈Γ ∥u(t) − γu ′ (t)∥ is found and applied to the reconstructed trajectory, so that the initial states are aligned.Note that even though the inverse symmetry reduction S −1 is exact, the selection of the discrete symmetry copy relies on the estimation of ũ′ t , which is a source of numerical errors.It is, therefore, important that the symmetryreducing transformation does not result in a loss of temporal resolution.
In ref. [45], which was influential for the present one, Bramburger et al. showed that deep neural networks can be utilized for discovering conjugacies between discrete-time dynamical systems, i.e. mappings.Although in principle the approach presented in [45] can be extended to continuous-time systems by means of Poincaré sections, in practice finding suitable Poincaré sections become very challenging in high-dimensional systems because the transversality between the dynamics and a Poincaré section hyperplane is not guaranteed beyond a close neighborhood of a periodic orbit.For this reason, in this work I formulate a conjugacy-based predictive modeling tool directly for the continuous-time system.The encoder-decoder pair (E, D) constitutes an autoencoder [43] which I approximate using neural networks from 30-dimensional symmetry-reduced state space to 3-dimensional latent space and backwards.Before describing the training through which the network parameters are determined, I must explain the linear time evolution in the latent space, which corresponds to the innermost block in Figure 3. Let η = (η 1 , η 2 , η 3 ) and R i (θ) be the 3 × 3 rotation matrix, action of which rotates η about η i by θ.If the leading Floquet exponents of the periodic orbit with the period T p are a complex conjugate pair λ 1,2 = µ ± iω, the corresponding latent-space evolution is given by where θ 0 = atan2(η 2 , η 1 ) and η1 = (1, 0, 0).The second case is the purely real Floquet exponents µ 1,2 with the latent-space dynamics given by where η3 = (0, 0, 1).Note that the term R 3 (−θ 0 )η − η1 that appears in both ( 9) and ( 8) measures η's deviation from the unit circle and if it is equal to 0, then the dynamics is along the unit circle on the η 1 η 2 -plane.Figure 4 illustrates the trajectories of ( 8) and ( 9) for some choice of parameters.Each curve shown in Figure 4 corresponds to a T p = 2π-long trajectory segment and the eigenvalues are chosen such that the qualitative differences between two cases are easy to see.
To model the dynamics in the neighborhoods of PO 1 and PO 2 I choose the latent-space dynamics given by ( 8) and ( 9), respectively, with Floquet exponents equal to those of the periodic orbits shown in Table 1.As a result, the leading Floquet eigenspace of the periodic orbits and the linear systems in ( 8) and ( 9) are conjugates of one another and the task of the autoencoder is to perform the coordinate transformations between these conjugate tubular neighborhoods.

Training
Both E and D are multilayer perceptrons with two hidden layers, each of which consist of 128-nodes with SiLU activation functions [46].Together, the encoder-decoder pair have 74785 adjustable parameters.To produce the training data, the trajectories initiated from (n) , where n = 1, 2, . . ., 1000 indexes the training data, u (n) p is a random state on the periodic orbit, and δu (n) is a random perturbation with the amplitude ∥δu (n) ∥ = 10 −3 ∥u (n) p ∥.These initial conditions were integrated for c/µ 1 + T p , where µ 1 is the real part of the leading Floquet exponent of the periodic orbit and T p is its period.Only the final T p -long segment of these trajectories are saved for training after symmetry reduction and the initial c/µ 1 -long parts are discarded as transients.After experimenting with different values, I took c = 3.6 for both periodic orbits.The training trajectories are visualized along with the periodic orbits in Figure 5 as three-dimensional projections onto the principal components of the periodic orbits.
After generating the data, each model is trained via Adam algorithm [47] to minimize a loss function.Let ξ p [k], where k ∈ {0, 1, 2, . . ., ⌊T p /dt⌋} be states on the periodic orbit, η p [k] be the corresponding latent states on the unit circle and (ξ l,i , ξ l,f ) be the l-th initial and final state pair sampled randomly from the n-th training trajectory as ξ l,i = ξ (n) [0] and ξ l,f = ξ (n) [m], where n ∈ {1, 2, . . ., 1000} and m ∈ {0, 1, 2, . . ., ⌊T p /dt⌋}.The loss function is where, N b is the batch size.As in Figure 3, ′ indicates the predicted quantity, i.e. ξ ′ l,f = D(L mδt E(ξ l,i )).Note that only the very last term in (10) corresponds to the prediction errors, however, when I trained models with this error term only, I found the transformations E and D to converge to a random mapping between the training trajectories and the linear systems defined in ( 8) and ( 9).To remedy this, I included the first term of the second sum in (10), which is the standard autoencoder loss, penalizing the deviation of D from the inverse of E. Finally, to ensure that the periodic orbits are indeed mapped to the unit circle in the latent space, I needed to include the first sum in (10), where the first term penalizes the deviations of the periodic orbit from the unit circle in the latent space and the second one is the autoencoder loss for the periodic orbit.I initiate the training only using the periodic orbit with a learning rate of 10 −4 for 100 epochs.Once this initial training was complete, I provided the training data in batches of N b = 100 and ran another 100-epoch training with the learning rate 10 −4 and finally, reduced the learning rate to 10 −5 and repeated.This very last training step was terminated before 100 epochs due to the validation losses not improving more than 10 −6 in two consecutive steps.Both models, results from which are presented in the next section, are trained by this procedure and the pytorch [48] implementation using lightning [49] framework can be found in the code repository [50] accompanying this paper.

Test trajectories
As the first set of results, Figure 6 shows two test trajectories, that were initiated in the vicinity of PO 1 by the same procedure as the training trajectories but were not contained in the training dataset along with the corresponding model predictions.The top panels (A-E) of Figure 6 show a trajectory that is initially very close to the orbit with the minimum relative distance approximately d 0 ≈ 0.017 and the bottom panels (F-J) illustrate another case where the initial state is further away with d 0 ≈ 0.29.These relative distances were determined as where ξ p (t) are the symmetry-reduced states on the periodic orbit and ξ 0 is the initial symmetry-reduced state.The leftmost panels (Figure 6A and F) show the test trajectories as projections onto the PO 1 's principal components, next to their projections in the threedimensional latent space (Figure 6B and G). Figure 6C and H show the same trajectories on    the plane of observables where the horizontal and vertical axes correspond to the instantaneous rates of power input P and the time derivative of the kinetic energy KE = P − ε, respectively.Here, I chose the (P, KE)-plane as opposed to the (P, ε)-plane which is more frequently encountered in the literature because the latter results in trajectories squashed around the diagonal leaving most of the figure panel as a white space; see [38] for the comparisons of two cases.Also plotted in each trajectory visualization Figure 6A-C and F-H are the corresponding model predictions (green) and the periodic orbit (dashed blue).The rightmost space-time plots in Figure 6 show the symmetry-reconstructed test trajectories and predictions, which are visually indistinguishable from one another.The analogous test trajectories and the corresponding model predictions for PO 2 are shown in Figure 7, where the panels are organized in the same order as Figure 6.These results clearly show that the models can generalize to predict trajectories that are not included in the training dataset.

Chaotic trajectories
To see whether the trained models can predict the chaotic time evolution, I tested them using states on the trajectory depicted in Figure 1, treating them as a series of measurements, from which the models can generate predictions.Because each model is intended for the neighborhood of one of the two periodic orbits, the prediction algorithm requires a way of determining whether a prediction can be made and, if yes, which model should be used.This is achieved according to the following rules.Given a state ξ(t), the algorithm computes the autoencoder error for each model and finds the one with the minimum ϵ AE .If this error is less than the threshold ϵ th = 10%, the model is used to predict the evolution from t to t + T p , where T p is the period of the associated periodic orbit, and the algorithm repeats from ξ(t + T p + δt).If the error is larger than the threshold, then the algorithm does not attempt a prediction and repeats from the next state ξ(t + δt). Figure 8 shows the relative errors of the predictions of the fields u(x, t) and their symmetry-reduced counterparts (top) as well as the relative errors of the observables P and ε (bottom) as a function of time.The episodes for which the 10% threshold is not satisfied, thus no prediction is made, are indicated by the gray shaded regions in Figure 8.Of course, the total duration of these episodes can be reduced by increasing the 10% threshold at the expense of prediction accuracy, or conversely, the accuracy of the predictions can be improved by reducing the threshold.After some experimentation, I chose 10% as it results in reasonably accurate predictions as well as several inaccurate ones which I believe are important to discuss here as examples of what can go wrong.Figure 9 illustrates two consecutive segments (top and bottom) of the chaotic trajectory along with the corresponding predictions from the models of PO 1 's neighborhood.Beginning and ending of these episodes are indicated by the vertical dashed lines in Figure 8, where the initial time of the second episode (Figure 9F-J) overlaps with the ending of the first one (Figure 9A-E).Similar to the Figure 6 and Figure 7, Figure 9 also shows different visualizations of the trajectories and predictions as projections onto the periodic orbits' principal components (A and F), in the latent space (B and G), on the (P, KE)-plane (C and H), and as    space-time plots (D, E, I, and J).While the spacetime plots in Figure 9D and E are visually indistinguishable, Figure 9I and J shows that soon after t = 250, the original and predicted trajectory differ significantly.Noting in Figure 8 that for this episode the relative errors in the symmetry-reduced coordinates are still on the order of 0.1 whereas they jump to 1 only for the symmetry-reconstructed predictions, one can conclude that the error is due to the symmetry-aligning transformation.Naturally, the observables P and KE are not affected by this error since they are integrals over the domain (4), thus invariant under the symmetries.While the trajectories in Figure 9H show instantaneous differences, the relative errors of the observables throughout this prediction window is on the order of 0.1 as shown in Figure 8. Two trajectories in the neighborhood of PO 2 along with the corresponding model predictions are shown in Figure 10 as, once again, projections onto the principal components of the periodic orbit (A, F), in the latent space (B, G), on the plane of observables (P, KE) (C, H) and as space-time plots (D, E, I, and J).The dotted vertical lines in Figure 8 indicate the beginning and endings of two episodes that are visualized in Figure 10.Notice that towards the end of the first episode around t ≈ 330, errors computed in both symmetry-reduced and the reconstructed representations shoot above 1.0.In this case, the model diverges to a region that is outside the attractor with observable values significantly different from those observed.This is also visible towards the end of the space-time plot of Figure 10E where the colors saturate.Similarly, in the second case shown in Figure 10F-J, the prediction errors    also significantly increase towards the end of the window as seen near the right-most dotted vertical line in Figure 8.

Control
To demonstrate how the understanding the latent-space dynamics of the presented models can be leveraged for applications, I would like to present a simple control method that stabilizes PO 1 .Let F (u) denote the RHS of the Kuramoto-Sivashsinky equation (1), i.e. u t = F (u), and consider the controlled system where In the equation above, h ∈ R >0 is a constant "gain", γ * = arg min γ∈Γ ∥u − γS −1 (S(u))∥ is the state-aligning symmetry operator, and η = E(S(u)) is the encoded state.Notice that the term in the innermost parentheses in (14) penalizes the latent states' norm's difference from 1 and becomes 0, if ∥η∥ = 1.  Figure 11A and B show the space-time plot of the u(x, t) obtained from the simulation of (13) with h = 0.4 and ϵ th = 0.1 starting from the same initial condition as Figure 1.As seen, the dynamics has now become relative periodic as opposed to the chaotic shown in Figure 1.For further illustration, Figure 11C shows the symmetry-reduced controlled trajectory (green) along with the chaotic trajectory segment (orange) for t ∈ [0, 200] as a projection onto the principal components of the PO 1 , which itself is shown as the dashed blue curve.Notice that the stabilized green trajectory and the dashed blue one mostly overlap in Figure 11C.Finally, the time periodicity can also be seen in the time-series of the observables ε and P , which are plotted green in Figure 11D and E, respectively, along with the corresponding time-series obtained from the chaotic trajectory in orange.

Discussion
With the results shown in Figure 9 and Figure 10, I attempted to demonstrate that the chaotic trajectories of the Kuramoto-Sivashinsky system can be predicted using models with linear latent-space dynamics.Among the four cases presented, three of them have visible discrepancies towards the end of the prediction window, which I would like to comment on and speculate possible remedies.
First, I would like to focus on Figure 10D and E where at the very end of the prediction window the original and the predicted trajectories differ significantly.This is due to the rapid growth of the predictions' amplitude, which is also manifested in Figure 10C by the final segment's exit from the area shown.Presumably, the very last segment of the corresponding latent-space trajectory lays in a region which was not included in the training data, which in turn results in its decoding into a trajectory outside the attractor.A similar divergence also takes place at the very end of the second case shown in the bottom row of Figure 10.In practice, one could discard such unphysical predictions by comparing their amplitude to the attractor limits.Another way of eliminating divergences could be further training the neural networks using data collected from the attractor or fine-tuning the transient times at the training stage.
Another discrepancy that is best seen in Figure 9J and Figure 10J are predicted peaks appearing in locations that are shifted in space, or mirrored locations with respect to the reference solutions.These issues are due to the difficulties with reconstructing discretesymmetry-reduced trajectories, which relies on the numerical estimation of the partial derivative u ′ t .As stated in Section 3.1, this is a source of numerical error since it is obtained via finite differences.In my initial tests of the symmetry aligning transformations by applying them to the symmetry-reduced trajectories ξ(t), errors such as those seen in Figure 9J and Figure 10J were not present.However, when trajectories were obtained by models via decoding the linearly-evolved latent-space curves, I observed momentary jumps in the decoded trajectories, which in turn resulted in selection of false symmetry copies.I would like to emphasize that if one is interested in predicting the domain-integrated observables, then this is no issue, since such observables are symmetry-invariant.For discrete symmetry reduction, I produced symmetry-reducing complex variables through an optimization over the attractor.In future iterations of the modeling approach presented here, I plan to explore carrying out such optimizations also locally in the state space, since the models that are local in the state space could feasibly employ symmetry-reducing transformations that are also local.Another idea would be abandoning discrete symmetry reduction altogether at the expense of having pre-periodic orbits become periodic only after two periods.This, however, would potentially introduce performance issues due to the exponential divergence of trajectories.
While I chose to present the prediction results with forecast windows equal to the periods of the orbits, I would like to note that this was a pedagogical choice rather than a necessity.One straightforward way of reducing the prediction errors is indeed would be reducing the prediction windows.In units of Lyapunov time t L ≈ 20.9, which I estimated via the Benettin algorithm [51], the periods of PO 1 and PO 2 are approximately 1.6t L and 2.1t L , respectively.It can be readily seen from Figure 8 that reducing the prediction window to, for instance, one Lyapunov time could eliminate the episodes with larger errors towards the end of the prediction windows indicated by the dashed and dotted lines.

Conclusion & Outlook
In this paper, I presented a reduced-order modelling tool that combines the knowledge of numerically exact periodic orbits with the neural networks to produce predictive models of spatiotemporal chaos in the Kuramoto-Sivashinsky system.My primary goal was to demonstrate that such models can be significantly simple and exhibit some degree of interpretability if they are built around the periodic orbits.To illustrate the usefulness of understanding models, I implemented an orbit stabilization method wherein the control input was determined in the latent space.
Although Kuramoto-Sivashinsky system shares several similarities with the moderate-Reynolds-number turbulent flows with which I motivated this study, there are several challenges that need to be addressed when applying the tools presented here to such problems.The first of these is that moving from one spatial dimension to two-and three-leads to a substantial increase in the number of degrees of freedom as it grows geometrically with the spatial dimension.When this becomes prohibitive for neural network tools, one could include a PCA projection between the symmetry reduction and encoder in Figure 3 as an initial dimensionality reduction step to ease the computational burden.While the results presented here showed that the two periodic orbits were sufficient to capture nearly 50% of the time evolution, it would be too optimistic to expect such a simple phenomenology in higher-dimensional settings.Thus, when moving towards two-and three-dimensional fluid flows, one should be ready to face the difficulties of finding many periodic orbits in higherdimensional settings.Such future efforts could presumably benefit from the adjoint-based periodic orbit search tools [52,53] which typically exhibit better convergence behavior in comparison to the Newton-based methods.
An obvious difficulty of the present method in comparison to the other data-driven model reduction tools is the need for finding periodic orbits, which renders it impossible to apply to observational data collected from systems without detailed models.One potential way of eliminating this need could be utilizing DMD to approximate nearly periodic dynamics borrowing ideas from Refs.[54,55], where it was shown that DMD can be used to approximate unstable (relative) periodic orbits in shear flow simulations.Another way of eliminating the dependence on detailed models could be developing clustering methods that incorporate dynamics by considering bundles of trajectories as opposed to states sampled on an attractor.These are the research directions which I hope to explore in the future.

A Symmetry reduction of the Kuramoto-Sivashinsky system
In Fourier space, the translations T δx u(x, t) = u(x − δx, t) correspond to the rotations R(θ)v k = e ikθ v k where θ = 2πδx/L.Let ϕ 2 (t) = ∠v 2 (t) be the complex phase of the second Fourier mode, it is straightforward to confirm that the transformation sets the phase of the v ′ k to π/2.Although its exact value is irrelevant at this stage, setting this to π/2 is going to be convenient later.Unlike fixing the phase of the first Fourier mode, however, the transformation (15) does not fully eliminate the translations, since R(π)v ′ 2 = v 2 also satisfies the phase condition, whereas the odd modes pick up a phase of e iπ .In other words, the reduced dynamics in the space of v ′ k have a π-rotation symmetry which can be expressed as This multiplicity is easily eliminated by the 2-to-1 transformation where ϕ ′ 1 = ∠v ′ 1 .Noting that ϕ ′ 1 → ϕ ′ 1 + π under the action of R(π), one can confirm that v ′′ k are invariant under this symmetry.The discrete symmetry reducing transformation ( 17) is a simple application of the general method for reducing cyclic symmetries which will be presented in a separate publication [56].The final remaining symmetry of the system is the reflection and in order to reduce it, one first needs to find its action on the translationinvariant modes v ′′ k .As stated earlier, the phase condition ϕ ′ 2 = π/2 which ensures that v ′ 2 is purely imaginary is convenient since it is invariant under the action of reflection Mv k = −v * k , where * indicates complex conjugation.Thus, the action of M on the modes v ′ k is the same as that on v k .The π-reducing transformation (17), however, introduces an additional phase factor and which yields the action of M on v ′′ k as In order to follow the same recipe as the reduction of the rotation-by-π symmetry, one first needs to express the symmetry action as a complex phase e iπ .This is achieved grouping the symmetry-invariant parts Re v ′′ where c k are coefficients and use its phase ϕ M = ∠ ŵ to transform w − k as ŵ− k = e iϕ M ŵ− k , which are invariant under M.For the results of this paper, I determined c k by maximizing the amplitude of ŵ (19) in order to avoid fast phase oscillations.These coefficients along with the python implementations of forward and backward transformations are openly available in the code repository [50] accompanying this paper.

Figure 1 :
Figure 1: A spatiotemporally chaotic trajectory visualized as a spacetime plot where the amplitude of u(x, t) is color coded.The colorbar shown on the right applies to all space-time plots of this paper.

Figure 2 :
Figure 2: (A) A chaotic trajectory (gray) spanning a time interval of t ∈ [0, 10 4 ] along with two periodic orbits PO 1 and PO 2 visualized as projections from the 30-dimensional symmetry-reduced state space onto the leading three principal components P C 1,2,3 of the attractor.(B, C) Space-time visualizations of PO 1 (B) and PO 2 (C) for a time interval t ∈ [0, 50].The periods T 1 = 32.80 and T 2 = 43.61 of the relative periodic orbits are indicated by the horizontal dashed lines in B and C.

Figure 3 :
Figure 3: Block diagram illustrating the proposed model architecture.

Figure 5 :
Figure 5: Training data visualized as projections onto the principal components of the PO 1 (A) and the PO 2 (B).

Figure 6 :
Figure 6: Top (A-E) and bottom (F-J) illustrates two tests where the model predicts forward time evolution of the initial states in the vicinity of the periodic orbit PO 1 .The test trajectories (orange, thick) are visualized as projections onto the leading three principal components of the periodic orbit in the symmetry reduced state space (A and F), in the latent space (B and G) and on the (P, KE)-plane along with the corresponding predictions (green) and the periodic orbit (dashed blue).(D, I) Space-time visualizations of the test trajectories.(E, J) Space-time visualizations of the predictions.

Figure 7 :
Figure 7: Top (A-E) and bottom (F-J) illustrates two tests where the model predicts forward time evolution of the initial states in the vicinity of the periodic orbit PO 2 .The test trajectories (orange, thick) are visualized as projections onto the leading three principal components of the periodic orbit in the symmetry reduced state space (A and F), in the latent space (B and G) and on the (P, KE)-plane along with the corresponding predictions (green) and the periodic orbit (dashed blue).(D, I) Space-time visualizations of the test trajectories.(E, J) Space-time visualizations of the predictions.

Figure 8 :
Figure 8: Prediction errors in the space of the scalar fields u(x, t) and symmetry-reduced states ξ(t) (top) as well as on the domain-integrated observables power input P and the rate of dissipation ε (bottom) as a function of time.Dashed and dotted vertical lines indicate the episodes shown in Figure 9 and Figure 10, respectively.

Figure 9 :
Figure 9: Predictions of trajectories in the vicinity of PO 1 .Top (A-E) and bottom (F-J) rows illustrates two different cases corresponding to the consecutive trajectory segments along with their model predictions as PCA projection (A, F), in the latent space (B, G), on the (P, KE)-plane (C, H), and as space-time plots (D, E, I, and J).

Figure 10 :
Figure 10: Predictions of trajectories in the vicinity of PO 2 .Top (A-E) and bottom (F-J) rows illustrates two different cases corresponding to distinct trajectory segments along with their model predictions as PCA projection (A, F), in the latent space (B, G), on the (P, KE)plane (C, H), and as space-time plots (D, E, I, and J).

Figure 11 :
Figure 11: Stabilization of the periodic orbit via control.The space-time plots (A and B) show the periodic time-evolution of the system starting from the same initial condition as Figure 1 for t ∈ [0, 400].Projection of the symmetry-reduced trajectories onto the principal components of the periodic orbit (C): Chaotic trajectory for t ∈ [0, 200] (orange), Controlled trajectory for t ∈ [0, 200] (green), and the periodic orbit PO 1 (dashed blue).Time series of the instantaneous rate of dissipation (ε) and power input (P ) of the chaotic (orange) solution and the solution of the controlled system (green).

Table 1 :
Leading non-marginal Floquet multipliers and exponents of the periodic orbits.