Neural Contractive Dynamical Systems

We are grateful to the reviewer for their interest in our work (which was excellently summarized in the review) and for valuable feedback.

• “Is training end-to-end? What are the loss functions used for training?”

  • Short answer: Training NCDS is end-to-end, but the latent NCDS is not fully end-to-end.

  • Detailed answer: Training NCDS is end-to-end, but the latent NCDS is not fully end-to-end. In the latter case, we first train the VAE (end-to-end), and then train the latent NCDS using the encoded data. In more detail, the VAE is trained using the standard evidence lower bound (ELBO):

    $\\mathcal{L} _{ELBO} = \\mathbb{E} _{q _{\\mathbf{\\xi }}(\\mathbf{z}|\\mathbf{x})}\\left[\\log(p _{\\mathbf{\\phi}}(\\mathbf{x}|\\mathbf{z}))\\right] - \\mathrm{KL}\\left(q _{\\mathbf{\\xi }}(\\mathbf{z}|\\mathbf{x})||p(\\mathbf{z})\\right) ,$

    where $\mathrm{KL}$ denotes the Kullback-Leibler divergence.

    Thereafter, the Jacobian network is trained according to the loss function $\mathcal{L}_{\text{Jac}}$, which measures the distance between the observed and approximated next state:

    $\\mathcal{L} _{\\text{Jac}} = \\|\\mathbf{z} _{t+1} - (\\mathbf{z} _{t} + \\hat{\\dot{\\mathbf{z}}} _t)\\|^2,$

    where $\\mathbf{z} _{t}$, $\\mathbf{z} _{t+1}$, and $\\hat{\\dot{\\mathbf{z}}}_t$ represent the current and next observed latent states, along with the calculated latent velocity.

  • Paper changes: we have updated the paper and added the loss functions and the algorithms in the Appendix A.4: Implementation details.

Algorithm 1: Neural Contractive Dynamical Systems (NCDS): Training in task space

Data: Demonstrations: $\\tau _n = \\left \\{(\\mathbf{x} _t, \\mathbf{R} _t) \\right \\}, n \\in [1, N], t \\in [1, T _n]$

Result: Learned contractive dynamical system

1. $\\mathbf{r} _{n,t} = \\operatorname{Log}(\\mathbf{R} _{n,t})$; // Obtaining skew-symmetric coefficients
2. $\\mathbf{p} _{n, t} = [\\mathbf{x} _{n, t}, \\mathbf{r} _{n, t}]$; // Create new state vector
3. \\operatorname{argmin}_\\xi \\mathcal{L} _{\\text{ELBO}}(\\xi^*; \mathbf{p} _{n,t}); // Train the VAE
4. $\\mathbf{z} _{n,t} = \\mu _ {\\mathbf{\xi}} ^{\\sim 1} (\\mathbf{p} _{n,t})$; // Encode all states using the trained VAE
5. $\\dot{\mathbf{z}} _{n,t} = \\frac{\\mathbf{z} _{n,t+1} - \\mathbf{z} _{n,t}}{\Delta t}$; // Compute latent finite differences
6. \\operatorname{argmin}_\\theta \mathcal{L} _{\\text{Jac}}(\\theta^*; \\mathbf{p} _{n,t}); // Train the Jacobian network

Algorithm 2: Algorithm Neural Contractive Dynamical Systems (NCDS): Robot Control Scheme

Data: Current state of the robot end-effector at time $t$: $[\\mathbf{x} _t, \\mathbf{R} _t]$

Result: Velocity of the end-effector at current time step $\\dot{\\mathbf{x}} _t$

1. $\\mathbf{r} _{t} = \\operatorname{Log}(\\mathbf{R} _{t})$; // Obtaining skew-symmetric coefficients
2. $\\mathbf{p} _{t} = [\\mathbf{x} _{t}, \\mathbf{r} _{t}]$; // Create new state vector
3. $\\mathbf{z} _t = \\mu _{\\mathbf{\\xi}} ^{\\sim 1} (\\mathbf{p} _t)$; // Compute the latent state
4. $\\hat{\\dot{\\mathbf{z}}} _t = f(\\mathbf{z} _t)$; // Compute the latent velocity
5. $J_{\\mu _{\\mathbf{\\xi}}}(\\mathbf{z} _t) = \\frac{\\partial\\mu _{\\mathbf{\\xi}}}{\\partial \\mathbf{z} _t}$; // Compute the Jacobian of the decoder
6. $\\dot{\\mathbf{x}} _t = J _{\\mu _{\\mathbf{\\xi}}}(\\mathbf{z} _t)\\hat{\\dot{\\mathbf{z}}} _t$; // Compute input space velocity

