Learning Stable Deep Dynamics Models for Partially Observed or Delayed Dynamical Systems

Abstract:
Learning how complex dynamical systems evolve over time is a key challenge in system identification. For safety critical systems, it is often crucial that the learned model is guaranteed to converge to some equilibrium point. To this end, neural ODEs regularized with neural Lyapunov functions are a promising approach when states are fully observed. For practical applications however, partial observations are the norm. As we will demonstrate, initialization of unobserved augmented states can become a key problem for neural ODEs. To alleviate this issue, we propose to augment the system’s state with its history. Inspired by state augmentation in discrete-time systems, we thus obtain neural delay differential equations. Based on classical time delay stability analysis, we then show how to ensure stability of the learned models, and theoretically analyze our approach. Our experiments demonstrate its applicability to stable system identification of partially observed systems and learning a stabilizing feedback policy in delayed feedback control.

heading: 1 Introduction
text: In this paper, we address the task of learning stable, partially observed, continuous-time dynamical systems from data. More specifically, given access to a data set {(t0, yl0), . . . , (tN , ylN )}Ll=1 of noisy, partial observations collected along L trajectories of an unknown, stable dynamical system, ż(t) = g (z(t)) , z(t) ∈ Rm x(t) = h(z(t)) , x(t) ∈ Rn ,m ≥ n, yi = x(ti) + i , i i.i.d.∼ N (0, σ2), (1)
we would like to learn a model for the dynamics of x(t). Moreover, it should be ensured that the model remains stable (we will be concerned with exponential convergence to 0) on unseen trajectories.
Learning such systems in a data-driven way is a key challenge in many disciplines, including robotics [Wensing et al., 2017], continuous-time optimal control [Esposito, 2009] or system biology [Brunton et al., 2016]. One powerful continuous-time approach to non-linear system identification are deep Neural ODEs (NODE), as presented by Chen et al. [2018]. Since neural networks are very expressive, they can be deployed in a variety of applications [Rackauckas et al., 2020]. However, because of that expressiveness, little is known about their system theoretical properties after training. Thus, there has been growing interest in regularizing such dynamics models to ensure favorable properties. In the context of ensuring stability of the learned dynamics, Kolter and Manek [2019] propose to jointly learn a dynamics model and a neural network Lyapunov function, that guarantees global stability via a projection method. Neural network Lyapunov functions have previously been employed by Richards et al. [2018] to estimate the safe region of a fixed feedback policy and by Chang et al. [2019]
∗Correspondence to andreas.schlaginhaufen@epfl.ch, wenkph@ethz.ch.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
to learn a stabilizing feedback policy for given dynamics. Moreover, Boffi et al. [2020] prove that neural Lyapunov functions can also be learned efficiently from data collected along trajectories.
Thus far, all of these approaches are working directly with a standard ODE dynamics model. If the system’s states are fully observed and the system is Markovian, this can be a valid choice. However, in many practical settings, partial observations and non-Markovian effects like hysteresis or delays are the norm. To address the limited expressivity of neural ODEs in the classification setting, Dupont et al. [2019] introduce Augmented Neural ODEs (ANODE). Here, a standard neural ODE is augmented with unobserved states, to extend the family of functions the model is able to capture. While Dupont et al. [2019] demonstrate that initializing the unobserved states at 0 is sufficient for the classification case, this is certainly not true when deploying ANODE as a dynamical system. In fact, our experiments in Section 4 demonstrate that learning this initial condition is a key problem in practice.
Inspired by state-augmentation methods in the time-discrete case, we thus propose to capture partial observability and non-Markovian effects via Neural Delay Differential Equations (NDDE). NDDEs were very recently proposed in the context of classification by Zhu et al. [2021] and in the context of closure models for partial differential equations by Gupta and Lermusiaux [2020]. While NDDEs offer an elegant solution to avoid the Markovianity of neural ODEs, again little can be said about stability outside of the training set. In fact, our experiments in Section 4 show that in a sparse observation and high noise setting, a NDDE model that is stable on training trajectories may become unstable along new unseen trajectories. We therefore extend the ideas of neural network Lyapunov functions, originally developed for stability analysis of non-linear ODEs, to time-delay systems and introduce a Lyapunov-like regularization term to stabilize the NDDE. In contrast to ODEs, NDDEs have an infinite-dimensional state space, which requires careful discretization schemes we introduce in this work. We then showcase the applicability of the proposed framework for the stabilization of the NDDE model and for the task of learning a stabilizing feedback policy in delayed feedback control of known open loop dynamics.
In summary, we demonstrate the applicability of NDDEs to the case of modeling a partially observed dynamical system. We then leverage classical approaches for stability analysis in the context of delayed systems to develop a novel, Lyapunov-like regularization term to stabilize NDDEs. Furthermore, we provide theoretical guarantees and code for our implementation.2

