Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds
We develop a methodology to construct low-dimensional predictive models from data sets representing essentially nonlinear (or non-linearizable) dynamical systems with a hyperbolic linear part that are subject to external forcing with finitely many frequencies. Our data-driven, sparse, nonlinear models are obtained as extended normal forms of the reduced dynamics on low-dimensional, attracting spectral submanifolds (SSMs) of the dynamical system. We illustrate the power of data-driven SSM reduction on high-dimensional numerical data sets and experimental measurements involving beam oscillations, vortex shedding and sloshing in a water tank. We find that SSM reduction trained on unforced data also predicts nonlinear response accurately under additional external forcing. Current data-driven modelling techniques perform reliably on linear systems or on those that can be linearized. Cenedese et al. develop a data-based reduced modeling method for non-linear, high-dimensional physical systems. Their models reconstruct and predict the dynamics of the full physical system.
/static/images/favicons/nature/apple-touch-icon-f39cb19454.png
Nature
Spectral submanifolds and their reduced dynamics

A recent result in dynamical systems is that all eigenspaces (or spectral subspaces) of linearized systems admit unique nonlinear continuations under well-defined mathematical conditions. Specifically, spectral submanifolds (SSMs), as defined by29, are the unique smoothest invariant manifolds that serve as nonlinear extensions of spectral subspaces under the addition of nonlinearities to a linear system. The SSM formulation and terminology we use here is due to29; the Methods section “Existence of SSMs” discusses the history of these results and further technical details.

We consider n-dimensional dynamical systems of the form

$$\dot{{{{{{{{\bf{x}}}}}}}}}={{{{{{{\bf{A}}}}}}}}{{{{{{{\bf{x}}}}}}}}+{{{{{{{{\bf{f}}}}}}}}}_{0}({{{{{{{\bf{x}}}}}}}})+\epsilon {{{{{{{{\bf{f}}}}}}}}}_{1}({{{{{{{\bf{x}}}}}}}},{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon ),\qquad {{{{{{{{\bf{f}}}}}}}}}_{0}({{{{{{{\bf{x}}}}}}}})={{{{{{{\mathcal{O}}}}}}}}({\left\vert{{{{{{{\bf{x}}}}}}}}\right\vert}^{2}),\qquad 0\le \epsilon \ll 1,$$

(1)

with a constant matrix \({{{{{{{\bf{A}}}}}}}}\in {{\mathbb{R}}}^{n\times n},\) and with class Cr functions \({{{{{{{{\bf{f}}}}}}}}}_{0}:{{{{{{{\mathcal{U}}}}}}}}\to {{\mathbb{R}}}^{n}\) and \({{{{{{{{\bf{f}}}}}}}}}_{1}:{{{{{{{\mathcal{U}}}}}}}}\times {{\mathbb{T}}}^{\ell }\to {{\mathbb{R}}}^{n}\), where \({{\mathbb{T}}}^{\ell }={S}^{1}\times \ldots \times {S}^{1}\) is the ℓ-dimensional torus. The elements of the frequency vector \({{{{{{{\boldsymbol{\Omega }}}}}}}}{\mathbb{\in }}{{\mathbb{R}}}^{\ell }\) are rationally independent, and hence the function f1 is quasiperiodic in time. The assumed degree of smoothness for the right-hand side of (1) is \(r\in {{\mathbb{N}}}^{+}\cup \left\{\infty ,a\right\}\), with a referring to analytic. The small parameter ϵ signals that the forcing in system (1) is moderate so that the structure of the autonomous part is still relevant for the full system dynamics. Rigorous mathematical results on SSMs are proven for small enough ϵ, but continue to hold in practice for larger values of ϵ as well, as we will see in examples. Note that eq. (1) describes equations of motions of physical oscillatory systems. It does not cover phenomenological models of phase oscillators, such as the Kuramoto model35.

The eigenvalues \({\lambda }_{j}={\alpha }_{j}+{{{{{{{\rm{i}}}}}}}}{\omega }_{j}\in {\mathbb{C}}\) of A, with multiplicities counted, are ordered based on their real parts, \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{j}\), as

$${{{{{{{\rm{Re}}}}}}}}{\lambda }_{n}\le {{{{{{{\rm{Re}}}}}}}}{\lambda }_{n-1}\le \ldots \ldots \le {{{{{{{\rm{Re}}}}}}}}{\lambda }_{1}.$$

(2)

Their corresponding real modal subspaces (or eigenspaces), \({E}_{j}\subset {{\mathbb{R}}}^{n}\), are spanned by the imaginary and real parts of the corresponding eigenvectors and generalized eigenvectors of A. To analyze typical systems, we assume that \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{j}={\alpha }_{j}\,\ne\, 0\) holds for all eigenvalues, i.e., x = 0 is a hyperbolic fixed point for ϵ = 0.

A spectral subspace \({E}_{{j}_{1},\ldots ,{j}_{q}}\) is a direct sum

$${E}_{{j}_{1},\ldots ,{j}_{q}}={E}_{{j}_{1}}\oplus {E}_{{j}_{2}}\oplus \ldots \oplus {E}_{{j}_{q}}$$

(3)

of an arbitrary collection of modal subspaces, which is always an invariant subspace for the linear part of the dynamics in (1). Classic examples of spectral subspaces are the stable and unstable subspaces, comprising all modal subspaces with \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{k} \, < \, 0\) and \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{k} \, > \, 0\), respectively. Projections of the linearized system onto the nested hierarchy of slow spectral subspaces,

$${E}^{1}\subset {E}^{2}\subset {E}^{3}\subset \ldots ,\qquad {E}^{k}:= {E}_{1,\ldots ,k},\quad k=1,\ldots ,n,$$

(4)

provide exact reduced-order models for the linearized dynamics over an increasing number of time scales under increasing k, as sketched in panel (a) of Fig. 2. This is why a Galerkin projection onto Ek is an exact model reduction procedure for linear systems, whose accuracy can be increased by increasing k. A fundamental question is whether nonlinear analogues of spectral subspaces continue to organize the dynamics under the addition of nonlinear and time-dependent terms in the full system (1).

Let us fix a specific spectral subspace \(E={E}_{{j}_{1},\ldots ,{j}_{q}}\) within either the stable or the unstable subspace. If E is non-resonant (i.e., no nonnegative, low-order, integer linear combination of the spectrum of A∣E is contained in the spectrum of A outside E), then E has infinitely many nonlinear continuations in the system (1) for ϵ small enough29. These invariant manifolds are of smoothness class CΣ(E), with the spectral quotient Σ(E) measuring the ration of the fastest decay exponent outside E to the slowest decay exponent inside E (see eq. (13) of the Methods section “Existence of SSMs”). All such manifolds are tangent to E for ϵ = 0, have the same quasiperiodic time dependence as f1 does and have a dimension equal to that of E.

Of these infinitely may invariant manifolds, however, there will be a unique smoothest one, the spectral submanifold (SSM) of E, denoted W(E, Ωt; ϵ). This manifold is Cr smooth if r > Σ(E) and can therefore be approximated more accurately than the other infinitely many nonlinear continuations of E. In particular, SSMs have convergent Taylor expansions if the dynamical system (1) is analytic (r = a). Then the reduced dynamics on a slow SSM, Ek, can be approximated with arbitrarily high accuracy using arbitrarily high-order Taylor expansions, without ever increasing the dimension of Ek, see panel (b) of Fig. 2. Such an approximation for dynamical systems with known governing equations is now available for any required order of accuracy via the open-source MATLAB® package SSMTool36. In contrast, reduced models obtained from projection-based procedures can only be improved by increasing their dimensions.