• “How are the sequential data used for training? Is each trajectory used as a whole as one data point, or cut into segments as multiple data points? Second-order?”

  • During training, data points from all the demonstrations undergo shuffling and are then fed into the network as a batch containing $B$ data points.

• “Multi-modality”

  • We do not quite understand the reviewer’s question here. The term "multimodality" in the context of a dynamical system could mean different things. We are happy to do our best to answer the raised question but request a bit more detail.

• “Limitation of baseline method choices: The baselines compared to are all focused on asymptotic stability guarantees. While these are the most relevant methods for comparison, it would be interesting to see how imitation learning methods without any stability guarantee works on the tasks, to demonstrate the necessity of stability guarantee.”

  • This is a fair point. We opted for the chosen baselines as our work was really driven by a desire to provide stability. We will try to conduct these additional baseline experiments within the discussion period and report back with results later.

• Reviewer’s comment: “Is there a choice of integrator that doesn't require adaptive step-sizes (perhaps a symplectic integrator)?”
  • Detailed answer: Given our limited experience with symplectic integrators, specifically designed for Hamiltonian systems, we find them less applicable to our current ODE-based system. Our research indicates that symplectic integrators excel in preserving certain properties, such as energy conservation, which are integral to Hamiltonian dynamics. Since our system does not exhibit these characteristics, opting for a symplectic integrator might introduce unnecessary complexity. On the contrary, the Runge-Kutta 4th order (RK4) method, while not inherently adaptive in step size, is a fixed-step-size numerical integration technique widely employed for solving ordinary differential equations (ODEs). This method strikes a balance between accuracy and computational efficiency, particularly in ODE systems with straightforward dynamics where stiffness is not a prominent concern. As highlighted in Figure 13 of the paper, RK4 exhibits a certain level of accuracy, although it is not as precise as adaptive step size methods such as the Dopri5 method. This observation underscores the importance of considering the specific characteristics and requirements of our ODE system when selecting an integration method. In cases where our system dynamics are uncomplicated and adaptive step sizes are not imperative, RK4 remains a viable and pragmatic choice.
• Reviewer’s comment: “It would be helpful if the training algorithm were stated explicitly in Section 3.”
  • Short answer: We agree that a training algorithm should be added to the paper.
  • Detailed answer: Here, we provide detailed steps of the training and robot control schemes for NCDS. The sequential steps for training the Variational Autoencoder (VAE) and Jacobian network are listed in Algorithm 1. Simultaneously, the procedural steps for employing NCDS to control a robot are shown in Algorithm 2.
  • Revision: A new subsection titled "Algorithms" has been included in the Appendix A.4, containing the respective algorithms for reference.

Algorithm 1: Neural Contractive Dynamical Systems (NCDS): Training in task space

Data: Demonstrations: $\\tau _n = \\left \\{(\\mathbf{x} _t, \\mathbf{R} _t) \\right \\}, n \\in [1, N], t \\in [1, T _n]$

Result: Learned contractive dynamical system

1. $\\mathbf{r} _{n,t} = \\operatorname{Log}(\\mathbf{R} _{n,t})$; // Obtaining skew-symmetric coefficients
2. $\\mathbf{p} _{n, t} = [\\mathbf{x} _{n, t}, \\mathbf{r} _{n, t}]$; // Create new state vector
3. \\operatorname{argmin}_\\xi \\mathcal{L} _{\\text{ELBO}}(\\xi^*; \mathbf{p} _{n,t}); // Train the VAE
4. $\\mathbf{z} _{n,t} = \\mu _ {\\mathbf{\xi}} ^{\\sim 1} (\\mathbf{p} _{n,t})$; // Encode all states using the trained VAE
5. $\\dot{\mathbf{z}} _{n,t} = \\frac{\\mathbf{z} _{n,t+1} - \\mathbf{z} _{n,t}}{\Delta t}$; // Compute latent finite differences
6. \\operatorname{argmin}_\\theta \mathcal{L} _{\\text{Jac}}(\\theta^*; \\mathbf{p} _{n,t}); // Train the Jacobian network