heading: 2 Model and background
text: The main model of this paper, NDDEs, mathematically belongs to the class of time-delay systems that come with some additional difficulties compared to ODEs, both on the theoretical as well as the numerical side. Thus, we first recall some preliminaries and notation on time-delay systems. Then we continue with the model architecture and stability of time-delay systems.

heading: 2.1 Time-delay systems
text: Suppose r > 0 and consider the infinite-dimensional state space Cr := C([−r, 0],Rn) of continuous mappings from the interval [−r, 0] to Rn. Throughout this paper, we endow Rn with the Euclidean norm ‖·‖2 and Cr with the supremum norm ‖φ‖r = sups∈[−r,0] ‖φ(s)‖2 for φ ∈ Cr. Further on, along a trajectory x ∈ C([−r, tf ],Rn) we make use of the notation xt(·) := x(t+ ·) ∈ Cr to denote the infinite-dimensional state at time t ∈ [0, tf ]. A subset B ⊆ Cr is referred to as invariant if γ+(B) = B, where γ+(B) := {xt(ψ) ∈ Cr : ψ ∈ B, t ≥ 0} denotes its positive orbit. For a locally Lipschitz function f : Cr → Rn, an autonomous time-delay system is defined by the family of initial value problems: {
ẋ(t) = f(xt) x(s) = ψ(s), s ∈ [−r, 0], ψ ∈ Cr. (2)
In contrast to ODEs, the dynamics are given by a Functional Differential Equation (FDE) and the initial condition by a function ψ ∈ Cr. As we are interested in autonomous dynamics, we will always set the initial time to zero. Furthermore, we denote by x(ψ)(t) ∈ Rn the solution and by xt(ψ) ∈ Cr the history state at time t, starting from the initial history ψ. As discussed in Section 2.2, we will mainly focus on the important special case of retarded delay differential equations with commensurate delays
ẋ(t) = f(x(t), x(t− τ), ..., x(t−Kτ)) = f(xτ−K(t)). (3)
2Code is available at: https://github.com/andrschl/stable-ndde
Since some discretization is always necessary for computational tracktability, this implicitly includes a numerical approximation of general FDEs if the number of delays is chosen sufficiently large. It holds that r = Kτ , where for convenience we introduced the short-hand notation xτ−K(t) := (x(t), x(t− τ), . . . , x(t−Kτ)). Note that while the instantaneous change (i.e., the vector field) in (3) depends only on a discrete set of observations xτ−K(t), the initial history ψ ∈ Cr has to be given on the entire interval [−r, 0] in order to have well-defined dynamics for t ≥ τ . As a consequence, in practice we need an interpolation of the initial history, and for numerical integration a specific DDE solver based on the method of steps [Alfredo Bellen, 2013] is required. Apart from this, existence and uniqueness of solutions to (2) and (3) follow in a similar fashion as for ODEs [Hale and Lunel, 1993, Diekmann, 1995].

heading: 2.2 Neural Delay Differential Equations
text: Model architecture The motivation of the model architecture is the following: We look for a general method to learn continuous non-Markovian time series which occur, for example, in partially observed dynamical systems. As already mentioned before, the temporal evolution of an ODE is uniquely determined by its current state, which makes NODEs inherently Markovian. Instead of augmenting NODEs with additional states, our approach is inspired by neural network based system identification of discrete-time dynamical systems: the latter copes with non-Markovian effects by augmenting the state space with past observations (i.e., literally memory states) in order to lift the problem back into a Markovian setting (see e.g. [Chen et al., 1990]). A continuous-time analog leads us to a FDE ẋ(t) = f(xt) where the current change is depending on the history xt(s) := x(t + s), s ∈ [−r, 0] up to some maximal delay r > 0. Since a neural network cannot represent a general non-linear functional f , we discretize the infinite-dimensional memory state xt as in Equation (3). This leads us to the NDDE model
ẋ(t) = fNNθ (x(t), x(t− τ), . . . , x(t−Kτ)) = fNNθ ( xτ−K(t) ) , (4)
which is illustrated in Figure 1. Here, fNNθ is a feedforward neural network and K the number of delays.
Predictions Given an initial history ψ ∈ Cr we integrate equation (4) to get the prediction at time t
x̂(t) = DDESolve(ψ, fNNθ , t0, t). (5)
However, in practice observations are subject to noise and cannot be sampled at an infinite rate. Hence, we need to approximate ψ by a smoothed interpolation. For this purpose we employ Gaussian Process (GP) regression. Given a set {(t0, y0), . . . , (tNhist , yNhist} ofNhist observations along the initial history, we fit for each scalar initial history component a zero-mean GP. As a kernel, we choose the Radial Basis Function (RBF) kernel
kγ,σk(t, t ′) = σ2k exp
( −|t− t
′| 2l2
) (6)
with length-scale l and kernel variance σ2k. This choice worked well in our experiments. Nevertheless, it is not crucial and other sufficiently smooth kernels such as Matérn 3/2 or Matérn 5/2 may be appropriate as well [Rasmussen and Williams, 2005]. For the smoothed interpolation of the initial history we are then using the posterior mean function,
ψ(t)i = ktT ( KTT + σ 2I )−1 Yi, 1 ≤ i ≤ n, (7)
where Yi = ((y0)i, . . . , (yNhist)i), T = (t0, . . . , tNhist), ktT = (k(t, t0), . . . , k(t, tNhist), and KTT = (k(tj , tk)) Nhist j,k=1. The kernel hyperparameters l, σ 2 k as well as the observation noise variance σ
2 are estimated from data by marginal likelihood maximization.
Training For training we proceed similar to NODEs and minimize the least squares loss
J = N∑ i=0 ‖yi − x̂(ti)‖22 (8)
along trajectories. While it is possible to utilize an interpolated continuous adjoint sensitivity method for calculating the loss gradients, differentiation through DDE solvers turned out to be significantly more efficient in our experiments. As discussed by Calver and Enright [2016], one reason is that jump discontinuities need to be accounted for that are later propagated in higher order derivatives along the solution of the adjoint state, and the DDE solver needs to be restarted accordingly. We therefore refrain from going into further details about adjoint methods and simply make use of the differentiable DDE solvers provided by Rackauckas and Nie [2017].
Approximation capabilities As opposed to neural ODEs, we are no longer learning an ordinary differential equation, but a retarded-type delay differential equation with constant delays. An interesting question is under which conditions a NDDE can model the time series corresponding to the partial observations h(z(t)) of the ODE system (1). As discussed in Appendix B, a sufficient condition for this is that the delay coordinate map,
E : Rm → R(K+1)n, z(t) 7→ xτ−K(t) = (h(z(t)), h(z(t− τ)), . . . , h(z(t−Kτ))) , (9)
is one-to-one. For dynamical systems confined on periodic or chaotic attractors, the delay embedding theorem by Takens [1981] indeed shows that this holds true for large enough K (for more details see Appendix B and the references therein). Although we do not assume that the system (1) is confined to such an attractor, our experimental results in Section 4 demonstrate approximation power and generalization capabilities of NDDEs when applied to dissipative systems.

heading: 2.3 Stability of time-delay systems
text: We discuss stability analysis for the general class of time delay-systems (2). We assume that the origin is an equilibrium, f(0) = 0, and slightly adjust the definition of exponential stability with respect to this equilibrium point as provided by Fridman [2014] to our needs:
Definition 1 For a fixed set of initial histories S ⊆ Cr and constants γ,M > 0, we call system (2) (γ,M)-exponentially decaying on S over the time horizon [0, tf ) if
‖x(s)‖2 ≤Me −γ(s−t) ‖xt(ψ)‖r for 0 ≤ t ≤ s < tf , ∀ψ ∈ S. (10)
For some invariant set B ⊆ Cr with 0 ∈ B the time delay system (2) is called (γ,M)-exponentially stable on B if it is (γ,M)-exponentially decaying on B over the time horizon [0,∞).
Here, γ measures the rate of decay and M is an upper bound on the transient overshoot. In cases where we do not care about the specific values of γ and M we simply call the system (2) exponentially stable. Note that if a time-delay system is (γ,M)-exponentially decaying on a set of initial histories S over [0,∞), then it is also (γ,M)-exponentially decaying over [0,∞) on γ+(S). Since γ+(S) is invariant by definition, (γ,M)-exponential decay on S over [0,∞) is equivalent to (γ,M)-exponential stability on γ+(S).
Razumikhin’s method A key method to prove exponential stability of non-linear ODEs are Lyapunov functions [Lyapunov, 1992]. However, directly applying ODE Lyapunov functions to time-delay systems leads to very restrictive results (e.g., a 1-dimensional system would not be allowed to oscillate [Fridman, 2014]). Nevertheless, along the same lines, two approaches geared towards stability analysis of non-linear time-delay systems exist. Whereas the method of Lyapunov-Krasovskii functionals [Krasovskii, 1963] is a natural extension of Lyapunov’s direct method to an infinite-dimensional state space, the idea of Razumikhin-type theorems [Razumikhin, 1956] is to make use of positive-definite Lyapunov functions V : Rn → R+ with finite domains familiar from the ODE case and to relax the decay condition. Namely, a negative derivative of V (x(t)) at time t is required only when we are about to leave the sublevel set V ≤η = {x ∈ Rn : V (x) ≤ η} of η := sups∈[−r,0] V (x(t + s)). The following theorem establishes sufficient conditions for exponential stability.
Theorem 1 ([Efimov and Aleksandrov, 2020]) Assume there exists a differentiable function V : Rn → R+, positive reals c1, c2, α, and a constant q > 1 such that along all trajectories starting in ψ ∈ S ⊆ Cr the following conditions hold for all t ∈ [0, tf ):
(i) c1 ‖x(t)‖22 ≤ V (x(t)) ≤ c2 ‖x(t)‖ 2 2
(ii) V̇ (x(t)) ≤ −αV (x(t)) whenever V (x(t+ s)) ≤ qV (x(t)) ∀s ∈ [−r, 0].
Then system (2) is (γ,M)-exponentially decaying on S over [0, T ) with decay rate γ = min(α, log qr )/2 and M = c2/c1. Moreover, if for an invariant set B ⊆ Cr with 0 ∈ B conditions (i), (ii) hold for all xt ∈ B then the time delay system (2) is exponentially stable on B.
A function V establishing stability of some invariant set B by satisfying conditions (i), (ii) in Theorem 1 is referred to as a Lyapunov-Razumikhin Function (LRF). However, due to the infinite dimension of Cr it is hard to verify the decay condition (ii) on the entire state space Cr. We thus focus on proving (γ,M)-exponential decay along trajectories starting within some fixed set of initial conditions S ⊂ Cr, which, as discussed before, is for an infinite time horizon equivalent to (γ,M)-exponential stability on B = γ+(S). Note that Theorem 1 establishes sufficient, but not necessary conditions for exponential stability. The problem is that often it is not strong enough to check the Razumikhin condition
qV (x(t))− V (x(t+ s)) ≥ 0 (11) only on the interval s ∈ [−r, 0], but we should take into account more of the past observations of V (x(t)). It can therefore be helpful to reinterpret problem (2) as one in the state space CrV with some rV > r and to apply Theorem 1 to that problem. However, in this new – larger – state space, only initial histories of the form
ψ̃(s) = { ψ(s− (r − rV )) , s ∈ [−rV , r − rV ] x(ψ) (s− (r − rV )) , s ∈ [r − rV , 0], with ψ ∈ Cr (12)
need to be considered for stability in Cr [Hale and Lunel, 1993]. Furthermore, Proposition 1, which we prove in Appendix A, shows that also exponential stability follows, albeit at the price of a larger bound on the transient overshoot.
Proposition 1 If f is Lf -Lipschitz, f(0) = 0, and ψ, ψ̃ defined as in (12), then∥∥∥ψ̃∥∥∥ rV ≤ ‖ψ‖r e Lf (rV −r). (13)
Centred around this idea, necessary and sufficient Razumikhin-type conditions for discrete-time delay systems are given by Gielen et al. [2013]. In the following, we therefore treat rV as a hyperparameter that has to be chosen for the respective problem at hand.

heading: 3 Learning stable dynamics
text: We now propose an approach, based on neural LRFs, to enforce stability of a parametric DDE. The key idea is to jointly learn a neural network Lyapunov-Razumikhin function and the dynamics model. Similarly to [Richards et al., 2018, Chang et al., 2019] we propose to enforce stability via the loss function. This is in contrast to Kolter and Manek [2019] who use a projection-based approach to ensure stability in the forward pass. The main reason for this design choice is that a projective approach based on LRFs leads to discontinuities in the forward pass, which are problematic for DDE solvers [Alfredo Bellen, 2013]. Moreover, incorporating the Lyapunov neural network into the forward pass renders the model slow during inference time and a loss function based approach offers the opportunity to actively stabilize an initially unstable system, as we demonstrate in Section 4.
Lyapunov neural network construction Except for the decay condition along solutions (condition (ii) in Theorem 1), an LRF has the same form as an ODE Lyapunov function. We thus employ the same Lyapunov neural network as proposed by Kolter and Manek [2019]. The construction is based on an Input-Convex Neural Network (ICNN) [Amos et al., 2016]. The ICNN x 7→ gNNφ (x) is convex by construction and any convex function can be approximated by such neural networks [Chen et al., 2019]. In order to satisfy the upper bound in condition (i) of Theorem 1, the
activation functions σ of gNNφ are required to additionally have slope no greater than one. To ensure strict convexity, and to make sure that the global minimum lies at x = 0, a final layer
V NNφ (x) = σ(g NN φ (x)− gNNφ (0)) + c ‖x‖ 2 2 (14)
is chosen. Here, c > 0 is a small constant. As for the activation function σ, having a global minimum at x = 0 requires σ(0) = 0. Furthermore, since we want to ensure Lipschitz continuity of the loss derivatives, we use a twice continuously differentiable smoothed ReLU version
σ(x) =  0 , x ≤ 0 x3 d2 − x 4 (2d)3 , 0 ≤ x ≤ d
x− d 2
, x > d.
(15)
This slightly differs from the original σ proposed by Kolter and Manek [2019], since they only needed a once continuously differentiable one. This construction ensures that V NNφ (x) = O(‖x‖ 2 2) as x → 0 and also V NNφ (x) = O(‖x‖ 2 2) as ‖x‖2 → ∞. We can therefore always find constants c1, c2 such that the conditions (i) in Theorem 1 are satisfied. In the next step, we explain how to employ this neural network architecture to learn neural LRFs and at the same time stabilize a parametric delay differential equation of the form (3).
Lyapunov-Razumikhin loss As stated before, V NNφ satisfies condition (i) in Theorem 1 by construction. The relaxed decay condition (ii) however has to be enforced during training. Since it is practically infeasible to check the Razumikhin condition (11) on the continuous interval [−rV , 0], we need some discretization that still allows for stability guarantees. As we will analyze in this section, this is satisfied by the the following loss with discretized Razumikhin condition
`LRF ( φ, θ,xτV−KV (t) ) =
ReLU ( V̇ NNφ,θ ( xτV−KV (t) ) + αV NNφ (x(t)) ) Θ ( qV NNφ (x(t))− max
1≤j≤KV V NNφ (x(t− jτV ))
) .
(16)
Here, Θ(·) denotes the unit step function with Θ(s) = 1 if s ≥ 0 and Θ(s) = 0 otherwise. Furthermore, for notational simplicity we choose τV ≤ τ and such that τ = l · τV for some integer l ∈ N. According to Theorem 1, a zero loss `LRF ( φ, θ,xτV−KV (t) ) along a trajectory of length [0, tf ) implies exponential decay along this trajectory. Moreover, if for a fixed set of initial histories Strain ⊆ Cr the loss (16) is zero along all trajectories starting in Strain and over a time horizon [0,∞), then the delay differential equation is stable on γ+(Strain). However, since we cannot check this for tf =∞, we choose tf large enough to ensure convergence to a sufficiently small region around the origin. Theorem 2 then also establishes exponential decay for trajectories starting not necessarily in – but close enough to – Strain. For its proof we refer to Appendix A.
Theorem 2 If the dynamics are Lf -Lipschitz and the LRF loss (16) is zero along trajectories starting in Strain ⊂ Cr over a time horizon [0, tf ), then the time-delay system is (γ,M)-exponentially decaying on Strain over [0, tf ). Moreover, if for another set of initial histories S ⊃ Strain and some ε > 0, the training set Strain is a δ-covering of S (in the ‖·‖r-norm) with δ = εe−(Lf+γ)tf , then the time delay system is (γ, M̄)-exponentially decaying on S \Bε(0) over the time horizon [0, tf ) and with M̄ = 2M + 1. Here, Bε(0) = {ψ ∈ Cr : ||ψ||r ≤ ε} denotes the ε-ball around the origin.
While for a zero loss exponential decay is guaranteed, the discretization of the Razumikhin condition might be introducing additional conservatism by requiring decay in V NNφ too often. However, Proposition 2 tells us that if the discretized Razumikhin condition holds, τV is small enough, and the current state lies outside of an ε-ball around the origin, then the continuous condition holds for some q̃ > q. Furthermore, q̃ converges quadratically to q as τV → 0. Remembering that the decay rate in Theorem 1 is γ = min(α, log qr )/2, it becomes apparent that by discretization we are requiring a slightly larger rate of decay, which can however be controlled by the choice of τV .
Proposition 2 Let KV and τV be such that rV := KV τV ≥ r and let x(·) be a solution of ẋ(t) = f(xt) passing through xt0 ∈ Cr. Assume ‖xt0‖rV +2r <∞ and ‖xt0‖rV > ε. Furthermore, let f be Lf -Lipschitz and differentiable. Then, if the discretized Razumikhin condition,
V NNφ (x(t0 − kτV )) ≤ qV NNφ (x(t0)) ,∀k ∈ {0, 1, ...,KV }, (17) is satisfied and τV is small enough, then the continuous Razumikhin condition,
V NNφ (x(t0 + s)) ≤ q̃(τV )V NNφ (x(t0)), (18)
holds for any s ∈ [−rV , 0] and some q̃(τV ) with q̃(τV ) = q +O(τ2V ) as τV → 0.
The condition that the current state lies outside some ε-ball around the origin may be replaced by an assumption on the solutions’ decay in ‖·‖rV . However, in practice we are usually satisfied with convergence to some small neighborhood of the origin, since, as already noted, we cannot choose an infinite time horizon and also as due to observation noise and data scarcity our model will always be subject to some modelling errors. We elaborate more on this issue and prove Proposition 2 in Appendix A.
Stabilizing NDDEs In order to stabilize an NDDE on a fixed set of initial conditions Strain ⊆ Cr we minimize the stabilizing loss (16) on a set of data points {x(1), . . . ,x(NLRF)} with x(i) ∈ Rn(KV +1) collected along trajectories starting in Strain. The resulting gradients are added up with those from the NDDE loss (8). While this enables us to stabilize the NDDE on unseen trajectories, we still need an efficient method to generate realistic initial histories. Especially in the setting of partially observed systems we do not know much more about initial histories than that they are contained within a bounded, Lipschitz subset of Cr. However, since for our NDDE model the initial history is given by a GP-mean function, it suffices to stabilize the NDDE for initial histories within the subset
{ψ ∈ Hk|ψ(t) = Nhist∑ i=1 cikl,σk(t, ti)} (19)
of the reproducing kernel Hilbert spaceHk corresponding to kl,σk(t, t′). Furthermore, boundedness of initial histories translates into a bound on the norm of the expansion coefficients ‖(c1, ..., cNhist)‖2 ≤ A and Lipschitz continuity can be accounted for by upper bounding the inverse length-scale 1/l ≤ B and the kernel variance σ2k ≤ C [Rasmussen and Williams, 2005]. To satisfy these constraints, we sample at each training iteration initial histories ψ ∈ Strain as follows: The expansion coefficients (c1, ..., cNhist) are sampled uniformly in an L2-ball, and 1/l, σk on bounded intervals [0, B], [0, C], respectively.
Note, that while another possibility would be to integrate the loss (16) as a continuous regularization term into the NDDE loss, the discontinuities in (16) turn out to be problematic for DDE solvers.
Delayed feedback control The stabilizing loss (16) is essentially applicable to any parametric DDE of the form (3). Equations of this form also occur in delayed feedback control. Assume we want to learn a stabilizing state feedback u(t) = πθ(x(t)) for a known open loop control system ẋ(t) = f(x(t), u(t − τ)) with input delay. In practice, such delays in the feedback loop are often introduced as a consequence of communication latencies and typically cause instability [Krstic, 2009]. The resulting closed loop system is a parametric DDE ẋ(t) = f(x(t), πθ(x(t − τ))), which can, for small enough delays, be stabilized in a data-driven way with our LRF loss (16). If the input delay exceeds some critical value, the system can no longer be stabilized by DDE methods and infinite-dimensional feedback taking into account the inputs history would be required [Krstic, 2009]. Experimental results for delayed feedback stabilization are provided in Section 4.
Choice of hyperparameters Our NDDE model (4) as well as the stabilizing loss (16) involve hyperparameters such as number and magnitude of delays, whose choice we discuss in the following. For our NDDE model, the number of delays K clearly controls the representational capabilities. In general it is sufficient to choose K large enough such that the delay coordinate map (9) is one-to-one. For periodic or chaotic attractors, Takens’ Embedding Theorem 3 (see Appendix B) provides a sufficient lower bound on K to ensure this. Moreover, in our experiments, larger values of K ease training and – perhaps surprisingly – do not hurt generalization performance. Of course, an overly large number of delays leads to long training time per iteration, thus slowing down training again. Except for the first experiment where we directly compare NDDEs to ANODEs, we fix a relatively large number of delays K = 10 throughout the paper. While Takens’ Theorem 3, is besides a periodicity condition, completely agnostic to the choice of the delay parameter τ , various heuristics such as Average Mutual Information or False Nearest Neighbours exist in practice (for an overview see [Wallot and Mønster, 2018]).
With regard to the stabilizing loss (16), Proposition 2 proves that the number of delays KV controls the conservatism we introduce through discretization of the Razumikhin condition. Furthermore, as discussed in (12), the maximal considered delay rV controls the conservatism inherent to Razumikhin’s Theorem 1 itself. Lastly, the parameters α and q are directly related to the rate of decay γ in Theorem 1 via γ = min(α, log q/rV ). We thus choose α ≈ log q/rV . Moreover, a too small choice of KV , rV or too large choice of α, q can be detected via a non-zero LRF loss (16).

heading: 4 Experiments
text: Learning partially observed dynamics We first compare the applicability of Vanilla NDDEs and ANODEs for the task of learning a partially observed harmonic oscillator,
ż(t) = d
dt ( z1(t) z2(t) ) = ( 0 1 −1 0 ) z(t), h(z(t)) = (1 0) z(t). (20)
We train the models over two training trajectories starting from z0,1 = (1, 0) and z0,1 = (0, 2) and with zero observation noise. For ANODEs we compare a model trained with given true augmented initial conditions (IC) against another model where we initialize the augmented states with zero and learn them via the adjoint method. Moreover, for the NDDE we compare a single delay model with K = 1 to a multiple delay model with K = 10. The resulting vector field plots for the ANODE models illustrated in Figures 2a-2c demonstrate that, whereas for true initial conditions the dynamics match the ground truth well, learning the augmented initial conditions turns out to be a key problem. In contrast, for our NDDE model the initialization is conveniently provided by the GP interpolation. This is also reflected in the learning curves in Figure 2d, where we see that both NDDE models yield a significantly lower train loss for fewer iterations compared to the ANODE models. Moreover, the NDDE with K = 10 achieves a better training score. For the rest of the experiments we therefore fix K = 10. For more information about the setup and additional experiments we refer to Appendix C.
Learning stable NDDEs Kolter and Manek [2019] show that in the ODE case, a neural network dynamics model trained on stable data may become unstable for long-term prediction. For NDDEs, we observed this to be a problem in the setting of sparse observations, a high noise level, and generalization over initial histories. In particular, we consider a partially observed damped double compound pendulum, where only the angles of deflection ϕ1 and ϕ2, but not the angular velocities are observed. This is a complex non-linear dynamical system which, for low friction, exhibits chaotic behavior [Shinbrot et al., 1992]. The governing equations are derived in Appendix C.
For observation noise of variance σ = 0.05 and training and test data along 4 trajectories, we compare the generalization performance of a Vanilla NDDE and a NDDE stabilized with LRF regularization. We repeat the training for 20 independent weight initializations and noise realizations. The resulting predictions illustrated in Figures 3a-3b demonstrate that while the median prediction is stable, the upper 0.95 quantile explodes for the Vanilla NDDE. In contrast, the stabilized NDDE remains stable
on all test trajectories. Moreover, whereas the test loss in Figure 3c explodes for the unstable NDDE, the train losses are approximately the same. Thus, the LRF loss guides us to a stable optimum without sacrificing training performance.
Stabilization with delayed feedback control As a first application for learning a stabilizing feedback policy of a known open loop system, we consider a friction-less inverted pendulum with an input delay τ = 0.03. The open loop dynamics are given by(
ẋ1(t) ẋ2(t)
) = ( x2(t) g
l sin(x1(t)) +
1
ml2 u(t− τ)
) . (21)
Here, the states are (x1, x2) = (ϕ, ϕ̇) where ϕ(t) is the angle of deflection with respect to the fully upright position, g indicates the acceleration of gravity, and l and m the length and mass of the pendulum. Furthermore, u(t) is the torque which is applied at the pivot point. The goal is to learn a stabilizing feedback policy
u(t) = π(x(t)) = k1x1(t) + k2x2(t). (22)
Similar to Chang et al. [2019], we initialize the parameters k1, k2 with the values from the Linear Quadratic Regulator (LQR) feedback policy calculated for the linearization of (21). For the training, we continuously generate new initial histories as follows: We sample ODE initial conditions on
a circle of radius π/2, assuming zero control for t < 0. Thus, the dynamics are described by an autonomous ODE along initial histories. As depicted in Figure 4a, the initially unstable state feedback can be stabilized by means of our Razumikhin loss. Furthermore, the speed of decay can be controlled by the choice of the hyperparameters α and q. Moreover, `LRF is zero along new test trajectories indicating that we indeed learned a valid LRF candidate for this set of initial histories.
As a second – more complex – experiment, we consider stabilizing a cartpole with a delayed input force acting on the cart. In contrast to the two-dimensional inverted pendulum, this is a fourdimensional non-linear system. Its states are x(t) = (ϕ(t), ϕ̇(t), ξ(t), ξ̇(t)), where ϕ again denotes the angle of deflection and ξ the position of the cart. For the exact equations, we refer to [Stimac, 1999]. For the control force acting on the cart we assume a delay of τ = 0.05 and aim at finding a stabilizing feedback policy πθ(x(t)). Similarly to the inverted pendulum experiment, Figure 5 shows that minimizing the LRF loss (16) enables us to find a stabilizing feedback policy from an initially unstable LQR feedback. Furthermore, the rate of decay can be controlled by the choice of the hyperparameter α and q.

heading: 5 Conclusion
text: In this paper, we demonstrated that NDDEs are a powerful tool to learn non-Markovian dynamics occuring when observing a partially observed dynamical system. Via state augmentation with the past history we avoid the estimation of unobserved augmented states, which we showed to be a major problem of ANODEs when applied to partially observed systems. Based on classical time-delay stability theory, we then proposed a new regularization term based on a neural network LyapunovRazumikhin function to stabilize NDDEs. We further showed how this approach can be used to learn a stabilizing feedback policy for control systems with input delays. Besides experiments showcasing the applicability of our approach, we also provide code and a theoretical analysis.

heading: Acknowledgments
text: This research was supported by the Max Planck ETH Center for Learning Systems. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement No 815943 as well as from the Swiss National Science Foundation under NCCR Automation, grant agreement 51NF40 180545.

Review summaries:
The paper presents a method to learn the dynamics of stable systems. The authors first exploit Neural Delay Differential Equations to model the partial observability and time delay of the system. They then employ Lyapunov-Razumikhin Function to ensure the stability of the learned dynamics. The authors evaluate their methods by presenting the sampled stable trajectories and loss plots of the learned model against other selected models.

This paper aims to learn stable deep dynamics models for partially observed or delayed dynamical systems. The authors build their method upon Neural Delay Differential Equations (NDDEs) [1], where the main idea is to augment the state with the past history, thus avoiding the estimation of unobserved augmented states. To stabilize the learned dynamics, the authors then propose to combine NDDEs with a Lyapunov-like regularization term.  The authors show that NDDEs are powerful tools to model partially observed dynamical systems where the dynamics are non-Markovian. Furthermore, they demonstrate that the proposed method outperforms Augmented Neural ODEs (ANODE) in partially observed systems and is able to learn a stabilizing feedback policy for control systems with input delays.  [1] Qunxi Zhu, Yao Guo, Wei Lin, "Neural Delay Differential Equations"

This paper extends prior work on stable deep neural network dynamics models to learn stable models with partially observed states, and to learn (simple) control policies with time-delayed feedback. To achieve this, the authors minimize a Lyapunov-Razumikhin loss term when training _neural delay differential equations_. They demonstrate the capabilities of their system on a few tasks.