The nearby coexisting stationary states in Fig. 2 happen to be contained in the SSM. In specific examples, however, these states may also be off the SSM, contained instead in one of the infinitely many additional nonlinear continuations, \(\tilde{W}(E,{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\), of the spectral subspace E. The Taylor expansion of the dynamics on \(\tilde{W}(E,{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\) and W(E, Ωt; ϵ) are, however, identical up to order Σ(E). Therefore, the reduced models we will compute on the SSM W(E, Ωt; ϵ) also correctly capture the nearby stationary states on \(\tilde{W}(E,{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\), as long as the polynomial order of the model stays below Σ(E). In large physical systems, this represents no limitation, given that Σ(E) ≫ 1.

Embedding SSMs via generic observables

If at least some of the real parts of the eigenvalues in (2) are negative, then longer-term trajectory data for system (1) will be close to an attracting SSM, as illustrated in panel (b) of Fig. 2. This is certainly the case for data from experiments that are run until a nontrivial, attracting steady state emerges, see, e.g., in panel (e) of Fig. 1. Measurements of trajectories in the full phase space, however, are seldom available from such experiments. Hence, if data about system (1) is only available from observables, the construction of SSMs and their reduced dynamics has to be carried out in the space of those observables.

An extended version of Whitney’s embedding theorem guarantees that almost all (in the sense of prevalence) smooth observable vectors \({{{{{{{\bf{y}}}}}}}}({{{{{{{\bf{x}}}}}}}})=({y}_{1}({{{{{{{\bf{x}}}}}}}}),,...,{y}_{p}({{{{{{{\bf{x}}}}}}}}))\in {{\mathbb{R}}}^{p}\) provide an embedding of a compact subset \({{{{{{{\mathcal{C}}}}}}}}\subset W(E,{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\) of a d-dimensional SSM, W(E, Ωt; ϵ), into the observable space \({{\mathbb{R}}}^{p}\) for high enough p. Specifically, if we have p > 2(d + ℓ) simultaneous and independent continuous measurements, y(x), of the p observables, then almost all maps \({{{{{{{\bf{y}}}}}}}}:{{{{{{{\mathcal{C}}}}}}}}\to {{\mathbb{R}}}^{p}\) are embeddings of \({{\mbox{}}}{{{{{{{\mathcal{C}}}}}}}}{{\mbox{}}}\)37, and hence the top right plot of Fig. 3 is applicable with probability one.

In practice, we may not have access to p > 2(d + ℓ) independent observables and hence cannot invoke Whitney’s theorem. In that case, we invoke the Takens delay embedding theorem38, which covers observable vectors built from p uniformly sampled, consecutive measured instances of a single observable. More precisely, if s(t) is a generic scalar quantity measured at times Δt apart, then the observable vector for delay-embedding is formed as \({{{{{{{\bf{y}}}}}}}}(t)=\left(s(t),s(t+{{\Delta }}t),...,s\left(t+(p-1){{\Delta }}t\right)\right)\in {{\mathbb{R}}}^{p}\). We discuss the embedding, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\subset {{\mathbb{R}}}^{p}\), of an autonomous SSM, W(E, Ωt0; 0), in the observable space \({{\mathbb{R}}}^{p}\) in more detail in the Methods section “Embedding the SSM in the observable space”.

Data-driven extended normal forms on SSMs

Once the embedded SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), is identified in the observable space, we seek to learn the reduced dynamics on \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\). An emerging requirement for learning nonlinear models from data has been model sparsity4, without which the learning process would be highly sensitive. The dynamics on \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), however, is inherently nonsparse, which suggests that we learn its Poincaré normal form39 instead. This classic normal form is the simplest polynomial form to which the dynamics can be brought via successive, near-identity polynomial transformations of increasing order.

Near the origin on a slow SSM, however, this simplest polynomial form is just the restriction of the linear part of system (1) to \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), as long as infinitely many nonresonance conditions are satisfied for the operator A40. The Poincaré normal form on \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) would, therefore, only capture the low-amplitude, linearized part of the slow SSM dynamics.

To construct an SSM-reduced model for non-linearizable dynamics, we use extended normal forms. This idea is motivated by normal forms used in the study of bifurcations of equilibria on center manifolds depending on parameters33,41. In that setting, the normal form transformation is constructed at the bifurcation point where the system is non-linearizable by definition. The same transformation is then used away from bifurcations, even though the normal form of the system would be linear there. One, therefore, gives up the maximal possible simplicity of the normal form but gains a larger domain on which the normal form transformation is invertible and hence captures truly nonlinear dynamics. In our setting, there is no bifurcation at x = 0, but we nevertheless construct our normal form transformation as if the eigenvalues corresponding to the slow subspace E were purely imaginary. This procedure leaves additional, near-resonant terms in the SSM-reduced normal form, enhancing the domain on which the transformation is invertible and hence the normal form is valid.

We determine the normal form coefficients directly from data via the minimization of a conjugacy error (see the Methods section). This least-square minimization procedure renders simultaneously the best-fitting normal form coefficients and the best fitting normal form transformation. As we will find in a specific example, this data-driven procedure can yield accurate reduced models even beyond the formal domain of convergence of equation-driven normal forms.

The simplest extended normal form on a slow SSM of an oscillatory system arises when the underlying spectral subspace E corresponds to a pair of complex conjugate eigenvalues. Writing in polar coordinates and truncating at cubic order,42 finds this normal form on the corresponding two-dimensional, autonomous SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), to be

$$\begin{array}{ll}\dot{\rho }=\,{\alpha }_{0}\rho +\beta {\rho }^{3},\!\!\!\\ \dot{\theta }=\,{\omega }_{0}+\gamma {\rho }^{2}.\end{array}$$

(5)

This equation is also known as the Stuart–Landau equation arising in the unfolding of a Hopf bifurcation43,44,45.

The dynamics of (5) is characteristically nonlinearizable when α0β < 0, given that a limit cycle coexists with the ρ = 0 fixed point in that case. Further coexisting steady states will arise when forcing is added to the system, as we discuss in the next section. We note that the cubic normal form on two-dimensional SSMs has also been approximated from data in46. That non-sparse procedure fits the full observer dynamics to a low-dimensional, discrete polynomial dynamical system, then performs an analytic SSM reduction and a classic normal form transformation on the SSM.

For higher accuracy, the extended normal form on an oscillatory SSM of dimension 2m is of the form

$$\begin{array}{c}{\dot{\rho }}_{j}={\alpha }_{j}({{{{{{{\boldsymbol{\rho }}}}}}}},{{{{{{{\boldsymbol{\theta }}}}}}}}){\rho }_{j},\\ \!\!\!\!\!{\dot{\theta }}_{j}={\omega }_{j}({{{{{{{\boldsymbol{\rho }}}}}}}},{{{{{{{\boldsymbol{\theta }}}}}}}}),\end{array}\,\,\,\,\,\,\,j=1,2,...m,\,\,\,\,\,\,{{{{{{{\boldsymbol{\rho }}}}}}}}\in {{\mathbb{R}}}_{+}^{m},\,\,\,\,\,\,{{{{{{{\boldsymbol{\theta }}}}}}}}\in {{\mathbb{T}}}^{m}.$$

(6)

If the linearized frequencies are nonresonant, then the functions αj and ωj only depend on ρ42. Our numerical procedure determines these functions up to the necessary order that ensures a required accuracy for the reduced-order model on the SSM. This is illustrated schematically for a four-dimensional slow SSM (m = 2) in the bottom right plot of Fig. 3.

Predicting forced dynamics from unforced data

With the normalized reduced dynamics (6) on the embedded SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), at hand, we can also make predictions for the dynamics of the embedded quasiperiodic SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\), of the full system (1). This forced SSM is guaranteed to be an \({{{{{{{\mathcal{O}}}}}}}}(\epsilon )\)Cr-close perturbation of \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) for moderate external forcing amplitudes. A strict proof of this fact is available for small enough ϵ > 029, but as our examples will illustrate, the smooth persistence of the SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\), generally holds for all moderate ϵ values in practice. Such moderate forcing is highly relevant in a number of technological settings, including system identification in structural dynamics and fluid-structure interactions, where the forcing must be moderate to preserve the integrity of the structure.