Algorithm 2: Algorithm Neural Contractive Dynamical Systems (NCDS): Robot Control Scheme

Data: Current state of the robot end-effector at time $t$: $[\\mathbf{x} _t, \\mathbf{R} _t]$

Result: Velocity of the end-effector at current time step $\\dot{\\mathbf{x}} _t$

1. $\\mathbf{r} _{t} = \\operatorname{Log}(\\mathbf{R} _{t})$; // Obtaining skew-symmetric coefficients
2. $\\mathbf{p} _{t} = [\\mathbf{x} _{t}, \\mathbf{r} _{t}]$; // Create new state vector
3. $\\mathbf{z} _t = \\mu _{\\mathbf{\\xi}} ^{\\sim 1} (\\mathbf{p} _t)$; // Compute the latent state
4. $\\hat{\\dot{\\mathbf{z}}} _t = f(\\mathbf{z} _t)$; // Compute the latent velocity
5. $J_{\\mu _{\\mathbf{\\xi}}}(\\mathbf{z} _t) = \\frac{\\partial\\mu _{\\mathbf{\\xi}}}{\\partial \\mathbf{z} _t}$; // Compute the Jacobian of the decoder
6. $\\dot{\\mathbf{x}} _t = J _{\\mu _{\\mathbf{\\xi}}}(\\mathbf{z} _t)\\hat{\\dot{\\mathbf{z}}} _t$; // Compute input space velocity

• Reviewer’s comment: “It would be nice to have the invariance of the contractivity stated formally.”
  • Short answer: We agree that the paper was not precise enough here.
  • Detailed answer: Let $\mu: \mathbb{R}^d \rightarrow \mathbb{R}^D$ be the injective decoder mapping from the low-dimensional latent space to the high-dimensional observation space. Further, let $z_t$ denote a contractive dynamical system in the latent space. Geometrically, the decoder spans a $d$-dimensional manifold embedded in $\mathbb{R}^D$ defined as $\mathcal{M} = f(\mathbb{R}^d)$. By design, the decoder is a diffeomorphism between $\mathbb{R}^d$ and $\mathcal{M}$; this construction is extensively discussed by Brehmer & Cranmer [1], which we build upon. By Theorem 1 in our paper, the above implies that $f(z_t)$ is contractive. The above further implies that the dynamical system $f(z_t)$ is along a low-dimensional manifold $\mathcal{M} \subset \mathbb{R}^D$, i.e. it does not cover all of $\mathbb{R}^D$. Such constructions are common in latent variable models. The control loop assumes that the initial robot configuration is in $\mathcal{M}$ such that it can be written $f(z_0)$. Following the dynamical system moves the robot along the manifold, and the resulting motion follows a contractive system. If the assumption is not satisfied, such that the initial configuration $x_0 \not\in \mathcal{M}$, then we have to rely on the encoder to approximate a projection onto $\mathcal{M}$, i.e. produce $z_0 = \mu^{\sim 1}(x_0)$. Note that the movement from $x_0$ to $f(z_0)$ need not be contractive, such that the distance to the nominal trajectory may increase. This is, however, a finite-time motion, after which the system is contractive.
  • Revisions: In Section 3.2, we have refined and formalized the correlation between injection and contraction.
• Reviewer’s comment: “In the discussion section (page 9), it is mentioned that the choice of integration scheme can siginificantly affect the behaviour of the learned model. This requires further discussion. For instance,”
  • “How significantly does the computation time affect the performance of the model?”
    • Short answer: Our findings show that as long as we keep the latent space low dimensional, the computation of the integration does not affect the overall performance of the system.
    • Detailed answer: We have extensively investigated this question in Section A.6.2: Dimensionality and execution time, and our findings show that as long as we keep the latent space low dimensional, the computation of the integration does not affect the overall performance of the system. It is important to acknowledge that while maintaining a low-dimensional latent space is favorable for real-time applications, the input space dimensionality introduces computational challenges. Specifically, computing the Jacobian of the injective generator at each integration step becomes more resource-intensive. Furthermore, to answer this question in a more extreme scenario we have conducted an experiment on higher-dimensional ($44D$) human motion dataset. The new experiment showcases the performance of injective generators in encoding complex, high-dimensional data into a low-dimensional latent space ($2D$). Certainly, the considerable dimensionality of the data space inevitably impacts the runtime, with each iteration of the control loop in Algorithm 2 requiring 121 milliseconds.
    • Revision: We have added a new subsection in the Appendix A.6.2 titled “Learning human motion”. We also updated Table 2 in the appendix with the new results from human motion experiments execution time.