We discuss the general extended normal form on \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\) in the Methods section “SSM dynamics via extended normal forms”. In the simplest and most frequent special case, the external forcing is periodic (ℓ = 1) and \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{\Omega }}t)\) is the embedding of the slowest, two-dimensional SSM corresponding to a pair of complex conjugate eigenvalues. Using the modal forcing amplitude f1,1 and modal phase shift ϕ1,1 in the general normal form (25)47, introduces the new phase coordinate ψ = θ − Ωt − ϕ1,1 and lets f = f1,1, α = α1, ω = ω1 to obtain the planar, autonomous, extended normal form on \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{\Omega }}t)\) as

$$\dot{\rho } =\alpha (\rho )\rho +f\sin \psi ,\\ \dot{\psi } =\omega (\rho )-{{\Omega }}+\frac{f}{\rho }\cos \psi$$

(7)

at leading order in ϵ. All stable and unstable periodic responses on the SSM are fixed points of system (7), with their amplitudes ρ0 and phases ψ0 satisfying the equations

$${{\Omega }}=\omega ({\rho }_{0})\pm \sqrt{\frac{{f}^{2}}{{\rho }_{0}^{2}}-{\alpha }^{2}({\rho }_{0})},\quad {\psi }_{0}={\tan }^{-1}\left[\frac{\alpha \left({\rho }_{0}\right)}{\omega \left({\rho }_{0}\right)-{{\Omega }}}\right].$$

(8)

The first analytic formula in (8) predicts the forced response curve (FRC) of system (1), i.e., the relationship between response amplitude, forcing amplitude and forcing frequency, from the terms α(ρ) and ω(ρ) of the extended normal form of the autonomous SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\). These terms are constructed from trajectories of the unforced system, thus eq. (8) predicts the behavior of a nonlinearizable dynamical system under forcing based solely on unforced training data. The stability of the predicted periodic response follows from a simple linear analysis at the corresponding fixed point of the ODE (7). The first formula in (8) also contains another frequently used notion of nonlinear vibration analysis, the dissipative backbone curve ω(ρ), which describes the instantaneous amplitude-frequency relation along freely decaying vibrations within the SSM.

As we will also show in examples, our unforced model-based predictions for forced periodic response (see the Methods section “Prediction of forced response from unforced training data”) are confirmed by numerical simulation or dedicated laboratory experiments on forced systems.

Examples

We now illustrate data-driven, SSM-based modeling and prediction on several numerical and experimental data sets describing non-linearizable physical systems. Further applications are described in48. Both the numerical and the experimental data sets were initialized without knowledge of the exact SSM. All our computations have been carried out by the publicly available MATLAB® package, SSMLearn, whose repository also contains further examples not discussed here. The main algorithm behind SSMLearn is illustrated in Fig. 3, with more detail given in the Methods section “Summary of the algorithm”.

To quantify the errors of an SSM-based reduced model, we use the normalized mean-trajectory-error (NMTE). For P observations of the observable vector yj and their model-based reconstructions, \({\hat{{{{{{{{\bf{y}}}}}}}}}}_{j}\), this modeling error is defined as

$${{{{{{{\rm{NMTE}}}}}}}}=\frac{1}{\parallel \underline{{{{{{{{\bf{y}}}}}}}}}\parallel }\frac{1}{P}\mathop{\sum }\limits_{j=1}^{P}\parallel {{{{{{{{\bf{y}}}}}}}}}_{j}-{\hat{{{{{{{{\bf{y}}}}}}}}}}_{j}\parallel\!.$$

(9)

Here \(\underline{{{{{{{{\bf{y}}}}}}}}}\) is a relevant normalization vector, such as the data point with the largest norm. When validating the reduced dynamics for a given testing trajectory, we run the reduced model from the same initial condition for the comparison. Increasing the order of the truncated normal form polynomials in eq. (6) generally reduces the NMTE error to any required level but excessively small errors can lead to overfitting. In our examples, we will be allowing model errors in the order of 1% − 4% to avoid overfitting.

As a first example, we consider a finite-element discretization of a von Kármán beam with clamped-clamped boundary conditions49, shown in panel (a) of Fig. 4. In contrast to the classic Euler-Bernoulli beam, the von Kármán model captures moderate deformations by including a nonlinear, quadratic term in the kinematics. We first construct a 33 degree-of-freedom, damped, unforced finite element model (i.e., n = 66 and ϵ = 0 in eq. (1)) for an aluminum beam of length 1 [m], width 5 [cm], thickness 2 [cm] and material damping modulus 106 [] (see the Supplementary Information for more detail).

Our objective is to learn from numerically generated trajectory data the reduced dynamics on the slowest, two-dimensional SSM, W(E1), of the system, defined over the slowest two-dimensional (d = 2) eigenspace E1 of the linear part. To do so, we generate two trajectories starting from initial beam deflections caused by static loading of 12 [kN] and 14 [kN] at the midpoint, as shown in panel (a) of Fig. 4. The latter trajectory, shown in panel (b) of Fig. 4, is used for training, the other for testing. Along the trajectories, we select our single observable s(t) to be the midpoint displacement of the beam.

The beam equations are analytic (r = a), and hence the SSM, W(E1), admits a convergent Taylor expansion near the origin. The minimal embedding dimension for the two-dimensional, W(E1), as required by Whitney’s theorem, is p = 5, which is not satisfied by our single scalar observable s(t). We therefore employ delay-embedding using \({{{{{{{\bf{y}}}}}}}}(t)=\left(s(t),s(t+{{\Delta }}t),\ldots ,s(t+4{{\Delta }}t)\right)\) with Δt = 0.0955 [ms]. By Takens’s theorem, this delayed observable embeds the SSM in \({{\mathbb{R}}}^{5}\) with probability one.

A projection of the embedded SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\in {{\mathbb{R}}}^{5},\) onto three coordinates is shown in panel (c) of Fig. 4. On \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), SSMLearn returns the 7th-order extended normal form

$$\dot{\rho } =\alpha (\rho )\rho ,\quad \alpha (\rho )=-3.02-5.79{\rho }^{2}+57.5{\rho }^{4}-191{\rho }^{6},\\ \dot{\theta } =\omega (\rho ),\,\,\quad \omega (\rho )=658+577{\rho }^{2}-347{\rho }^{4}-387{\rho }^{6},$$

(10)

to achieve our preset reconstruction error bar of 3% on the test trajectory (NMTE = 0.027), shown in panel (d) of Fig. 4.

We now use the model (10), trained on a single decaying trajectory, to predict the forced response of the beam for various forcing amplitudes and frequencies in closed form. We will then compare these predictions with analytic forced response computations for the forced SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{\Omega }}t)\), obtained from SSMTool36 and with numerical simulations of the damped-forced beam. The periodic forcing is applied at the midpoint node; the Taylor expansion order in SSMTool for the analytically computed dynamics on \({{{{{{{{\mathcal{M}}}}}}}}}_{\epsilon }({{\Omega }}t)\) is set to 7, as in (10). Panel (e) of Fig. 4 shows the FRCs (green) and the backbone curve (blue) predicted by SSMLearn based on formula (8) from the single unforced trajectory in panel (b) of Fig. 4. To obtain the relevant forcing amplitudes f in the delay-observable space, we have followed the calibration procedure described in the Methods section “Prediction of forced response from unforced training data” for the forcing values \(\left|\epsilon {{{{{{{{\bf{f}}}}}}}}}_{1}\right|=15,45,95\) [N] at the single forcing frequency Ω = 103.5 [Hz]. Recall that coexisting stable (solid lines) and unstable (dashed lines) periodic orbits along the same FRC are hallmarks of non-linearizable dynamics and hence cannot be captured by the model reduction techniques we reviewed in the Introduction for linearizable systems.