We appreciate the reviewer for taking the time to provide thoughtful and constructive feedback.

• Reviewer’s comment: “This work focuses on learning contractive dynamics. However, suppose the system $\dot{x} = f(x)$ is not contractive. What can you say about the solution to the optimization problem (i.e. the feasibility of eqns (4) and (5)) in this scenario?”
  • Short answer: Good question. When the system dynamics $\dot{x}=f(x)$ are not contractive, the dynamics described by equations may integrate to arbitrary (potentially unstable) dynamics. The system will most likely follow the trends of the training data when evaluated near the data, and then exhibit arbitrary behavior everywhere else (thereby violating stability guarantees).
  • Detailed answer: If the system dynamics $\dot{x}=f(x)$ are not contractive, i.e. the Jacobian of the function $f$ is not uniformly negative definite, the trajectory described by equations (4) and (5) may not exhibit stable behavior. It is worth noting that these equations do not depend on the Jacobian being uniformly negative definite. The optimization of equations 4 and 5 without considering the contraction formulation in equation 3 means the system might be only able to reconstruct the trajectories if started on the data support regions (not guaranteed) but the global behavior of the system outside of the data support will not be stable. In Fig. 13 in the updated paper, you can see the differences between contractive and unconstrained dynamical systems. As the figures show, when the resulting vector field is not contractive, the reproduced motion follows the trend of the data only on the demonstrations region and not on the outside regions. Overall, this system does not show any stable behavior.
  • Revision: We have added a new subsection in Appendix A.6.3: Unconstrained Dynamical System discussing these results.

• Reviewer’s comment: “In a similar vein, it seems like this work could easily be extended to (neural) controller synthesis in the latent space. Can the authors comment on this?”
  • Our interpretation of this comment is the potential extension of the current work to incorporate a controlled formulation of the dynamical system, i.e. $f(x,u)$, where $u$ is the control input**.** If this is correct, then it is plausible to envision formulating the system within the latent space. The latent space provides a compact yet semantically rich representation of the system dynamics that captures essential features of the underlying system in a more abstract and expressive manner. Extending the model to integrate a neural controller within this latent space framework is a neat idea for control strategies at a higher level of abstraction. The method could then explore the synergy between the latent space representation and the dynamics of the controlled system. This is a very interesting idea, but not one we have explored

• Reviewer’s comment: “Can the authors comment on the practical effectiveness of the model in greater detail? As mentioned in the discussion section (p9), the cost of numerical integration can be extensive. Have the authors come across examples where this has been an impediment to performance?“
  • Short answer: No such examples have been encountered.
  • Detailed answer: numerical integration is one of four steps affecting performance, where the others are: (1) encoding the current state, (2) computing the Jacobian, and (3) calculating input space velocity. Each iteration in our robot control scheme (Algorithm 2) takes up to 10 milliseconds (100 Hz). Current experiments involve 7 DOF robots, encoding either pose (6-dimensional) or joint configuration (7-dimensional), with a 2-dimensional latent space. With respect to the cost of the numerical integration, then this depends mostly on the latent dimension and less so on the input dimension (Table 3 in the appendix). Finally, it is worth noting that our Python implementation is optimized for correctness rather than performance, so there is substantial room for performance improvements.

• “How is the collision-avoidance on the entire manipulator handled? The learned dynamical system seems to model the end-effector pose, but collision-avoidance should be handled across the body of the robot. One approach to handle this is to pull the dynamical system to the C-space of the robot and define body-points for the collision-avoidance, as done in (Zhi 2022, L4DC). This is a relevant reference and should be reviewed, as it also takes a diffeomorphic learning approach.”
  • Short answer: Thank you for bringing this to our attention; we indeed overlooked this aspect in the paper.
  • Detailed answer: The paper mentioned by the reviewer [2] is indeed interesting: it ensures collision avoidance in the C-space by leveraging Diffemorphic Transforms (DTs). Specifically, it employs smooth obstacle gradients obtained from continuous maps based on occupancy information to construct a vector field. Then, it uses the natural gradient of this vector field to smoothly warp around the obstacle when approaching its boundary. This is similar to obstacle avoidance approach proposed in Beik-Mohammadi et al. [3] where the C-space obstacle avoidance is accomplished by employing a pull back metric derived from the VAE decoder. This approach formulates the C-Space obstacle avoidance in the latent space of the VAE instead of the task space. To clarify, as mentioned in the paper our current modulation matrix is formulated in the input space of the VAE, therefore, the ability to perform multiple limb obstacle avoidance can be achieved when using joint space demonstration trajectories. However, to handle obstacle avoidance across the robot body we need to transform the joint space information into the task space and define body-points for multiple-limb collision-avoidance. Of course, this will be computationally expensive. Certainly, adopting the obstacle avoidance concepts introduced in [2] or [3] eliminate the necessity for a meticulous computation of the modulation matrix, particularly for multiple limb obstacle avoidance. Moreover, these methods allow for the transition of the obstacle avoidance task into the latent space, mitigating any impact on real-time performance due to dimensionality of the input space. While these approaches show potential in refining our obstacle avoidance capabilities, it is important to note that obstacle avoidance is not the central focus of our work. Nevertheless, they play a valuable role in the broader context of exploring potential enhancements in future works.
  • Revisions: We have clarified C-space obstacle avoidance in the Section 3.4.