The data-based prediction for the FRCs agrees with the analytic FRCs for low forcing amplitudes but departs from it for higher amplitudes. Remarkably, as the numerical simulations (red) confirm, the data-based FRC is the correct one. The discrepancy between the two FRCs for large amplitudes only starts decreasing under substantially higher-order Taylor series approximations used in SSMTool (see the Supplementary Information). This suggests the use of the data-based approach for this class of problems even if the exact equations of motion are available.

As a second example, we consider the classic problem of vortex shedding behind a cylinder8. Our input data for SSM-based reduced modeling are the velocity and pressure fields over a planar, open fluid domain with a hole representing the cylinder section, as shown in panel (a) of Fig. 5. The boundary conditions are no-slip on the circular inner boundary, standard outflow on the outer boundary at the right side, and fixed horizontal velocity on the three remaining sides50. The Reynolds number for this problem is the ratio between the cylinder diameter times the inflow velocity and the kinematic viscosity of the fluid.

Available studies8,50,51 report that, at low Reynolds number, the two-dimensional unstable manifold, Wu(SS), of the wake-type steady solution, SS, in panel (b) of Fig. 5 connects SS to the limit cycle shown in panel (c) of Fig. 5. Here we evaluate the performance of SSMLearn on learning this unstable manifold as an SSM, along with its reduced dynamics, from trajectory data at Reynolds number equal to 70. For this SSM, we again have d = 2 and r = a, as in our previous example. There is no external forcing in this problem, and hence we have ϵ = 0 in eq. (1). In contrast to prior studies that often consider a limited number of observables8,51,52, here we select the full phase space of the discretized Navier-Stokes simulation to be the observable space for illustration, which yields n = p = 76, 876 in eq. (1). We generate nine trajectories numerically, eight of which will be used for training and one for testing the SSM-based model.

The nine initial conditions of our input trajectory data are small perturbations from the wake-type steady solution along its unstable directions, equally spaced on a small amplitude circle on this unstable plane. All nine trajectories quickly converge to the unstable manifold and then to the limit cycle representing periodic vortex shedding.

We choose to parametrize the SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}={W}^{u}(SS)\), with two leading POD modes of the limit cycle, which have been used in earlier studies for this problem. The training trajectories projected onto these two POD modes are shown in panel (d) of Fig. 5. To limit the modeling error (9) to less than NMTE = 1%, SSMLearn requires a polynomial order of 18 in the SSM computations. For this order, our approach can accommodate the strong mode deformation observed for this problem51, manifested by a fold of the SSM over the unstable eigenspace in panel (f) of Fig. 5. Panel (g) of Fig. 5 shows the strongly nonlinear geometry of \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) projected to the observable subspace formed by the velocities and the pressure of a probe point in the wake.

To capture the SSM-reduced dynamics with acceptable accuracy, we need to compute the extended normal form up to order 11, obtaining

$$\begin{array}{c}\dot{\rho }=\alpha (\rho )\rho =0.0584\rho -0.479{\rho }^{3}+1.27{\rho }^{5}+6.80{\rho }^{7}-58.9{\rho }^{9}+108{\rho }^{11},\\ \!\!\!\!\!\!\!\!\!\dot{\theta }=\omega (\rho )\,=0.553+0.441{\rho }^{2}-3.38{\rho }^{4}+55.5{\rho }^{6}-321{\rho }^{8}+626{\rho }^{10}.\end{array}$$

(11)

To describe a transition qualitatively from a fixed point to a limit cycle, the reduced-order dynamical model should be at least of cubic order51. Capturing the qualitative behavior (i.e., the unstable fixed point and the stable limit cycle), however, does not imply a low NMTE error for the model. Indeed, the data-driven cubic normal form for this example gives a reconstruction error of NMTE = 117% normalized over the limit cycle amplitude, mainly arising from an out-of-phase convergence to the limit cycle along the testing trajectory. In contrast, the \({{{{{{{\mathcal{O}}}}}}}}(11)\) normal form in eq. (11) reduced this error drastically to NMTE = 3.86% on the testing trajectory, as shown in panel (e) of Fig. 5.

We show in Section 1.2.3 of the Supplementary Information that for comparable accuracy, the Sparse Identification of Nonlinear DYnamics (SINDy) approach of4 returns non-sparse nonlinear models for this example. Similarly, while the DMD13 can achieve highly accurate curve-fitting on the available training trajectories with a high-dimensional linear model, that model only captures linearizable dynamics near the origin. As a consequence, its trajectories grow without bound over longer integration times and hence fail to capture the limit cycle.

As a third example, we consider fluid oscillations in a tank, which exhibit highly nonlinear characteristics53. To describe such non-linearizable softening effects observed in the sloshing motion of surface waves, Duffing-type models have been proposed54. While amplitude variations observed in forced experiments can be fitted to forced softening Duffing equations, nonlinear damping remains a challenge to identify55.

The experiments we use to construct an SSM-reduced nonlinear model for sloshing were performed in a rectangular tank of width 500 [mm] and depth 50 [mm], partially filled with water up to a height of 400 [mm], as shown in panel (a) of Fig. 6. The tank was mounted on a platform excited harmonically by a motor. The surface level was detected via image processing from a monochrome camera. As an observable s(t) we used the horizontal position of the computed center of mass of the water at each time instant, normalized by the tank width. This physically meaningful scalar is robust with respect to image evaluation errors55.

We identify the unforced nonlinear behavior of the system from data obtained in resonance decay experiments56. In those experiments (as in Fig. panel (a) of 6, but with a shaker instead of a motor), once a periodic steady state is reached under periodic horizontal shaking of the tank, the shaker is turned off and the decaying sloshing is recorded. We show such a decaying observable trajectory (orange line) in panel (b) of Fig. 6, with the shaker switched off slightly before zero time. This damped oscillation is close, by construction, to the two-dimensional, slowest SSM of the system. We use three such decaying observer trajectories (two for training and one for model testing) for the construction of a two-dimensional (d = 2), autonomous, SSM-based reduced-order model for s(t). For delay embedding dimension, we again pick p = 5, the minimal value guaranteed to be generically correct for embedding the SSM by Takens’s theorem. The delay used in sampling s(t) is Δt = 0.033 [s]. For this input and for a maximal reconstruction error of 2%, SSMLearn identifies a nearly flat SSM in the delayed observable space–see panel (c) of Fig. 6–with a cubic extended normal form

$$\dot{\rho }=-0.063179\rho -0.041214{\rho }^{3},\quad \dot{\theta }=7.8144-1.5506{\rho }^{2}.$$

(12)

This lowest-order, Stuart–Landau-type normal form, cf. (5), already constitutes an accurate reduced-order model with NMTE = 1.88% on the testing data set, see panel (b) of Fig. 6. The amplitude-dependent nonlinear damping, α(ρ), provided by this model is plotted in panel (d) of Fig. 6 with respect to the physical amplitude.

In another set of experiments with the setup of panel (a) Fig. 6, steady states of periodically forced sloshing were measured in sweeps over a range of forcing frequencies under three different shaker amplitudes. As in the previous beam example, we identify the corresponding forcing amplitude, f, in (7) at the maximal amplitude response of each frequency sweep. Shown in panels (e, f) of Fig. 6, the closed-form predictions for FRCs from eq. (8) (solid lines) match closely the experimental FRCs (dots). Given the strong nonlinearity of the FRC, any prediction of this curve from a DMD-based model is bound to be vastly inaccurate, as we indeed show in Section 1.3 of the Supplementary information.

The phase ψ0 of the forced response relative to the forcing has been found difficult to fit to forced Duffing-type models55, but the present modeling methodology also predicts this phase accurately using the second expression in (8). The blue curve in panel (e) of Fig. 6 shows the backbone curve of decaying vibrations, which terminates at the highest amplitude occurring in the training data set. This plot therefore shows that the closed-form FRC predictions obtained from the SSM-based reduced model are also effective for response amplitudes outside the training range of the reduced model.

Existence of SSMs

In the context of rigid body dynamics, invariant manifolds providing generalizations of invariant spectral subspaces to nonlinear systems were first envisioned and formally constructed as nonlinear normal modes by58 (see59 for a recent review of related work). Later studies, however, pointed out the nonuniqueness of nonlinear normal modes in specific examples (60,61).

In the mathematics literature,62 obtained general results on the existence, smoothness and degree of uniqueness of such invariant manifolds for mappings on Banach spaces. These results use a special parameterization method to construct the manifolds even in evolutionary partial differential equations that admit a well-posed flow map in both time directions (see63 for a mechanics application). The results have been extended to a form applicable to dynamical systems with quasiperiodic time dependence64. An extensive account of the numerical implementation of the parametrization method with a focus on computing invariant tori and their whiskers in Hamiltonian systems is also available65,29 Discussed the existence of the SSM, W(E, Ωt; ϵ), depending on its absolute spectral quotient,

$${{\Sigma }}(E)={{{{{{{\rm{Int}}}}}}}}\,\left[\frac{\mathop{\max }\limits_{\lambda \in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{S})}| {{{{{{{\rm{Re}}}}}}}}\lambda | }{\mathop{\min }\limits_{{\lambda }_{e}\in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{E})}| {{{{{{{\rm{Re}}}}}}}}{\lambda }_{e}| }\right],$$

(13)

where Spect(A∣S) is the stable (unstable) spectrum of A if the SSM is stable (unstable). For a stable SSM, Σ(E) is the integer part of the quotient of the minimal real part in the spectrum of A and the maximal real part of the spectrum of A restricted to E.

Based on Σ(E), we call a d-dimensional spectral subspace E non-resonant if for any set \(\left({m}_{1},\ldots ,{m}_{d}\right)\) of nonnegative integers satisfying \(2\le \mathop{\sum }\nolimits_{j = 1}^{d}{m}_{j}\le {{\Sigma }}(E)\), the eigenvalues, λk, of A satisfy

$$\mathop{\sum }\limits_{j=1}^{d}{m}_{j}{{{{{{{\rm{Re}}}}}}}}{\lambda }_{j}\,\ne \,{{{{{{{\rm{Re}}}}}}}}{\lambda }_{k},\quad {\lambda }_{k}\in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}})-{{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{E}).$$

(14)

This condition only needs to be verified for resonance orders between 2 and Σ(E)64. In particular, a 1: 1 resonance between E1 and E2 is allowed if \(\dim {E}_{1}=\dim {E}_{2}=1\), in which case each strongly resonant spectral subspace gives rise to a unique nearby spectral submanifold.

If E violates the nonresonance condition (14), then E can be enlarged to a higher-dimensional spectral subspace until the nonresonance relationship (14) is satisfied. In the absence of external forcing (ϵ = 0), the nonresonance condition (14) can also be relaxed with the help of the relative spectral quotient,

$$\sigma (E)={{{{{{{\rm{Int}}}}}}}}\,\left[\frac{\mathop{\max }\limits_{\lambda \in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{S})-{{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{E})}| {{{{{{{\rm{Re}}}}}}}}\lambda | }{\mathop{\min }\limits_{{\lambda }_{e}\in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{E})}| {{{{{{{\rm{Re}}}}}}}}{\lambda }_{e}| }\right],$$

(15)

to the form

$$\mathop{\sum }\limits_{j=1}^{d}{m}_{j}{\lambda }_{j}\,\ne\, {\lambda }_{k},\quad {\lambda }_{k}\in {{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}})-{{{{{{{\rm{Spect}}}}}}}}({{{{{{{\bf{A}}}}}}}}{| }_{E}),\qquad 2\le \mathop{\sum }\limits_{j=1}^{d}{m}_{j}\le \sigma (E).$$

(16)

This is indeed a relaxation because condition (16) is only violated if both the real and the imaginary parts of eigenvalues involved are in the exact same resonance with each other. In contrast, (14) is already violated when the real parts are in resonance with each other.

If \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{1} < 0\) in eq. (2) and all Ek subspaces are nonresonant, then the nested set of slow spectral submanifolds,