• References:
  • [1] Johann Brehmer and Kyle Cranmer. Flows for simultaneous manifold learning and density estimation. In Neural Information Processing Systems (NeurIPS, pp. 442–453, 2020. URL https://proceedings.neurips.cc/paper_files/paper/2020/file/051928341be67dcba03f0e04104d9047-Paper.pdf.
  • [2] Weiming Zhi, Tin Lai, Lionel Ott, Fabio Ramos . Diffeomorphic Transforms for Generalised Imitation Learning. Proceedings of The 4th Annual Learning for Dynamics and Control Conference, in Proceedings of Machine Learning Research. 168:508-519 Available from https://proceedings.mlr.press/v168/zhi22a.html.
  • [3] Beik-Mohammadi Hadi, Hauberg Soren, Arvanitidis Georgios, Neumann Gerhard, Rozo Leonel. Reactive motion generation on learned Riemannian manifolds. The International Journal of Robotics Research. 2023;42(10):729-754. doi:10.1177/02783649231193046. https://journals.sagepub.com/doi/10.1177/02783649231193046

We appreciate the reviewer's engagement with our work and the valuable feedback provided.

• "Theorem 1 states that the contractivity is perserved under a diffeomorphism. There is a lack of analysis and formal guarentees around the assumption that an injective function acting on a lower dimensional contractive system produces a higher-D contractive system.“
  • Short answer: We agree that the paper was not precise enough here. Shortly put, the injective decoder is a diffeomorphism between the latent space and the image of the decoder, which is sufficient to ensure that the resulting system is contractive. We detail this argument below.
  • Detailed answer: Let $\mu: \mathbb{R}^d \rightarrow \mathbb{R}^D$ be the injective decoder mapping from the low-dimensional latent space to the high-dimensional observation space. Further, let $z_t$ denote a contractive dynamical system in the latent space. Geometrically, the decoder spans a $d$-dimensional manifold embedded in $\mathbb{R}^D$ defined as $\mathcal{M} = f(\mathbb{R}^d)$. By design, the decoder is a diffeomorphism between $\mathbb{R}^d$ and $\mathcal{M}$; this construction is extensively discussed by Brehmer & Cranmer [1], which we build upon. By Theorem 1 in our paper, the above implies that $f(z_t)$ is contractive. The above further implies that the dynamical system $f(z_t)$ is along a low-dimensional manifold $\mathcal{M} \subset \mathbb{R}^D$, i.e. it does not cover all of $\mathbb{R}^D$. Such constructions are common in latent variable models. The control loop assumes that the initial robot configuration is in $\mathcal{M}$ such that it can be written $f(z_0)$. Following the dynamical system moves the robot along the manifold, and the resulting motion follows a contractive system. If the assumption is not satisfied, such that the initial configuration $x_0 \not\in \mathcal{M}$, then we have to rely on the encoder to approximate a projection onto $\mathcal{M}$, i.e. produce $z_0 = \mu^{\sim 1}(x_0)$. Note that the movement from $x_0$ to $f(z_0)$ need not be contractive, such that the distance to the nominal trajectory may increase. This is, however, a finite-time motion, after which the system is contractive.
  • Revision: In Section 3.2, we have clarified and formalized the correlation between injection and contraction.

• “Just as a thought experiment, what happens to coordinate points in the higher-D data space where there does not exist a coordinate in the lower-D latent space?”
  • This matches the last remarks of the detailed answer above. The model provides a contractive dynamical system on top of a $d$-dimensional manifold embedded in $\mathbb{R}^D$. When the initial robot configuration is not on this manifold, we first have to project it there. We approximate this with our encoder $\mu^{\sim 1}$ and note that the finite-time motion from the initial configuration to the one on the manifold need not follow a contraction. The robot video in the supplementary material includes an example of this. Here, the operator often pushes the robot away from the manifold spanned by the decoder after which the robot first moves back on the manifold and thereafter follows a stable behavior.

Thank you for your interest in our work. We answer your question below.

“I had one question on the experimental results, namely the learned vector fields in Figure 5. According to the second sentence of Section 4.1, no latent structure was used. Thus, the learned vector field is contractive in the 2-norm in the sense of Definition 1. Accordingly, any two trajectories of the vector field satisfy, i.e., any two trajectories converge to one another exponentially fast without any overshoot (see Theorem 3.9 in [1]). Pictorially, there appear to be trajectories that start close to one another and then grow further apart before converging together near the equilibrium point. This seems to violate the condition that there is no overshoot (in the 2-norm). Are the authors applying extra invertible transformations on top of the learned dynamics to yield vector fields which are contractive in more general Riemannian metrics?”

Short Answer: No, the observed behavior is not due to the extra invertible transformations on top of the learned dynamics.

Detailed Answer: The observed behavior in Figure 5 is a direct result of our hyperparameter tuning, specifically the regularization parameter, a critical factor influencing the system's contraction rate. As a side note, Equation 2 in the paper characterizes the contraction. In this example, we set the regularization term (in Equation 3 of the paper) to $10^{-4}$, thereby restricting the maximum negative eigenvalue of the Jacobian of the contractive system to this small value, therefore, estimating a negative definite Jacobian. Consequently, the system exhibits a less pronounced contractive behavior, strategically designed to enhance its ability to generalize across multiple demonstrations. This deliberate choice prevents the system from abruptly converging to a single narrow region, achieving adaptability and robust performance. The rationale behind this line of thinking is extensively discussed in Appendix A.6.3 on Regularization. Moreover, we successfully experimented with increasing this parameter, resulting in a demonstrably strong contractive dynamical system, as illustrated in our new experiment detailed in Appendix A.6.2 (Figure 12). However, given the inherent nature of the task, which involves generalizing to multiple trajectories, a decision was made to opt for a smaller regularization. It is worth noting, that the regularization in Equation 3 in the paper affects all the diagonal values of the Jacobian matrix (also the values that are already negative), therefore, a more targeted approach may improve the results.

Additionally, a comparable behavior can be seen in other contractive dynamical system methodologies when subjected to similar experiments. A noteworthy instance can be found in Figure 3 in [1], illustrating the quantitative performance of the CDSP algorithm on the LASA dataset. In this context, trajectories occasionally exhibit minor divergence before converging at the so-called target, again that is not a target and only a part of the contractive region. Furthermore, this behavior can be seen in the contours of the contractive vector field in CVF approach [2] in Figure 2. Similarly, you can observe the same behavior in the C-GMR behavior in Figure 4 [3]. We believe that such experimental insight in the cited papers could be possibly explained by very low contraction rates in their learned systems.

[1] Ravichandar, H.c., Dani, A. Learning position and orientation dynamics from demonstrations via contraction analysis. Auton Robot 43, 897–912 (2019). https://link.springer.com/article/10.1007/s10514-018-9758-x

[2] Sindhwani, Vikas & Tu, Stephen & Khansari, Mohi. (2018). Learning Contracting Vector Fields For Stable Imitation Learning. https://arxiv.org/abs/1804.04878

[3] C. Blocher, M. Saveriano and D. Lee, Learning stable dynamical systems using contraction theory, 2017 14th International Conference on Ubiquitous Robots and Ambient Intelligence (URAI), Jeju, Korea (South), 2017, pp. 124-129, https://mediatum.ub.tum.de/doc/1364708/1364708.pdf