$$W({E}^{1},{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\subset W({E}^{2},{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\subset W({E}^{3},{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\subset \ldots ,$$

gives a hierarchy of local attractors. All solutions in a vicinity of x = 0 approach the reduced dynamics on one of these attractors exponentially fast, as sketched in panel (b) of Fig. 2 for the ϵ = 0 limit. As we will see, non-linearizable dynamics tend to emerge on W(Ek, Ωt; ϵ) due to near-resonance between the linearized frequencies within Ek and the forcing frequencies Ω. The specific location of nontrivial steady states in W(Ek, Ωt; ϵ) is then determined by a balance between the nonlinearities, damping and forcing.

A resonant Ek subspace can be enlarged by adding the next \(k^{\prime}\) modal subspaces to it until \({E}^{k+k^{\prime} }\) in the hierarchy (4) becomes non-resonant and hence admits an SSM, \(W({E}^{k+k^{\prime} },{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\). This technical enlargement is also in agreement with the physical expectation that all interacting modes have to be included in an accurate reduced-order model. Finally, we note that SSMs are robust features of dynamical systems: they inherit smooth dependence of the vector field in (1) on parameters29.

For discrete-time dynamical systems of the form

$${{{{{{{{\bf{x}}}}}}}}}_{k+1}=\tilde{{{{{{{{\bf{A}}}}}}}}}{{{{{{{{\bf{x}}}}}}}}}_{k}+{\tilde{{{{{{{{\bf{f}}}}}}}}}}_{0}({{{{{{{{\bf{x}}}}}}}}}_{k})+\epsilon {\tilde{{{{{{{{\bf{f}}}}}}}}}}_{1}({{{{{{{{\bf{x}}}}}}}}}_{k},{{{{{{{{\boldsymbol{\phi }}}}}}}}}_{k};\epsilon ),\qquad {{{{{{{{\boldsymbol{\phi }}}}}}}}}_{k+1}={{{{{{{{\boldsymbol{\phi }}}}}}}}}_{k}+\tilde{{{{{{{{\boldsymbol{\Omega }}}}}}}}},$$

(17)

the above results on SSMs apply based on the eigenvalues μk of \(\tilde{{{{{{{{\bf{A}}}}}}}}}\). One simply needs to replace λk with \(\log {\mu }_{k}\) and \({{{{{{{\rm{Re}}}}}}}}{\lambda }_{k}\) with \(\log | {\mu }_{k}|\) in formulas (13)-(16)29.

We close by noting that in a neighborhood of an SSM, an invariant family of surfaces resembling the role of coordinate planes in a linear system exists66. This invariant spectral foliation (ISF) can, in principle, be used to generate a nonlinear analogue of linear modal superposition in a vicinity of a fixed point. Constructing the ISF from data has shown both initial promise and challenges to be addressed.

Embedding the SSM in the observable space

Originally conceived for autonomous systems, the Takens delay embedding theorem38 has been strengthened and generalized to externally forced dynamics32. By these results, the embedding for a d-dimensional compact SSM subset, \({{{{{{{\mathcal{C}}}}}}}}\subset W(E,{{{{{{{\boldsymbol{\Omega }}}}}}}}t;\epsilon )\), in the delay observable space as \({{{{{{{\mathcal{M}}}}}}}}({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\) is guaranteed for almost all choices of the observable s(t) if p > 2(d + l) and some generic assumptions regarding periodic motions on \({{{{{{{\mathcal{M}}}}}}}}({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\) are satisfied37.

Of highest importance in technological applications is the case of time-periodic forcing (ℓ = 1), with frequency \({{{{{{{\boldsymbol{\Omega }}}}}}}}={{\Omega }}\in {\mathbb{R}}\) and period T = 2π/Ω. In this case, the Whitney and Takens embedding theorems can be applied to the associated period-T sampling map (or Poincaré map) \({{{{{{{{\bf{P}}}}}}}}}_{{t}_{0}}:{{\mathbb{R}}}^{n}\to {{\mathbb{R}}}^{n}\) of the system based at time t0. This map is autonomous and has a time-independent SSM that coincides with the d-dimensional SSM, \({{{{{{{\mathcal{M}}}}}}}}({{\Omega }}{t}_{0})\), of the full system (1). In this case, by direct application of the embedding theorems to the discrete dynamical system generated by \({{{{{{{{\bf{P}}}}}}}}}_{{t}_{0}}\), the typically sufficient embedding dimension estimate is improved to p > 2d for Whitney’s and Takens’s theorem.

Technically speaking, the available data will never be exactly on an SSM, as these embedding theorems assume. By the smoothness of the embeddings, however, points close enough to the SSM in the phase space will be close to \({{{{{{{\mathcal{M}}}}}}}}({{{{{{{\boldsymbol{\Omega }}}}}}}}t)\) in the observable space under the embeddings. Moreover, as slow SSMs attract nearby trajectories exponentially, the distance of observable data from the embedded slow SSM will shrink exponentially fast. Therefore, even under uncorrelated noise in the measurements, mean-squared estimators are suitable for learning slow SSMs from data in the observable space, as we illustrate in the Supplementary Information.

After a possible coordinate shift, the trivial fixed point of the autonomous limit of system (1) will be mapped into the y = 0 origin of the observable space. To find an embedded, d-dimensional SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\in {{\mathbb{R}}}^{p}\), attached to this origin for ϵ = 0, we focus on observable domains in which \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) is a graph over its tangent space \({T}_{{{{{{{{\bf{0}}}}}}}}}{{{{{{{{\mathcal{M}}}}}}}}}_{0}\) at the origin y = 0. Such domains always exist and are generally large enough to capture non-linearizable dynamics in most applications (but see below). Note that \({T}_{{{{{{{{\bf{0}}}}}}}}}{{{{{{{{\mathcal{M}}}}}}}}}_{0}\) coincides with the image of the spectral subspace E in the observable space.

To learn such a graph-style parametrization for \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) from data, we define a matrix \({{{{{{{{\bf{U}}}}}}}}}_{1}\in {{\mathbb{R}}}^{n\times d}\) with columns that are orthonormal vectors spanning the yet unknown \({T}_{{{{{{{{\bf{0}}}}}}}}}{{{{{{{{\mathcal{M}}}}}}}}}_{0}\). The reduced coordinates \({{{{{{{\boldsymbol{\eta }}}}}}}}\in {{\mathbb{R}}}^{d}\) for a point \({{{{{{{\bf{y}}}}}}}}\in {{{{{{{{\mathcal{M}}}}}}}}}_{0}\) are then defined as the orthogonal projection \({{{{{{{\boldsymbol{\eta }}}}}}}}={{{{{{{{\bf{U}}}}}}}}}_{1}^{T}{{{{{{{\bf{y}}}}}}}}\). We week a Taylor-expansion for \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) near the η = 0 origin, denoting by η2:M the family of all monomials of d variables from degree 2 to M. For example, if d = 2 and M = 3, then \({{{{{{{{\boldsymbol{\eta }}}}}}}}}^{2:3}={({\eta }_{1}^{2},{\eta }_{1}{\eta }_{2},{\eta }_{2}^{2},{\eta }_{1}^{3},{\eta }_{1}^{2}{\eta }_{2},{\eta }_{1}{\eta }_{2}^{2},{\eta }_{2}^{3})}^{T}\). As a graph over \({T}_{{{{{{{{\bf{0}}}}}}}}}{{{{{{{{\mathcal{M}}}}}}}}}_{0}\), the manifold \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) is approximated as y = V1η + Vη2:M, where the matrices V1 and V contain coefficients for the d-variate linear and nonlinear monomials, respectively. Learning \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) from a data set of P observations y1, …, yP then amounts to finding the \(({{{{{{{{\bf{U}}}}}}}}}_{1}^{* },{{{{{{{{\bf{V}}}}}}}}}_{1}^{* },{{{{{{{{\bf{V}}}}}}}}}^{* })\) matrices that minimize the mean-square reconstruction error along the training data:

$$\begin{array}{ll}({{{{{{{{\bf{U}}}}}}}}}_{1}^{* },{{{{{{{{\bf{V}}}}}}}}}_{1}^{* },{{{{{{{{\bf{V}}}}}}}}}^{* })=&\arg \mathop{\min }\limits_{{{{{{{{{\bf{U}}}}}}}}}_{1},{{{{{{{{\bf{V}}}}}}}}}_{1},{{{{{{{\bf{V}}}}}}}}}\mathop{\sum }\limits_{j=1}^{P}\parallel {{{{{{{{\bf{y}}}}}}}}}_{j}-{{{{{{{{\bf{V}}}}}}}}}_{1}{{{{{{{{\bf{U}}}}}}}}}_{1}^{T}{{{{{{{{\bf{y}}}}}}}}}_{j}-{{{{{{{\bf{V}}}}}}}}{({{{{{{{{\bf{U}}}}}}}}}_{1}^{T}{{{{{{{{\bf{y}}}}}}}}}_{j})}^{2:M}{\parallel }^{2},\\ &{{{{{{{{\bf{U}}}}}}}}}_{1}^{T}{{{{{{{{\bf{U}}}}}}}}}_{1}={{{{{{{\bf{I}}}}}}}}.\end{array}$$

(18)

The simplest solution to this problem is U1 = V1 with the additional constraint \({{{{{{{{\bf{V}}}}}}}}}_{1}^{T}{{{{{{{\bf{V}}}}}}}}={{{{{{{\bf{0}}}}}}}}\), which represents a basic nonlinear extension of the principal component analysis67.

The above graph-style parametrization of the SSM breaks down for larger y values if \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) develops a fold over \({T}_{0}{{{{{{{{\mathcal{M}}}}}}}}}_{0}\). That creates an issue for model reduction if a nontrivial steady state on \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) falls outside the fold, as the limit cycle does in our vortex shedding example. In that case, alternative parametrization methods for \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\) can be used to enhance the domain of the SSM-reduced model. These methods include selecting the columns of U1 to be the leading POD modes of the nontrivial steady state, or enlarging the embedding space with (further) delayed observations. In these cases, the columns of V1 are still orthonormal vectors spanning \({T}_{{{{{{{{\bf{0}}}}}}}}}{{{{{{{{\mathcal{M}}}}}}}}}_{0}.\)

In both panels (c) of Figs. 4, 6, the SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), is nearly flat in the delay-embedding space. This turns out to be a universal property of delay embedding for small delays and low embedding dimensions (see the Supplementary Information).

For ϵ > 0 small (i.e., for moderate forcing), the autonomous SSM, \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), already captures the bulk nonlinear behavior of system (1). Indeed, for this forcing range, the reduced dynamics on the corresponding SSM can simply be computed as an additive perturbation of the autonomous dynamics on \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\)47,68,69 (see section “Predicting forced dynamics from unforced data”).

SSM dynamics via extended normal forms

For an autonomous SSM \({{{{{{{{\mathcal{M}}}}}}}}}_{0}\), the reduced dynamics is governed by a vector field

$$\dot{{{{{{{{\boldsymbol{\eta }}}}}}}}}={{{{{{{\bf{r}}}}}}}}({{{{{{{\boldsymbol{\eta }}}}}}}})$$

(19)

with a flow map \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{r}}}}}}}}}^{t}({{{{{{{\boldsymbol{\eta }}}}}}}})\). We can generically assume that the Jacobian Dr(0) is semisimple, i.e., Dr(0)B = BΛ, where \({{{{{{{\boldsymbol{\Lambda }}}}}}}}\in {{\mathbb{C}}}^{d\times d}\) is a diagonal matrix containing the eigenvalues of Dr(0). Classic normal form theory would seek to simplify the reduced dynamics (19) in a vicinity of η = 0 via a nonlinear change of coordinates, η = h(z), so that the transformed vector field \(\dot{{{{{{{{\bf{z}}}}}}}}}={{{{{{{\bf{n}}}}}}}}({{{{{{{\bf{z}}}}}}}})\) with flow map \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{n}}}}}}}}}^{t}({{{{{{{\bf{z}}}}}}}})\) has a diagonal linear part and has as few nonlinear terms in its Taylor expansion as possible. In our present setting, the origin is assumed hyperbolic, in which case the classic normal form is simply \(\dot{{{{\bf{z}}}}}={{\mathbf{\Lambda}}}{{{\bf{z}}}}\) under appropriate non-resonance conditions that generically hold40. The corresponding normal form transformation h(z), however, is only valid on a small enough domain in which the dynamics is linearizable.

To capture non-linearizable behavior, we employ extended normal forms motivated by those used to unfold bifurcations33. In this approach, we construct normal forms that do not remove those polynomial terms from (19) whose removal would result in small denominators in the Taylor coefficients h(z) and hence decrease its domain of convergence. Instead, we seek a normal form for (19) of the form

$$\begin{array}{l}{{{{{{{\bf{n}}}}}}}}({{{{{{{\bf{z}}}}}}}};{{{{{{{\bf{N}}}}}}}})={{{{{{{\boldsymbol{\Lambda }}}}}}}}{{{{{{{\bf{z}}}}}}}}+{{{{{{{\bf{N}}}}}}}}{{{{{{{{\bf{z}}}}}}}}}^{2:N},\\ {{{{{{{\bf{h}}}}}}}}({{{{{{{\bf{z}}}}}}}};{{{{{{{\bf{H}}}}}}}})={{{{{{{\bf{B}}}}}}}}({{{{{{{\bf{z}}}}}}}}+{{{{{{{\bf{H}}}}}}}}{{{{{{{{\bf{z}}}}}}}}}^{2:N}),\,\,\,\,\,\,{{{{{{{{\bf{h}}}}}}}}}^{-1}({{{{{{{\boldsymbol{\eta }}}}}}}};{{{{{{{{\bf{H}}}}}}}}}_{{{{{{{{\boldsymbol{\star }}}}}}}}})={{{{{{{{\bf{B}}}}}}}}}^{-1}{{{{{{{\boldsymbol{\eta }}}}}}}}+{{{{{{{{\bf{H}}}}}}}}}_{\star }{({{{{{{{{\bf{B}}}}}}}}}^{-1}{{{{{{{\boldsymbol{\eta }}}}}}}})}^{2:N},\end{array}$$

(20)

where the matrices N, H and H⋆ contain the coefficients for the appropriate d-variate monomials. To identify near-resonances, we let S2:N be the matrix of integers whose columns are the powers of the d-variate monomials from order 2 to N. We then define a matrix Δ2:N containing all relevant integer linear combinations of eigenvalues as follows:

$${({{{{{{{{\boldsymbol{\Delta }}}}}}}}}^{2:N})}_{j,k}={({{{{{{{\rm{Im}}}}}}}}{{{{{{{\boldsymbol{\Lambda }}}}}}}})}_{j,j}-\mathop{\sum }\limits_{s=1}^{d}{({{{{{{{\rm{Im}}}}}}}}{{{{{{{\boldsymbol{\Lambda }}}}}}}})}_{s,s}{({{{{{{{{\bf{S}}}}}}}}}^{2:N})}_{s,k}.$$

(21)

Following the approach used in universal unfolding principles41, we collect in a set S the row and column indices of the entries of Δ2:N for which near-resonances occur, i.e., for which the corresponding entry of Δ2:N is smaller in norm than a small, preselected threshold. (The default threshold is 10−8 in SSMLearn.) The entries of H and H⋆ with indices contained in S are then set to zero but the corresponding monomial terms are retained in n(z; N). Conversely, coefficients of non-near-resonant entries of H and H⋆ are selected in a way so that the corresponding non–near-resonant monomials vanish from the normal form n(z; N). As a result, the matrix N is sparse, containing only the coefficients of essential, near-resonant monomials.

For example, if d = 2, N = 3 and the eigenvalues of Dr(0) form a complex pair λ = α0 ± iω0 with \({\omega }_{0}={{{{{{{\mathcal{O}}}}}}}}(1)\), then we have

$${{{{{{{{\bf{S}}}}}}}}}^{2:N} =\left[\begin{array}{lllllll}2&1&0&3&2&1&0\\ 0&1&2&0&1&2&3\end{array}\right],\,\,{{{{{{{{\boldsymbol{\Delta }}}}}}}}}^{2:N}\\ =\left[\begin{array}{lllllll}-{\omega }_{0}&{\omega }_{0}&3{\omega }_{0}&-2{\omega }_{0}&0&2{\omega }_{0}&4{\omega }_{0}\\ -3{\omega }_{0}&-{\omega }_{0}&{\omega }_{0}&-4{\omega }_{0}&-2{\omega }_{0}&0&2{\omega }_{0}\end{array}\right].$$

(22)

Only two elements of Δ2:N are (near-) zero, and hence the reduced dynamics in extended normal form will require learning the following coefficients:

$${{{{{{{{\bf{H}}}}}}}}}_{\star } =\left[\begin{array}{lllllll}{h}_{20}&{h}_{11}&{h}_{02}&{h}_{30}&0&{h}_{12}&{h}_{03}\\ {\bar{h}}_{02}&{\bar{h}}_{11}&{\bar{h}}_{20}&{\bar{h}}_{03}&{\bar{h}}_{12}&0&{\bar{h}}_{30}\end{array}\right],\,\,\,\,\\ {{{{{{{\bf{N}}}}}}}} =\left[\begin{array}{lllllll}0&0&0&0&{h}_{21}&0&0\\ 0&0&0&0&0&{\bar{h}}_{21}&0\end{array}\right].$$

(23)

The corresponding cubic polar form (5) is then obtained from the relations z = (ρeiθ, ρe−iθ) and h21 = β + iγ.

For a data-driven construction of the extended normal form (20), we first obtain an estimate for the Jacobian Dr(0) from linear regression. This determines the matrix B and the types of monomials arising in h−1 and n. Next, we note that the flow map \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{r}}}}}}}}}^{t}\) of the SSM-reduced dynamics and the flow map \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{n}}}}}}}}}^{t}\) of the extended normal form are connected through the conjugacy relationship \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{n}}}}}}}}}^{t}={{{{{{{{\bf{h}}}}}}}}}^{-1}\circ {{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{r}}}}}}}}}^{t}\circ {{{{{{{\bf{h}}}}}}}}\). We find the nonzero complex coefficients of h−1 and n by minimizing the error in this exact conjugacy over the available P data points, represented in the η coordinates. Specifically, we determine the nonzero elements of H⋆ and N as

$$({{{{{\bf{H}}}}}}_{\star }^{* },{{{{{{\bf{N}}}}}}}^{* })= \arg {\mathop{\min }\limits_{{{{{{\mathbf{H}}}}}}_{\star }}},{{{{{\mathbf{N}}}}}}{\mathop{\sum }\limits_{j=1}^{P}}\Vert \frac{d}{dt}{{{{{{\bf{h}}}}}}}^{-1}({{{{{{\boldsymbol{\eta }}}}}}}_{j};{{{{{{\mathbf{H}}}}}}_{\star }})-{{{{{{{\bf{n}}}}}}}}({{{{{{{{\bf{h}}}}}}}}}^{-1}({{{{{{{{\boldsymbol{\eta }}}}}}}}}_{j};{{{{{{\bf{H}}}}}}_{\star }});{{{{{{{\bf{N}}}}}}}}){\Vert}^{2},\\ {({{{{{{{\bf{N}}}}}}}})}_{s,k}=0,\,\,\forall (s,k)\in S;{({{{{{{\mathbf{H}}}}}}_{\star }})}_{s,k}=0,\forall (s,k)\,\notin\, S.$$

(24)

Once h−1 is known, we obtain the coefficients H of h via regression.

As initial condition for the minimization problem (24), we set all unknown coefficients to zero. This initial guess assumes linear dynamics, which the minimization corrects as needed. We can compute the time derivative in (24) reliably using finite differences, provided that the sampling time Δt of the trajectory data is small compared to the fastest timescale of the SSM dynamics. For larger sampling times, one should use the discrete formulation of SSM theory, as discussed in section “Existence of SSMs” and29. In that formulation, the conjugacy error must be formulated for the 1-step prediction error of the normal form flow map \({{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{n}}}}}}}}}^{{{\Delta }}t}({{{{{{{\bf{z}}}}}}}})\). The matrix defined in eq. (21) also carries over to the discrete time setting, with Λ defined as the diagonal matrix of the logarithms of the eigenvalues of \(D{{{{{{{{\boldsymbol{\varphi }}}}}}}}}_{{{{{{{{\bf{r}}}}}}}}}^{{{\Delta }}t}({{{{{{{\boldsymbol{0}}}}}}}})\).

Prediction of forced response from unforced training data

Forced SSMs continue to be embedded in our observable space, provided that we also include the phase of the forcing among our observables32. (In the simplest case of periodic forcing, this inclusion is not necessary, as we pointed out Section “Embedding SSMs via generic observables”). The quasiperiodic SSM-reduced normal form of system (1) in the observable space takes the general form

$$\begin{array}{c}{\dot{\rho }}_{j}={\alpha }_{j}({{{{{{{\boldsymbol{\rho }}}}}}}},{{{{{{{\boldsymbol{\theta }}}}}}}}){\rho }_{j}-\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}\in {K}_{j}^{\pm }}{f}_{j,{{{{{{{\bf{k}}}}}}}}}\sin \left(\langle {{{{{{{\bf{k}}}}}}}},{{{{{{{\boldsymbol{\Omega }}}}}}}}\rangle t+{\phi }_{j,{{{{{{{\bf{k}}}}}}}}}\mp {\theta }_{j}\right),\\ \!\!\!\!\!\!{\dot{\theta }}_{j}={\omega }_{j}({{{{{{{\boldsymbol{\rho }}}}}}}},{{{{{{{\boldsymbol{\theta }}}}}}}})+\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}\in {K}_{j}^{\pm }}\frac{{f}_{j,{{{{{{{\bf{k}}}}}}}}}}{{\rho }_{j}}\cos \left(\langle {{{{{{{\bf{k}}}}}}}},{{{{{{{\boldsymbol{\Omega }}}}}}}}\rangle t+{\phi }_{j,{{{{{{{\bf{k}}}}}}}}}\mp {\theta }_{j}\right),\end{array}\,\,\,\,\,\,\,j=1,2,...m,\,\,\,\,\,\,{{{{{{{\bf{k}}}}}}}}\in {{\mathbb{Z}}}^{\ell },\,\,\,\,\,\,{{{{{{{\boldsymbol{\Omega }}}}}}}}\in {{\mathbb{R}}}_{+}^{\ell },$$

(25)

where the terms fj,k and ϕj,k are the forcing amplitudes and phases for each mode of the SSM and for each forcing harmonic 〈k, Ω〉, while \({K}_{j}^{\pm }\) are the set containing the indexes k of the resonant forcing frequencies for mode j (see the Supplementary Information). The normal form (25) will capture non-linearizable dynamics arising from resonant interactions between the eigenfrequencies of the spectral subspace E (which may also contain internal resonances) and the external forcing frequencies in Ω. One can use numerical continuation70 to find nontrivial co-existing steady states (such as periodic orbits and invariant tori) in eq. (25) under varying forcing amplitudes and forcing frequencies.

To predict forced response from the SSM-based model trained on unforced data, the forcing amplitude f relevant for eq. (7) in the observable space needs to be related to the forcing amplitude \(\left|\epsilon {{{{{{{{\bf{f}}}}}}}}}_{1}\right|\) relevant for system (1) in the physical phase space. This involves (1) employing a single forcing amplitude-frequency pair \(\left(\left|\epsilon {{{{{{{{\bf{f}}}}}}}}}_{1}\right|,{{\Omega }}\right)\) in the experiment (2) measuring the periodic observable response y(t) (3) computing the corresponding normalized reduced and normalized response amplitude ρ0 (4) substituting ρ0 into the first formula in (8) and (5) solving for f in closed form. This f can then be used to make a prediction for the full FRC and response phase via (8) in the experiment for arbitrary Ω forcing frequencies. The predicted FRC may have several connected components, including isolated responses (isolas) that are notoriously difficult to detect by numerical or experimental continuation68.

Summary of the algorithm

The data-driven model reduction method used in this paper is available in the open-source MATLAB® package SSMLearn. User input is the measured trajectory data of the autonomous dynamical system (ϵ = 0), the SSM dimension d, the polynomial orders or approximation (M, N) for the SSM and for the extended normal form, as well as the type of the dynamical system (discrete or continuous). If the number of observables is not sufficient for manifold embedding, the data is automatically augmented with delays to reach the minimum embedding dimension p = 2d + 1. If the manifold learning returns poor results (due to, e.g., insufficient closeness of the data to the SSM), then the starting value of p can be increased until a good embedding is found. Then, the algorithm learns the SSM geometry in observable space and, after unsupervised detection of the required normal form, identifies the extended normal form of the reduced dynamics. The level of accuracy can be increased with larger polynomial orders, keeping in mind that excessive orders may lead to overfitting.

SSMLearn also offers all the tools we have used in this paper to analyze the reduced dynamics and make predictions for forced response from unforced training data. In particular, it contains the MATLAB®-based numerical continuation core COCO70. which can compute steady state and help with the design of nonlinear control strategies. In principle, there are no restrictions on the dimensions of the reduced-order model, yet the larger the SSM is, the more computationally expensive the problem becomes.

Qualitative or partial a priori knowledge of the linearized dynamics (e.g., some linearized modes and frequencies) helps in finding good initial conditions for trajectories to be used in SSMLearn. For example, the resonance decay method56 (which we exploited in our sloshing example), targets a specific 2-dimensional, stable SSM in laboratory experiments. This method consists of empirically isolating a resonant periodic motion on the SSM based on its locally maximal amplitude response under a forcing frequency sweep. Discontinuing the forcing will then generate transient decay towards the equilibrium in a close proximity of the SSM. For noisy data, filtering or dimensionality reduction can efficiently de-noise the data67, provided that the polynomial orders used for the description of the SSM and its reduced dynamics are not excessively large (see the Supplementary Information). For higher-dimensional SSMs, it is desirable to collect diverse trajectories to avoid bias towards specific motions. Good practice requires splitting the data sets into training, testing and validation parts.

Algorithm 1

SSMLearn

Input parameters: SSM dimension d, polynomial approximation orders (M, N), selection among discrete or continuous-time dynamics

Input data: measured unforced trajectories

Output: SSM geometry, extended normal form of reduced dynamics, predictions for forced response.

1 Embed data in a suitable p-dimensional observable space with p > 2d.

2 Identify the manifold parametrization in reduced coordinates.

3 Estimate the normalized reduced dynamics after an automated identification of the required type of extended normal form.

4 Run analytics and prediction of forced response on the SSM-reduced and normalized model.