Input-to-State Stable Coupled Oscillator Networks for Closed-form Model-based Control in Latent Space
Authors: Maximilian Stölzle
Cosimo Della Santina
Venue: NeurIPS 2024 spotlight (2024)
Keywords: Dynamical Systems
Control Theory
Robotics
Decision and Control
Deep Autoencoders
Primary area: robotics
TL;DR: We leverage input-to-state stable coupled oscillator networks for conducting model-based control in latent space.

Abstract: Even though a variety of methods have been proposed in the literature, efficient and effective latent-space control (i.e., control in a learned low-dimensional space) of physical systems remains an open challenge.
We argue that a promising avenue is to leverage powerful and well-understood closed-form strategies from control theory literature in combination with learned dynamics, such as potential-energy shaping.
We identify three fundamental shortcomings in existing latent-space models that have so far prevented this powerful combination: (i) they lack the mathematical structure of a physical system, (ii) they do not inherently conserve the stability properties of the real systems, (iii) these methods do not have an invertible mapping between input and latent-space forcing.
This work proposes a novel Coupled Oscillator Network (CON) model that simultaneously tackles all these issues. 
More specifically, (i) we show analytically that CON is a Lagrangian system - i.e., it possesses well-defined potential and kinetic energy terms. Then, (ii) we provide formal proof of global Input-to-State stability using Lyapunov arguments.
Moving to the experimental side, we demonstrate that CON reaches SoA performance when learning complex nonlinear dynamics of mechanical systems directly from images.
An additional methodological innovation contributing to achieving this third goal is an approximated closed-form solution for efficient integration of network dynamics, which eases efficient training.
We tackle (iii) by approximating the forcing-to-input mapping with a decoder that is trained to reconstruct the input based on the encoded latent space force.
Finally, we leverage these three properties and show that they enable latent-space control. We use an integral-saturated PID with potential force compensation and demonstrate high-quality performance on a soft robot using raw pixels as the only feedback information.

Input-to-State Stable Coupled Oscillator Networks for Closed-form Model-based Control in Latent Space
Maximilian Stölzle
Cosimo Della Santina Department of Cognitive Robotics Delft University of Technology
Department of Cognitive Robotics
Delft University of Technology M. W.
C. DellaSaat ina@t udelft. n1
Abstract
Even though variety of methods have been proposed in the literature, efficient and effective latent-space control (i.e., control in a Iearned low-dimensional space) of physical systems remains an open challenge. We argue that a promising avenue to leverage powerful and well-understood closed-form strategies from control theory literature in combination with learned dynamics, such as potential-energy shaping. We identify three fundamental shortcomings in existing satent-space models that have so far prevented this powerful combination: (i) they lack the mathematical structure of a physical system, (ii) they do not inherently conserve the stability properties of the real systems, (iii) these methods do not have an invertible mapping between input and latent-space forcing. This work proposes novel Coupled Oscillator Network (CON) model that simultaneously tackles all these issues. More specifically, (i) we show analytically that CON is a Lagrangian system i.e., it possesses well-defined potential and kinetic energy terms. Then, (ii) we provide formal proof of global Input-to-State stability using Lyapunov arguments. Moving to the experimental side, we demonstrate that CON reaches SoA performance when learning complex nonlinear dynamics of mechanical systems directly from images. An additional methodological innovation contributing to achieving this third goal is approximated closed-form solution for efficient integration of network dynamics, which eases efficient training. We tackle (iii) by approximating the foreing-to-inpel mapping with a decoder that is trained to reconstruct the input based on the encoded latent space force. Finally, we leverage these three properties and show that they enable latent-space control. We an integral-saturated PID with potential force compensation and demonstrate high-quality performance on a soft robot using raw pixels as the only feedback information.
Introduction
Learning how environment evolves around us from high-dimensional observations (i.e., world models (1) is essential for achieving both artificial and physical intelligence [2]. For example, world models are required for effectively planning an artificial/robotic agent's 's actions in complex and unstructured environments (3]. However, learning such dynamics directly high-dimensional observation space usually intractable. Seminal works have shown that we can leverage autoencoders to compress the state information into a low-dimensional latent space [4]5] in which it is much more feasible to learn the dynamics [6] to using these learned models generate low-level intelligence.
10]. However, strong limitations still persist when it comes
One outstanding challenge is how to perform closed-loop control in the learned latent space ie., how to generate control inputs based on high dimensional sensory input such that desired movement is generated. Prior works have explored, among other approaches, Reinforcement Learning (RL)
11 13 14], Model Predictive Control (MPC) [7 15 16 17], Linear-quadran Regulators (LQRRs) [8 19 20 and gradient-based optimization [2] for planning and control towards target evolution that given observation space. However, all existing latent-space control strategies have
38th Conference on Neural Information Processing Systems (NeurIPS 2024).

---

Detutee
Purrel porlacied
SA D. diace
Euroeler
Demucler Besusitracee
Fylt
Emper wat
repurt uiz)
NAA
Coepled Gerilaaor Netaork
Predicte abwersatiie
M
Takegg
cet ther way
lukia di alfe)
Fu00
62%)
(a) Coupled Oscillator Network (CON)
(b) Learning patent-space dynamics with CON
Figure 1: Panel (a): The proposed CON network consists n damped harmonic oscillators that are coupled through the neuron-like connection tank W x++) and the non- non-diagonal stiffness K - and damping coefficients D-d, respectively. The state the network is captured by the positions xx() and velocities =()) of the oscillators. The time-dependent input mapped through the (possibly nonlincar) function g(u) a forcing T acting on the oscillators. Panel (b): Exploiting Coupled Oscillator Networks (CONs) for learning latent dynamies from pixels: We encode the initial observation o(fa) and the input u(t) into latent space where leverage the Coupled Oscillator Network (CON) predict future latent states. Finally, decode both the latent-space torques (()) and the predicted latent states 2t)..
shortcomings, such as limited planning horizon and slow control rates (MPC and gradient-based approaches), sample inefficiency (RL), or they pose a requirement for learning linear dynamics (6]
(LQR), which is not even possible for systems that are inherently non-linearizable [22]. One
interesting avenue is to leverage model-based control approaches, such as potential shaping [3] 24]
10), for effective and computationally efficient control in latent space (25] For these techniques to be feasible, the dynamical model needs to fulfill four characteristics: (i) the dynamies need to have the mathematical structure of physical systems, (ii) conserve the stability properties of real systems, (iii) the latent state needs to be relatively low low-dimensional, and (iv) there needs to exist a well-defined, invertible mapping between the input and the forcing in latent space. However, existing model structures that are used for learning latent dynamics (26) do not meet all of these criteria. Relevant examples are Multilayer Perceptrons (MLPs), Neural ODEs (NODEs) (27 28), many variants of Recurrent Neural Networks (RNNs) (e.g.. LSTMs 29). Gated Recurrent Units (GRUs) (0), etc.), and physics-informed informed neural networks (e.g.. Lagrangian Neural Networks (LNNs) [1 32,
Hamiltonian Neural Networks (HNNs)) B3) For example, MLPs do not have a physical interpretation and do not provide an invertible mapping of the forcing generated by the input, NODEs are usually
not easily stabilizable (34), most RNNs require a relatively nigh-dimensional latent space (i.e., many
hidden states), and energy- energy-shaping control approaches based on LNNs 10 do not come with any formal stability guarantees.
have been shown exhibit of-te-of-the-art performance on time sequence modeling tasks while being parameter-efficient thus fulfilling our requirement (iii). Consequently, we believe that they are a promising option for control-oriented dynamics learning in latent space. Still, these models do not fulfill the remaining requirements that we have listed above. Despite being an interpretable combination of harmonic oscillators, they do not have the structure of a physical system i.c., they do not possess a well-define energy function. Moreover, only local stability [35]3]] has been shown, with sufficient conditions that appear to be very stringent. Finally, in addition to training an encoder that maps inputs to latent-space forcing, we propose also training a decoder that learns to reconstruct inputs based on aatent-space forcing. This enables us to easily switch between inputs and forcing, which is essential when implementing control strategies.
In recent years, oscillatory networks B5 36 37 38
39
We resolve all the above-mentionel challenges by proposing Coupled Oscillator Networks (CONs), a new formulation of a coupled oscillator network that is inherently Input-to-State Stability (ISS) stable, for learning the dynamics of physical systems and subsequently exploiting its structure for model-based control in latent space. The network consists of damped, harmonic oscillators connected through elastic springs, damping elements, and a neuron-like like coupling force and can be excited by a nonlinear actuation term. We identify transformation into set of coordinates from which we can derive the networks' kinetic and potential energy. This allows to leverage Lyapunov arguments [40] for proving the global asymptotic stability of the unforced system and ISS stability for the forced system under relatively mild assumptions on the network parameters. Even though we constrain the dynamics to a very specific structure, we demonstrate (a) the CON network achieves similar performance as NODES when learning the dynamics of unactuated, mechanical systems with two

---

orders of magnitude fewer parameters and (b) that the proposed model achieves, for the complex task of learning the actuated, highly nonlinear dynamics of continuum soft robots [4] 42 from pixels. 60 % lower prediction enror than Coupled Oscillatory Recurrent Neural Network (coRNN) [35] and reaches the SoA performance across all techniques that we tested. Finally, we show some initial results that the proposed CON model is also able to learn the latent dynamics of Partial Differential Equations (PDEs), in this case containing caactionddiffsion
directly
9 dynamics.
Subsequently, we derive approximate closed- solution, that is, in parameter regimes in which the linear, decoupled dynamics dominate transient, more accurate than numerical integrators with comparable computational requirements and which increases training speed by 2x with a small decrease in prediction accuracy. Finally, as we can derive the system's 's potential energy, we can leverage potential shaping [23] 24 to derive controller that combines
integral-saturatd PID controller with feedforward term compensating potential forces. As the feedback acts on well- shaped potential field, tuning the feedback gains becomes very simple and out-of-the-boo and the controller exhibits a faster response time and 26 % lower trajectory tracking Root Mean Squared Error (RMSE) than a pure feedback controller based
latent NODE (27 model.
The proposed methodology is particularly well-suited for learning the latent dynamics of mechanical systems with continuous dynamics, dissipation, and single, attractive equilibrium point. Examples of such systems include many soft robots, deformable objects with dominant elastic behavior, Lagrangian systems immersed in a dominant potential field, or locally other mechanical systems such as robotic manipulators, legged robots, etc. For these systems, we can fully leverage the structural prior of the proposed latent dynamies, including the integrated stability guarantees. If the system actuated, the learned dynamics can be subsequently exploited for model-based control, as demonstrated in Sec.
The code associated with this paper is available on GitHuut
2 Input-to-State Stable (ISS) Coupled Oscillator Networks (CONs)
Formulation.
36 37 38 are one- dimensional, potentially damped, harmonic oscillators, which are described by their state y/ = [2, E."" E IR², where z; and é; are the position and velocity of the oscillator, respectively. Then, the oscillator's dynamics are defined by the following Equation of Motion (EOM)
The integral component to (coupled) oscillatory RNNs 35
my #1()) de dt)(() + Ky sr.(1) F,((), with mi, Ka; d, E IR+
(I)
Here, my is the mass, is the stiffness, and d, is the damping coefficient of the damped harmonic oscillator. (t() E IR is (possibly time-dependent) external forcing term acting on the mass.
Even though the state is extremely low dimensional and the number of parameters is small, this single, damped harmonic oscillator can already exhibit a variety of (designable) behaviors: The expressions
d; let us determine the natural frequency and the damping factor, ra and C 2 mn
Wa.s
V
respectively and allow us to design the transient behavior. For example, Wa,i lets us isolate spectrum of the input signal (.(1) (37] and ç, determines the damping regime: underdamped (wi., < 1), critically damped (wn.; 1), overdamped (wn.; >> 1). Furthermore, as (damped) harmonic oscillators are omnipresent in nature (and especially in physical systems), they have been intensively studied and well understood (e.g.. characteristics, closed-form solutions, etc.). In this work, we will exploit some of these properties and knowledge to learn stable (latent) dynamics efficiently.
By intercoupling damped harmonic oscillators, we can drastically increase the expressiveness of the dynamical system (35 [37] 38 while preserving some of the intuition and understanding we have for these systems. In this work, we propose a stalstable CON consisting of damped harmonic oscillators that are coupled through both linear and nonlinear terms. The networks' state is defined T
as y = T &
E 22tn and its dynamics can be formulated as a 2nd -rder Ordinary Differential
Equation (ODE)
i(t)=
dz
((v()),(0)) = (
*(t) Kx(() Da(() tanh x(t) + +b]
(2)
where K, E R are the linear stiffness and damping matrices, respectively. The neuron-inspired term Lanh (W =1{)= + b) with W E R=*xn_ & E R" provides nonlinear coupling between the harmonic oscillators. The network is excited by the time-dependent input u(t) E I through the possibly nonlinear mapping g : 1T + I2". Specitically, we consider in this work a formulation where an input dependent matrix
attps:///ittae cowpttddcccomuuoeeeeeeoo phipppcc///ee www.///cccevinggeeeeeeeee iss coupled oscilletor aetworks fron- pixele
3

---

Buu) E R=xm projects the imput u(t) to time-dependent forcing the oscillators: g(u) (u)) Here (u)) could, for example, be parametrized by MLP.
We specifically designed the network architecture such that (i) the system exhibits a unique and isolated equilibrium and (ii) we can derive expressions for the kinetic and potential energies. These two features allow us (a) prove Global Asymptotic Stability (GAS) and ISS stability using an established procedure based on striet Lyapunov arguments [44 45 and (b) implement model-based controller based on potential shaping.
One key insight of this work is that in the coordinates (()), cannot derive potential the hyperbolic force (aah(Wa(tt)) +b) is not symmetric. Therefore, propose coordinate transformation into W-coordinates:
2 (t)
W x(t) W a(t)
yo (1)
E RRn. The coordinate transformation valid if its Jacobian full-rank, which
2 (t)
the case if rank rank(W) = na. In W-coordinates, the dynamics can be rewritten as
dite
S. (t)
U (f)
dr
fm(u(t), Mw (g(u(()) - Ku (1) Du i (t) (r. (t) + b]]
(3)
with Kw K W. Da. DW and My W
A difference of this formulation compared to prior work 35 36 37
38
is that (i) the forcing produced by the
input term g(u) is fully separated from the forcing produced by the clastic coupling terms Kw. and (ii)
the generalized force is symmetric, which we prove in Appendix [A.1] allowing us to define a potential energy expression, which we can later on leverage for stability analysis and control.
equilibria Fw = F oTT²² E IR the unforced network are given by the roots of the characteristic equation tanh(Tw + b) + Kw Jow = O.
Lemma Ler Kn [
Then, dynamics defined have a single, isolated equilibrium yo
0"]T.
Proef. The proof is straightforward and provided in Appendix A.2
Next, we introduce a mapping into the tilde coordinates fw = cquilibrium (w) can now be stated as
Jw. The residual dynamics (w.rt. the
fhe (t) (v) M (o(u(t)) Kw (Iw + Z. (1)) D. (t) tanh(# (1) bn)]
i (t)
(4)
the following, we will write ||A| || to denote the induced norm of matrix and (m))),, (M((A) to refer to its minimum and maximum Eigenvalue respectively.
Global Asymptotic Stability (GAS) for the unforced system. We first consider the unforced system with 7 = g(u) = 0, E (tu, t and strive to prove global asymptotic stability (40) for the attractor ž. We propose a strict Lyapunov candidate with skewed level sets (45)
(()) =
2 Ve Pv f +
Σ
tanh(T + a by) do
E) tanh(Tw. + b.) do,
2 v fm
Σ
(leosh (Tw.s Sw.. by) lcosh(ā b.) tanh(āw. + bi) Sw.)), (5)
with Py
I MT
Kw Mw]
M
E 2m222n
lcosh(-) = log(cosh(.)). and μ 0.
Lemma 2. The scalar finction (iw) defined in (5) is continuously differentiable and verifies the condition
(0) = 0. Furthermore, let Mw, Kw > O. Now, if we choose 0 < y Ama(dd (m()K) = pvv. then
[|
(ū)) o4g 222 \10). Additionally, then (ū(ū)) is radially unbounded as || = Vū
Proof. We provide the proof in Appendix Lyapunov candidate he positive-definite
A.3 and demonstrate that bounds on μ are required for
Theorem 1. Let Mw, Kw and Dw be positive definite and suppose the system be unforced: ( = -.. Then,
Tw = is globally asymptotically stable for the system dynamics defined (4) such that V ((jw) < 0. "w R (0).

---

Prouf. First, we show that yh- 0 is equilibrium of (
Fo(0.0) M (-Kw Jw tanh( + b)]
0 M (-Kw Iw + Kw Tww
0 [
(6)
Lemma states that we can always choose js such that 5 is strict Lyapunov function (c.g., Lipschitz continuous,
zero-valued at Fh = O, positive definite, and radially unbounded) [40 45]. We now evaluate the time-derivative
of Va (fw) in the case of an unforeed system (i.c., g(u(t)) 0):
Ny
@y dyw
fw((ww)
The Pv DN (tanh((æw + Ew + b) tannl( + b))* iw
you
1 Kw
Vw
2 D
Dw r Mw
Uw [aalh(@w + Ew b)" wee Kw) Tww.
Vw Pv Nw μ (Lanh((ww + Iw + Lanh((@ + 6))" Tw:
S - 5w Pv Nw s -m (Pv) 11
where we exploited the force balance at equilibrium Kw Fy = tanh(āw + b) and Lemma for defining the upper bound Vz (Uw). Lemma states that Pv 0 for μ < Mé-- Similarly, Lemma requests that < a < jay. Indeed, both conditions can always be fulfilled by choosing E (0, min(fv, Pvy /v)) With Pv > 0 + Am (Py) > [46]. we can state that V (ūw) 4 Dw E BR2 \ {0} and conclude that the unforced system is globally asymptotically stable around Dw
Global Input-to-State Stability (ISS) for the forced system. We now take the forcing g(u) into account again and demonstrate that the system states remain proportionally bounded to the initial conditions and function of the supremum of the input forcing.
Theorem Suppose Ma, Kw; Dw > O, 8 and that we choose < min(wv, Mb Then, is globally Input-to-State Stable (ISS) such that the solution Th (1) verifiex
1 lla S 3 (to)ll,, to)
sup Io((((())ll)) case'tt
AA to
(8)
- (athe) (M()))) 11 Am (P)) K. 82 Am (Pm) (Pp)
where B(r,t) KCC. v(r)
We first demonstrate the ISS-Lyapunov function is bounded from both sides by K. functions. Subsequently, derive the energy dissipation of the forced system and establish attracting regions function of the norm of the forcing. Outside these regions, it ensured that the system has minimal rate decay and, therefore, converges exponentially fast into the attracting region.
Proof. The proof is provided Appendix A.S
0.75
0.50
Velocity [a//]]
0.25
Positica
0.00
0.25
ground-truth
-0.2
CON Thi5 US-
-0.50
CON Euler 0.05x
CON Enler di 0.06s
CFA
CON CEA-CON
de 0.1s
-0.4
CFA CON & 0.1s
0
10
20
30
80
10
20
30
40
Time 8)
Tine (+)
(a) Positions
(b) Velocities
Figure 2: Analysis of approximation error of CFA-CON: we compare the ground-truth solution a s rollout of the CON network consisting of three oscillators (n 3) with the CFA-CON executed at time step of = 0.1 and solution generated by integrating the ODE at time step of 0.05 s with the Fuler method.
5

---

3 An approximate closed-form solution for the rollout of CON
To predict future system states, we need to integrate the ODE in Eq. 22. with the solution given by y(fx++1) = yt,, + S1h+ [v(('), u((')) di'. Unfortunately, a closed-form solution for the nonlinear dynamics fly, u) does not (yet) exist. Therefore. we traditionally need to revert nrigh-order) numerical ODE solvers that computationally very expensive and introduce additional memory overhead [47]. This considerably increases the training time of models involving such continuous-time dynamies. While the computational time can be reduced by increasing the (minimum) time step of the integrator. this comes at the expense of an integration error, and we lose (part of) the theoretical guarantees and practical characteristics that the nominal ODE provides. In this work. we take an alternative approach by splitting the problem into (i) decoupled linear dynamics that can be cheaply and precisely integrated using a closed- form solution and (ii) the residual, coupled nonlinear dynamics. w which we integrate numerically at a slower time scale:
((t)
F (x)(t) d'i(t) + g(((t) - (K x)x(()))))))))))))<)))))) - (D d) a(t) tanh (W/ xx1) + bb)
(..a(lV):dccoppel linear dynamics
aad (v. coupled, nonlinesr dynumics
where x diag K1 ... n. d (Dui... Dan) are the diagonal components of the stiffness and damping matrices, respectively, and F E R is constant, external forcing term the oscillators.
For short-time-intervaal öt, we now approximate (
with F(ts) = - -f (v(ta), aitt)). For a scalar -oddrr lincar ODE of form ús (x(t),,, ) F(tx)),
#(tx + 5t) = fa.a(v((ta + 6t), F(x)))
(10)
well-known, closed-form solution exists. We exploit this characteristic by formulating the approximate solution as
48
utta + 5t) 2 fCFA (CrA-cCoNy(tt), u(ta)) = v(ta)= + S.a +5t
fial(tt), dt'
(11)
and denote SCFA-CON R" X I + 1e" as the Closed-Forts Approximation of the Coupled Oscillator Network (CFA-CON) model. The implicit assumption bchind (TO) is that fa.a(l) > (y, u) (i.c. the lincar, decoupled dynamics dominate the nonlincar, coupled, time-varying dynamies). We refer the interested reader to Appendix for derivation and implementation details, where we summarize the integration procedure in Al Algorithm]] I We also provide qualitative results for the integration accuracy in Fig. and quantitative results for the integration accuracy and computational speed-up w.r.t. numerical integrators in Appendix B
4 Learning control-oriented latent dynamics from pixels
We now move towards learning latent dynamical models based on CON and CFA-CON. CONs are an ideal fit for learning latent dynamics as they guarantee that the latent states stay bounded.
We assume to have access to observations in the form of images R R Xm% where c. denotes the number of channels. Please note that this could also be other high-dimensional observations such as LiDAR scans, point clouds, etc. We now leverage an encoder-decoder architecture to map these high-dimensional observations into compressed latent space: The encoder dimensional latent representation E R of the images. The decoder V/ R + B X approximates the inverse operation by reconstructing an image ô R x based on the latent representation. To promote the learning of smooth and monotonic mapping into latent space, we specifically choose to implement the autoencoder here as 3 Variational Autoencoder (VAE) [ 49]. Instead just statically reconstructing the image ö(t4). ). we are interested in predicting future observations 0(lx41). where E 1... N. For this, we train a 2nd 2"doorder dynamical model that is, when integrated, able to predict future latent representations z(fa+1). This
R X xe. + R with Ti, < ho Wh identifies a low-
requires us to define latent state E(t) [x (t) 4 [t)] E 2=² consisting of the latent representation and latent velocity (t) E R"=
È E fy (t)). where we defined { Yo. and = = Tw. To ensure stability, we make use of the Cholvesky decomposition to ensure that Ma. K2r and D., always remain positive definite (see Theorem 2) It is important to note that we train the encoder, decoder, and dynamical model all jointly. Please refer to Appendix for more implementation details.
We now rely on CON with n. oscillators to provide us with the latent state derivative
Training. It important to remember that because are using a B- VAE 49 the image encoding becomes stochastic, and the encoder neural network actually oulputs p= (o), 2 log(o=) E R=* After executing the reparametrization trick as 2(tx) ~ ((xx(ta) (tx)), we formulate the loss function, evaluated on cach trajectory consisting of N time-steps, as
MSB(o(())) (=(te))
PKLL ((ns (e)) (l)) ((+)) N Kallbakk Leitter divergence
MSB(a(4)) ((EL)
MSB(= R((R
N Danamtic inigg Poosadaaciion Ioss
State COOOEOEAAAA
LEaEE dytraties conpiitces los

---

Model
RMSE M
20
RMSES
20
RMSE
D-P+ 26
RMSE CS
RMSE PCC-N-82
RMSE
PCCNS-3 4
RNN
0.2739
= 0.0057
0.2378 + 0.0352
0.1694
0.0004
0.1011
+ 0.0009
0.1373 + 0.0185
0.2232
0.0075
GRU [30]
0.0267
0 0033
0.1457 1: 0.0078
0.1329
00005
0.1125-
1: 0.0110
0.0951 1. 0.0021
0.2148
0.0196
ooRNN 130
0.0265-
0.0002
0.1333 1: 0.0044
0.1324
0.0016
0.2537
: 0.0018
0.2501
1: 0 0899
0.2474 0.0018
NODE 22]
0.0264
0.0010
0.1260
+ 0.0013
0.1324
0.0024
0.2415
0.0021
0.1867
+: 0.00561
0.3373
0.UD65
MRCH-NODE
0.0328
0.0034
0.1650
+ 0.0205
0.1710
- 0.0111
0.2494
4 0.0028
0.1035
: 0.0012
0. 0.1900 = 0.0024
CON 5 tour)
0.0303
0.0053
0.1303
1: 0.0064
0.1323
1: 0.0018
0.1993
4 0.0646
0.0996
1: 0.00012
0.1792
0. 0038
CON-M (our)
0.0303-
0.0p53
0.1303
0.00G4
0.1323
0.0018
0.1063
+ 0.0027
D.1008
+ 0.000G
0.1785 0.0023
CEN-OON lour)
0.0313 0.0026
0.1352
+ 0.0075
0.1307
0.0012
0.1482
0.0211
0.1124
+ 0.0020 0.1803 0.0003
Table 1: Benchmarking of CON and CFA-CON learning latent dynamics against baseline methods. The first three datasets, based on (26), contain samples of a mass-spring with friction (M-SP+ +F,++ single pendulum with friction (S-P+ F), and a double pendulum with friction (D-P + 1") (all without system inputs). The CS dataset considers continuum soft robot consisting of one segment with three constant planar strains. The PCC-NS-2 and PCC-NS-3 datasets contain trajectories of continuum soft robot made of two and three piecewise constant curvature segments. respectively. We choose the latent dimensions of the models as
4. n2 4, and R2 = 12 for the M-SP + F, S-P + F, and D-P + F datascts, and n2 = 8, n2 12, and n; = 12 for the PCC-NS-2, 2, PCC-NS-3, 3, and CS soft robotic datascts. We report the mcan and standard deviation over three different random seeds.
Dataset
RNN
GRU
30
coRNN
NODE
MECH-NODE
CON-S (our)
CON-M (our)
CFA-CON (our)
M.SPoF
4
88
248
40
3368
3244
34
34
34
D-P+F
12
672
1968
348
4404
4032
246
248
246
PCC-NS-2
320
928
152
3856
3062
676
7048
7048
Table 2: Number of trainable parameters for the various latent dynamic models and the examples of the M-SP+F (n. unactuated). D-P+F (n. = 12, unactuated) and PCC-NS-2 (n. = 8, actuated) datasets. The number of trainable parameters for all models and datasets can be found in Appendix D
where s(tx) is predicted by E(tx)= Sh Je(E(t') u((')) dt', and E(to) = (to) ='(tn)] Here, (to) is given by the encoder, and (fo) is approximated using finite differences image-space (see Appendix C.5 for more details). ,,, Ao, X₂ E are loss weights.
Models. We train the CON with the input-to-forcing mapping g()) = Buu) u. where B(u) parametrized by a MLP with a hyperbolic tangent activation function applied between layers. We report results for two variants of the CON model: for the medium-sized CON-M and small-sized CON-S, the MLP consists of five and two layers with hidden dimension of 30 and 12, respectively. The model CFA-CON uses the same architecture as CON-M. We compare against several popular latent space model architectures: The
NODE model uses MLP with an hyperbolic activation functions and predicts { = INODE (roDDE((()).ud)) u(t)) To make the comparison fair, we parametrize the NODE's MLP in the same fashion as for CON-M -M.-M. The MECH-NODE integrates prior knowledge towards learning 2nd 2dd-order mechanical ODEs and, therefore. predicts a(t) MAECH-NNOEE (E (t), u(t)). Furthermore, we consider multiple autoregressive models; RNN, GRU, and coRNN and let them parameterize the following transition function: {(f++1) far(((t)), u(fx)) As common the relevant literature (26), we allow the autoregressive models to perform multiple time step transitions before predicting the next sample. For the autoencoder, we use vanilla Convolutional Neural Network (CNN). More details can be found in Appendix
Datasets. We consider in total six datasets that are based on simulations of unactuated mechanical systems. and actuated continuum soft robots. The first three, mechanical dataset are based on the work of Botev al. 26 and contain video sequences of a mass-spring system with friction (M-SP+F), single pendulum with friction (F)P+F), and a double pendulum with friction (D-P (D-P+F). Continuum soft robots have theoretically infinite Degree of Freedom (DOF), evolve with highly nonlinear and often time-dependent dynamical behaviors, and are
notoriously difficult to model from first principles [50]. For that reason, it is very interesting proposition if
we could learn latent-space dynamical models directly from video [51] and later leverage them for control [52) Therefore. we generate three datasets based on the Piecewise Constant Strain (PCS) soft robot model. CS considers one segment with constant strain and is modeled using three configuration variables. PCC-NS-2 and PCC-NS-3 only consider bending deformations and contain soft robots with two and three segments, respectively.
For all datasets, we render images with a resolution of 32x32px of the system's 's state. More information on the datasets can be found in Appendix C.T
We tune all hyperparameters for each model and dataset separately using Optuna 53
Results. Unactuated mechanical datasets: The results in Tab. show that the NODE model slightly outperforms the CON network on the M-SP+F and S-P+F datascts. However, as the datasets do not consider system inputs, we can remove the input mapping from all models (c.g., RNN, GRU, CoRNN, CON, and CFA- CON). With that adjustment, the CON network has the fewest parameters among all models, particularly two orders of magnitude less than the NODE model. Therefore, we find very impressive that the CON network roughly par with the NODE model. For the D-P+F dataset, we can conclude that the CFA-CON model offers the best performance across all methods. Finally, most of the time, the CON & CFA-CON networks outperform the other baseline methods that have more trainable parameters. Actuated continuum soft robot datasets: The results in Tab. [ show that CON-M matches the performance of the ofte-of-thhe-- methods
7

---

0.5
SIRCE-NO SIRCE-TODB
SIDCH- SODIESS
0.4
0.4
MECCI-NGDE
RIMSE
0.3
F3RS PON
RMSE
0.3
CPA. CON
0.1
2000 4000
6000
8000
10000
12000
Model parameters
(a) RMSE vs. latent dimension ns
(b) R.SSEE vs. model parameters
Figure 3: Evaluation prediction performance of the various models vs, the dimension of their latent representation n, and the number of trainable parameters the dynamics model, respectively, on the PCC-NS-2 dataset. All hyperparameters are tuned for each model separately for n. The error bar denotes the standard deviation across three random seeds.
across all experiments. In the case of PCC-NS-3, CON-M even decreases the RMSE error by % w.r.t. the closest baseline method (MECH-NODE), Impressively, the performance is not reduced (but instead often even improved) compared other models that offer much larger design space for learning the dynamics (c.g. NÓDE). Furthermore, CON-S and CFA-CON often only exhibit slightly lower performance than CON-M, even though they have significantly fewer parameters and consider approximated solution, respectively. Supplementary results (c.g.. more evaluation metrics) can be found in Appendix D We also conduct the PCC-NS-2 dataset analysis concerning the effect of the latent dimension the performance (see Fig. 3 For this experiment, all hyperparameters were tuned for n, 8, and we observe that the CON models have much-improved consistency and smaller variance w.r.t the baseline methods when the latent dimensionality is increased.
The results for a dataset based on the -orode reaction-diffusion PDE 9 are presented in Apx. D.I
5 Exploiting the dynamic structure for latent-space control
We consider the problem setting of guiding the system towards a desired observation of by providing a sequence of inputs u((x) such that at time iv. the actual observation olvv) matches o". A relatively simple way would be to encode the desired observation into latent space 2d (oo") and then to design a feedback controller (c.g., PID) in latent space: g(u) PID(2" z(t), Unfortunately, several challenges appear: first, is not clear how the latent-space force r = g(u) can mapped back an actual input (()) as the inverse input-to-forcet mapping g -1 is generally not known. Furthermore, relying entirely on PID controller has several well-knou known drawbacks, such as poor and slow transient behavior, stcady-state errors (in case the integral gain is chosen to be zero), and instability for high proportional and integral gains. We take inspiration from potential shaping strategies [23] [24], which are widely used for effectively controlling (elastic) robots, and, therefore, combine a feedforward term compensating the latent-space potential forces with an integral-saturated. PID-like feedback term. For mapping the desired forcing 7 back to an input u(!). we train an forcing decoder RR" - R
that approximates g Specifically, we consider here the structure u 17(T) F((T) 7, where E B parameterized by MLP.
The latent-space control law is given by
7(1) g(u)= Kn 2d + tanh(2" dl +)) + Kp (2 -2) Ka Ki lanh(u(=" =(()))))))())))))) ( (13)
Peodforwazd term: compensation potential
Feedback tena: P-sutt-1>
where Kp. Ki, KA E R are the proportional, integral, and derivative control gains, respectively. As integral terms can often lead to instability when applied to nonlinear systems [54]. we adopt an integral term saturation (55) with the associated dimensionless gain v E R. which ensures that the integral error added at each time step is bounded to the interval (-1, 1). Subsequently, 7 is decoded to the input as u(t) = n)() =
E((T) r. For training this decnder, we add a reconstruction loss to (12) L. (tx) = > MSF((att)) ū(ta)) An, MSE (n(tx), 7 (((((t(x))) ())))))).
Experimental setup. We train CON model with two latent variables (n, = 2)
PCC-NS-2 dataset. Analog to the input encoder mapping Buu), the forcing decoder mapping E(r) is parametrized by an MLP consisting of five layers with hidden dimension 30 and hyperbolic tangent activation function. The CON model achieves RMSE of (. 1626 on the test set. We benchmark two controllers on the simulated continuum soft robot consisting of two segments: (i) a pure P-sail-D controller (i.e.. the feedback term in (13) that leverages
the smooth mapping into the latent representation enabled by the CON dynamic model and the 3-VAE, and (ii)
a P-satl D FF (i.e.. (13 that exploits the structure of the CON dynamics by compensating for the potential forces. The stable closed-loop system dynamics made the control gain tuning very easy, and selected
8

---

Beal-sorrld Inspiration
Rondaood Datasst Emsggy
Potential shaping Tit Ky tambls"
Physical system
Target obsrrvation
Yowee-pring + trutton
(35-P) 13
Encoder
Latent-space forcing decoder Blyl
Double Peodiiu Fiieiio (DP+P
Current observation
Encoder
Continum Soft Robor (PCC-NS-))
Ate
Rowtim Diffuria (8-D3
P-snt L-D feedbock
-iig (2 - a) Na tmmbiots" =(()))
(a) Samples of used datasets
(b) Blockscheme of model-based control in latent space
Figure 4: Panel (a): Samples some datasets used as part the experimental verification, specifically for the results reported in Tab. [ The real-world Seaction-Diffusiio image is adopted from [43 Panel (b): Model-based control latent space by exploiting the physical structure of the CON model.
Kp 1, K; = 2, Ka 0.02, e = 1 for the P-satl-D controller and K 0, K; 2, Ka 0.05, v = for the P-satI-D+FF controller, respectively.
Furthermore, we compare the control performance of our model-based controllers with baseline control strategy based on the MECH-NODE (n, P-satI-D feedback controller as for the CON model generate the control action 7(t) in latent space. As the
that achieves error of 0.1104 the test set. First, we utilize the same
MECH-NODE uses MLP to parameterize the function & = fc (E, u), cannot easily map r(t) into input (t).. Therefore, we linearize the latent space dynamics w.r.l to the input as (E, u) = fc (E, + A(E) where A(E) afs 3u (g, 0) is computed using autodiff. Then, u(f) A () =(t). After tuning the control gains, we choose Kp 0.001, Ki 0.02, Ka le-5,
We train all models (c.g., NODE,NODE, CON) three different random seed and choose the best model instance. Please refer to Appendix D for more details the model selection. We generate trajectory of 7 setpoints, where z" (tz) ~ (-5m, 5)) rad/m m E R2 is a sampled configuration of the soft robot. Then, we render an image o (4,) that represents the target observation for the controller and encode it into latent space to retrieve 2d E R At time step k, we render an image o(tx) of the robot's 's current configuration g(ts) and encode the image. Subsequently, we evaluate the control law and apply the decoder u((4) = n (1x ))}. which is finally passed to the simulator that integrates the ground-truth dynamics to the next time-step f&+1 considering the actuation ). ). The controller runs at 100 Hz, and we simulate the ground tround-truth dynamics with a Dopriš ODE integrator at a time-step of le- s.
Results. As an evaluation metric, we consider the RMSE between the actual and the reference trajectory. The P-sarl-D applied to the MECH-NODE model (baseline) achieves an RMSE of 2.88 rad/ / w.r.t. to the desired configuration d (but unknown to the algorithm). The P-satl-D CON controller. which does not exploit the learned latent dynamics for control, exhibits an RMSE of 4.08 rad/ m w.r.t. to the desired configuration q" The P-satl-D++F controller exhibits an RMSE of 2.12 2.12 m w.r.t. to desired configuration qd We also visualize the closed-loop trajectories Fig. and as sequences stills in Apx. structured latent space generated by the 3-VAE allows the P-satl-D controller effectively regulate the system towards the setpoint. although the response time rather slow. The P-satl-D4 FD controller is able exploit the structure of the CON model through its potential shaping feedforward term. With that, CON P-satl D+ exhibits a faster response time and
conclude that the nicely
lower RMSE than the MECH-NODE P-tat-D baseline. We provide results for the control of actuated damped harmonic oscillator Apx.
6 Conclusion and Limitations
Conclusion. In this work, we propose new formulation a coupled oscillator that inherently input-to- state stable. Additionally, we identify a Elosed-form approximation, that able simulate the network dynamics more accurately compared numerical ODE integrators with similar computational costs. When learning latent dynamics with CON, we observe that the performance on par or slightly better compared SoA methods such RNNs, NODES, etc., even though constrained the solution space ISS-stable coupled oscillator structure. Funthermore, we point out that the performance of the CON models is compared to the baselines and improved when not specifically tuned for a given dimension. Furthermore, as number of model parameters, similar
consistent across latent dimensions
seen in Tab. 8 the closed-form approximation achieves, with the accuracies and double the training speed w.r.t. the continuous time model. Finally, demonstrate that even simple PID-like patent-space controller effectively regulate the system setpoint. By exploiting the

---

10
Confgitaaion (nvd)
welak::
22 2L 3
-0.02
4
-0.00
-0.00
(a) Config. for P-satl-D D with MIECH-NODEE
(b) Latent for P-satl-D with MECH-NODE
10
Conggiaaton lookha)
0.0
-02
-05
16
0
10
15
20
55
02
15
20
33
55
Tine
(c) Configuration for P-satI-D with CON.
(d) Latent for P-satI-D with CON.
20
Configattin (ndd/nn)
10
0.2
0.0
-0.6
10
Time 3
Time
(e) Configuration for P-satI-D+FF with CON.
(I) Latent for P-satI-D+FF with CON.
Figure 5: Latent-space control of continuum soft robot (simulated using two piecewise constant curvature segments) at following sequence of setpoints: The upper two rows show the performance of a pure P-satI-D feedback controller operating in latent space learned with the MECH-NODS NODE and CON models, respectively. The lower row displays the results for a latent space controller based on the CON model that additionally also compensates for the learned potential forces.
network structure and compensating for potential forces, regulation performance can be greatly improved, and response time decreased by more than 55 %.
Limitations. While we think our proposed method shows great potential and opens interesting avenues for future research, there exist certain limitations. For example, the proposed method of learning (latent) dynamics implicitly assumes that the underlying system adheres the Markov property (c.g., the full of the system observable), that system with mechanical structure asymptotically stable equilibrium. This is, for example, the continuum soft robots, deformable objects, and elastic structures) with continuous dynamics, convex elastic behavior. dissipation, and whose demendepn effects (e.g.. viscoelasticity, hysteresis) negligible. Even if these conditions are not met globally, the method can be applied model the local bchavior around an asymptotic cquilibrium point of the system (c.g., robotic manipulators, legged robots) with added stability benefits for -iitof-distriuution samples. Alternatively, the method could be extended to relax some of these assumptions, e.g., by allowing for multiple equilibria, zero damping, or by incorporating additional terms to capture discontinuous dynamies (e.g.. stick-slip models) or period motions (e.g.. limit cycles such as the Van der Pol oscillator). The proposed method might not be suitable for some physical systems, such nonholonomic systems, partially observable systems, or systems with non non-Mart Markovian properties. Examples of such systems include mobile robots and systems with hidden states or delayed observations.
approximate it. and that it has isolated. globally
for many mechanical systems (e.g.. some
Furthermore, the approximated closed- fored-form solution shows best integration situations where linear, decoupled dynamics dominate the transient. For dominant nonlinear, coupled forces, the performance CFA- CON degrades, and it might better to revert numerical integration the CON ODE. Finally, control works exceptionally well in the setting where the latent dimension equals input dimension. We hypothesize that this enables the method to identify diffeomorphism between the input and the latent-space forcing. Still not investigated
performance could degrade if nts > m (or né m that matter).
10

---

Acknowledgments and Disclosure of Funding
This work was supported under the European Union's Horizon Europe Program from Project EMERGE Grant Agreement No. 101070918.
References
[1] David Ha and Jürgen Schmidhuber. World models. arXiv preprint arXiv: 1803. 100p./0033.022 10/22, 2018.
[2] Danijar Hafner, Jurgis Pasukonis. Jimmy Ba, and Timothy Lillicrap. Mastering diverse domains through world models. arXiv preprint arXiv:2301.04 04.044 2023. [
[3] Yutaka Matsuo, Yamn LeCun, Maneesh Sahani, Doina Precup. David Silver, Masashi Sugiyama, Eiji Uchibe, and Jun Morimoto. Deep learning, reinforcement learning. and world models. Neural Networks,
152:267-275, 2022. |
[4] Cheng-Yuan Liou, Wei-Chen Cheng, Viun-Wei Liou, and Daw-Ran Liou, Autoencoder for words. Neuro- computing, 139:84-96 96, 2014.
[5] Dicderik P Kingma and Max Welling. Auto-encoding variational bayes. In Proceedings of the International Conference on Learning Representations, 2014.
[6] Manuel Watter, Jost Springenberg, Joschka Boedecker, and Martin Riedmiller. Embed to control: A locally linear latent dynamics model for control from raw images. Advances in neural information processing systems, 28, 2015. 2
[7] Ian Lenz, Ross A Knepper, and Ashutosh Saxena. Deepmpe: Learning deep latent features for model predictive control. In Robotics: Science and Systems, volume 10, page 25, Rome, Italy, 2015. !
[8] Niklas Wahlström, Thomas B Schön, and Marc Peter Deisenroth. Learning deep dynamical models from image pixels. IFAC-PPpersOOmLinee 1855514405,,1044 2015.
[9] Kathleen Champion, Bethany Lusch, Nathan Kutz, and Steven I. Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45):22445-
22451, 2019.
[10] Yaofeng Desmond Zhong and Naomi Leonard, Unsupervised learning of lagrangian dynamics from images for prediction and control. Advances in Neural Information Processing Systems, 33:10741-10752,, 2020.
[11] Herke Van Hoof, Nutan Chen, Maximilian Karl, Patrick van der Smagt, and Jan Peters. Stable reinforcement learning with autoencoders for tactile and visual data. In 2016 IEEE/RSJ international conference on intelligent robots and systems (IROS), pages 3928-3934. IEEE, 2016. D
[12] Carles Gelada, Saurabh Kumar, Jacob Buckman, Ofir Nachum, and Mare G Bellemare. Deepmdp: Learning continuous latent space models for representation learning. In International conference on machine learning, pages 2170-217220 2170- PMLR, 2019. []
[13] Danijar Hafner, Timothy Lillicrap. Jimmy Ba, and Mohammad Norouzi. Dream to control: Learning behaviors by latent imagination. In International Conference on Learning Representations, 2019. ||
[14] Wilko Sehwarting, Tim Seyde, Igor Gilitschenski, Lucas Liebenwein, Ryan Sander, Sertac Karaman, and Danicla Rus. Deep latent competition: Learning to race using visual control policies in latent space. In Conference on Robot Learning, pages 1855-1870. PMLR, 2021.
[15] Danijar Hafner. Timothy Lillicrap. Ian Fischer. Ruben Villegas, David Ha, Honglak Lee, and James Davidson. Learning latent dynamics for planning from pixels. In International conference on machine learning, pages 2555-2565. PMLR, 2019. [
[16] Lukas Hewing. Kim P Wabersich, Marcel Menner, and Melanie N Zeilinger. Learning- based model predictive control: Toward safe learning in control. Anmual Review of Control, Robotics, and Autonomous
Systems, 3:269-2966 2020. [
[17] John Irvin Alora, Luis A Pabon, Johannes Kühler, Mattia Cenedese, Ed Schmerling, Melanie N Zeilinger. George Haller, and Marco Pavone. Robust nonlinear reduced-order model predictive control. In 2023 62nd IEEÉ Conference on Decision and Control (CDC), pages 4798-4805. IEEE, 2023. [
[18] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, and J Nathan Kutz. Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control. PloS one, 500))70,
2016. 0
[19] Giorgos Mamakoukas, Maria Castano, Xiaobo Tan, and Todd Murphey. Local koopman operators for data-driven control of robotic systems. In Robotics: science and systems, 2019. [
[20] David A Haggerty. Michael J Banks, Ervin Kamenar, Alan B Cao, Patrick C Curtis, Igor Mezić, and Elliot W Hawkes. Control of soft robots with inertial dynamics. Science robotics, 8(81) x681)cadd6864, 2023. [
[21] Jun Yamada, Chia-Man Hung, Jack Collins, loannis Havoutis, and Ingmar Posner. Leveraging scene embeddings for gradient-base based motion planning in latent space. In 2023 IEEE International Conference on Robotics and Automation (ICRA), pages 5674-5680 IEEE, 2023. [
11

---

[22]
Mattia Cenedese. Joar Axás, Bastian Bäuerlein, Kerstin Avila, and George Haller. Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds. Nature communications, 13(1))872,
2022. 2
[23]
Anthony M Bloch, Dong Eui Chang, Naomi Ehrich Leonard, and Jerrold E Marsden. Controlled lagrangians and the stabilization of mechanical systems. ii potential shaping. IEEE Transactions on Automatic Control,
(6(10):1556-157) 2001.
[24] Romeo Ortega, José Guadalupe Romero, Pablo Borja, and Alejandro Donaire. PID passivity- control
of nonlinear systems with applications. John Wiley & Sons, 2021. 2
[25] Marco Lepri, Davide Bacciu, and Cosimo Della Santina. Neural autoencoder-based structure-preserving model order reduction and control design for high-dimensional physical systems. IEFE Control Systems Letters, 2023. 2
[26] Aleksandar Botev, Andrew Jacgle, Peter Wimsserger, Daniel Hennes, and Irina Higgins. Which priors mat- ter? benchmarking models for learning latent dynamics. In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 1). 2021.
49
[27] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018. 2
30 31
[28] Aleksei Sholokhov, Vuying Liu, Hassan Mansour, and Salch Nabi. Physics-informed neural ode (pinode); embedding physics into models using collocation points. Scienrific Reports, 1311), 2023. 2
[29] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):1735 1780, 1997. 2
[30] Kyunghyun Cho. Bart Van Merriënboer. Caglar Gulcehre. Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using encoder-decoder for statistical machine translation. arXiv preprint arXiv: 1406.1078. 1078, 2014.
2
7
26
29
30
31
[31] Michael Lutter. Christian Ritter, and Jan Peters. Deep lagrangian networks: Using physics as model prior
for deep learning. In International Conference on Learning Representations. 2018. 2
[32] Miles Cranmer, Sam Greydamus. Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. La- grangian neural networks. arXiv preprint arXiv:2003.04633 2020. P
[33] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks, Advances in neural
information processing systems, 32, 2019. 2
[34] Alistair White, Niki Kilhertus, Maximilian Gelbrecht, and Niklas Boers. Stabilized neural differential
equations for learning dynamics with explicit constraints. Advances Neural Information Processing
Systems, 36: 16:99-9112 2023. -
[35] T Konstantin Rusch and Siddhartha Mishra. Coupled oscillatory recurrent neural network (cornn): An accurate and (gradient) stable architecture for learning long time dependencics. In International Conference on Learning Representations, 2020.
3
4
31
[36] T Konstantin Rusch and Siddhartha Mishra. Unicornn: A recurrent model for learning very long time
dependencies. In International Conference on Machine Learning, pages 9168-9178. PMLR, 2021.
[37] Andrea Ceni, Andrea Cossu, Maximilian Stölzle, Jingyue Liu, Cosimo Della Santina, Davide Bacciu, and Claudio Gallicchio. Random oscillators network for time series processing. In International Conference Artificial Intelligence and Statistics, pages 1807-4815 PMLR, 2024. 2
[38] Samuel Lamthaler, T Konstantin Rusch, and Siddhartha Mishra. Neural oscillators are universal. Advances in Neural Information Processing Systems, 36, 2024. 3
[39] T Konstantin Rusch and Daniela Rus. Oscillatory state-space models. arXiv preprint ar arX0v:2410003943,
2024. 2
[40]
Hassan K Khalil. Nonlinear systems third edition. Patience Hall, 115, 2002.
1
[41] John Irvin Alora, Mattia Cenedese, Edward Schmerling, George Haller, and Marco Pavone. Data-driven spectral submanifold reduction for nonlincar optimal control of high-dimensional robots. In 2023 IEEE International Conference on Robotics and Automation (ICRA), pages 2627-2633. IEEE, 2023. [
[42] Maximilian Stölzle and Cosimo Della Santina. iiston-driven pneumatically- actuated soft robots: Modeling and backstepping control. IEF. Control Systems Letters, 6: 6:183711842, 2021. 37
[43] Irving R Epstein and Bing Xu. Sraction-diffusion processes at the nano-and microscales. Nature nanotechnology. 11(4):312-31 119,,,19222201 2016. 3
[44] Davide Calzolari, Cosimo Della Santina. and Alin Schäffer Exponential convergence rates of nonlinear mechanical systems: The 1-dof case with configuration--eeeeneet inertia. IEEE Control Systems
Letters, 5(2):445-455 4500, 2020. 4
[45] Xuwei Wu. Christian Ott, Alin Albu-Schäffer, and Alexander Dietrich. Passive decoupled multitask controller for redundant robots. IEEE Transactions on Control Systems Technology. 31(1):1-16, 2022. 45
[46] Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013. 5
[6
18
[47] P Kidger. On neural differential equations. PhD thesis. University of Oxford, 2021.
6
27
12

---

[48] Pieter Pas. The damped harmonic oscillator, Mar 2023.
6 20
[49] Irina Higgins, Loic Matthey, Arka Pal. Christopher P Burgess. Xavier Glorot, Matthew M Botvinick. Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with constrained variational framework. ICI. (Poster), 3, 2017. [
[50] Costanza Armanini, Frédéric Boyer, Anup Teejo Mathow, Christian Duriez, and Federico Renda. Soft
robots modeling: A structured overview. IEEE Transactions on Robotics, 39(3):1728-1744 2023. |
[51] Thomas George Thuruthel and Fumiya lida. Multi- sensor fusion for learning rich models for interacting soft robots. In 2023 IEEE International Conference on Soft Roboticx (RoboSoft), pages 1-6.
IEEE, 2023. 7
[52] Elijah Almanzor, Fan Ye, Jialei Shi, Thomas George Thuruthel, Helge A Wurdemann. and Fumiya Tida. Static shape control of soft continuum robots using deep visual inverse kinematic models. IEÉE Transactions on Robotics, 2023. 7
[53] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A
next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pages 2623 1631, 2019. 7 27
[54] Maximilian Stölzle, Danicla Rus, and Cosimo Della Santina. An experimental study of model-based control for planar handed shearing auxcties robots. Experimental Robotics, pages 153-167, Cham, 2024. Springer Nature Switzerland. 24
[55] Pietro Pustina, Pablo Borja, Cosimo Della Santina, and Alessandro De Luca. P-sati-d shape regulation of soft robots. IEEE Robotics and Automation Letters, 8(1):1-8, 2022. g
[56] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
17 18
[57] Mattia Gazzola, IHH Dudte, AG McCormick, and Lakshminarayanan Mahadevan. Forward and inverse problems in the mechanics of soft filaments. Royal Sociery open science, 5(6):171628, 2018.
24
[58] Federico Renda, Vito Cacucciolo, Jorge Dias, and Lakmal Seneviratne. Discrete cosserat approach for soft robot dynamics: A new piece-wise constant strain model with torsion and shears. In 2016 IEEE/RS/ International Conference on Intelligent Robots and Systems (IROS), pages 5495-5502. IEEE, 2016. 24
[59] Cosimo Della Santina, Christian Duriez, and Daniela Rus. Model-based control of soft robots: A survey of the state the art and open challenges. IEEE Control Systems Magazine, 43(3):30 65, 2023. 24
[60] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclau- rin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Vilne, and Qiao Zhang. JAX:
composable transformations of Python+NumPy programs, 2018. [4] 27
[61] Maximilian Stölzle, Sonal Santosh Baberwal, Daniela Rus, Shirley Coyle, and Cosimo Della Santina. Guid- ing soft robots with mooorr-imagrry brain signals and impedance control. In 2024 IEEE 7th International Conference on Soft Robotics (RoboSoft), pages 276-283. IEEE, 2024. 24
[62] G. Bradski. The OpenCV Library. Dr. Dobb's 's Journal of Software Tools, 2000. 24
[63] Robert J Webster III and Bryan A Jones. Design and kinematic modeling of constant curvature continuum robots: A review. The International Journal of Robotics Research, (133):66611168 2010. 24
[64] Emanuele Riccardo Rosi, Maximilian Stölzle, Fabio Solari, and Cosimo Della Santina. Sensing soft robots' shape with cameras: an investigation on kinematics aware slam. In 2022 IEEE 5th International
Conference on Soft Robotics (RoboSoft), pages 795 801. JEEE, 2022. 24
[65] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization arXiv preprint 1eXivv1607.0444 2016. 25
[66] Vincent Dumoulin and Francesco Visin. A guide convolution arithmetic for deep learning. arXiv preprint wrXiv:1603.0728. 2016. 25
[67| Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7115)5110, 2008. 25
[68] Jonathan Heek, Anselm Levskaya, Avital Oliver, Marvin Ritter, Bertrand Rondepierre, Andreas Steiner, and Marc van 7ec. Flax: A neural network library and ecosystem for JAX, 2023. 27
[69] R. Dormand and J. Prince. A family of embedded Runge-Kutta formulae. J. Comp. Appl. Math, 6:19 26, 1980. 27
[70] Ilya Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with warm restarts. In International Conference on Learning Representations, 2016. 27
[71] Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv: I4Xiv:1112.6988 1412.6980, 2014. 27
[72] Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Conference on Learning Represemtations, 2018. 27
[73] Shuhei Watanabe. Tree- parzen estimator: Understanding its algorithm components and their roles for better empirical performance. arXiv preprint arXiv 2023. 27
13

---

[74] Liam Li, Kevin Jamieson, Afshin Rostamizadeh, Ekaterina Gonina, Jonathan Ben-Tzur, Moritz Hardt, Ben- jamin Recht, and Ameet Talwalkar. A system for massively parallel hyperparameter tuning. Proceedings of Machine Learning and Systems, 2:230 246. 2020. [77
[75] Guilin Liu, Fitsum A Reda, Kevin Shih, Ting- Ting-Chun Wang, Andrew Tao, and Bryan Catanzaro. Image inpainting for irregular holes using partial convolutions. In Proceedings of the European Conference on Computer Vision (ECCV), pages 85 100, 2018. 27
[76] Jiahui Yu, Zhe Lin, Jimei Yang, Xiaohui Shen, Xin Lu, and Thomas S Huang. Generative image inpainting with contextual attention. Proceedings of the IEEE conference on computer vision and pattern recogmition, pages 5505-5514, 2018. 27
[77] Maximilian Stölzle, Takahiro Miki, Levin Gerdes, Martin Azkarate, and Marco Hutter. Reconstructing oc-
cluded elevation information terrain maps with self-supervised learning. IEEE Robotics Automation Letters, 7(2):1997-1706 2022. 27
[78] Zhou Wang. Alan Bovik, Hamid R Sheikh, and Eero Simoncelli. Image quality assessment: from error visibility to structural similarity. JEEE transactions on image processing, 13(4)x600-61 2004, 27
[79] DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju. Jake Bruce, Peter Buchlovsky. David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Claudio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Mattco Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Miloš Stanojević, Wojciech Stokowiec, Luyu Wang, Guangyao Zhou, and Fabio Viola. The DeepMind JAX Ecosystem, 2020. 28
14

---

A Appendix on proof of Input-to-State Stability (ISS)
the following, we will write ||A] to denote the induced norm of matrix and (A), (A)(A) to refer to its minimum and maximum Eigenvalue respectively.
(f)]l:
sup [lo(u(e))] Musegg
2.0
(to)l. to) ExC
(soll 40) EKC
0.5
0.0
-1.0
-15
-15
(a) GAS: g(u)
(b)
Figure 6: Illustration of global asymptotic stability for the unforced system with g(u) and input-to-state stability for the forced system, where the black dashed line denotes the input u(t).
We introduce some expressions often used throughout this section: The gradient of Vy, w.r.t. the residual coordinate Tw given
IV
(aanh (Jw + dw + b)
tanhh(Tx +
(ihw Pv iw +
(14)
0
0"
Next, the Hessian of the Lyapunov candidate can be derived as
Hy (v(
8VV₂ Oya
[K + S2 (āw)
1 MM Mr
E 222n22n
(15)
M
where
Szech (#)) = diag(sech( (Tw + Tw + b)) E Raxn >
Viw E R..
(16)
Furthermore, the Schur complement of Pv is given by
Spy My 12 MM x Kw Mw
(17)
A.1 Potential force and energy
Lemma 3. Let Jw IR" be generalized coordinates and Rw. b E I2" constants. Then, the potential force
system 4
fun (āw) = Kw (āw + Rww + tanh(Bw + Iw + b),
(18)
stems from the potential
Um (5w) = Σ /" tanh( to +b+ do - - E/mm tanh(#... b.) R.
(19)
Proof. First, we take the derivative of (((w):
(w + zow tanh(ā + for + b) fur.
(20)
83w
The Hessian of the potential given by
8u Hu. (ā..) = di3
Kw + Sauch E nxn Diw
(21)
As Kw x 0 + Kw KT, we can easily show that the potential force is symmetric:
Hi: K.t + (Zx)) Kw + (āx) Hu
(22)
L
15

---

A.2 Proof of Lemmal I single & isolated equilibrium
Lemma restated. 5w
Let Kw > 0. Then, the dynamics defined in have single, isolated equilibrium 07]².
Proof. We regard the characteristic equation function: beg (xw) tanh(a + b) + Kw 2w. For there exist multiple equilibria, heg (3w) 0 would need be true for multiple However, take the partial derivative of (Tw) w.r.t 2w and see that
Oheg Kw + Sunch(2in O, VIw R"
with Sarch (Iw) diag(sech(Ew + 6]) ==**
(23)
Swach (xw) > O the zero line once.
Vie E IR" and Kw > O. Therefore, loag (Fw) is continuously increasing and can only cross
A.3 Proof of Lemma Validity of strict Lyapunov candidate (ūu(w)
Lemma 4. Suppose Mw Kw and p positive definite.
y (M ( Amilww M|
pvv. Then Spy, as defined
Proof. The minimum Eigenvalue of Sny is bounded by
Am (Spr) Am (m((Ww n²² I Kw Mwll.
||
(24)
Am (Mw) - n²² Aw (Kw)
Based on the assumption Mw > 0, Kw O. we can state
[| (Kw)
> Therefore, the critical case for the
lower bound on Am (Spx) is μ Am [|
Am(
pV. Hence,
Xm (Sp)) > (M Au (Mw) Am (Kw)
1l Aw (Kw)
(25)
|M 1|
Consequently, the Eigenvalue sensitivity theorem (46 demands that Spy > 0.
Lemma 5. Ler M. > 0, Kw > O, and a < v Am Amm(M)) Am = pay. Then, it follows that Py > and Hv(āw) > Sw E RD.
[M
Proof. By inspecting the expressions for Pv > and Hy(āw) > 0 Y E TR" in Equations 5 and 15 respectively, it can be easily seen that Hy (Tw) > Py Viw E R. As Lemma Astttes that the Schur complement of By is positive definite, follows that (Iw) 2 Py >
Lemma 6. Suppose Tw, Rw, b 1tt and n E N Then,
( E/ tanh(āw., +o bb) Σ Aw..
tanh(āw.. do
(26)
positive semi-definite function.
Proof. Proving Iv (. is equivalent showing that the scalar function Ry (v, tanh(g a) do - f tanh(a) do 2 0 V R, where we set = Bo.s and Rw.a + b..
We strive to find the critical points (i.c., minimas and maximas) F of (r) and, for this, analyze where the first derivative of hv.u. (r)) is zero
Obe ar
(F) tanh(F + a) tanh(a) =0,
(27)
which is case only 0. Next, we compute the second derivative at as
oev
sech (r) 1.
(28)
Thus, (r) is convex and its global minimum F takes the value hy, (0) (0) -. As result, loy (5w) is also positive semi-definite.
16

---

A.3.1 Proof of Lemma
Lemma 2 restated. The scalar function Vu(Uw) defined in ( is continuously differentiable and verifies
condition (0) Furthermore, let Mw, Kw Now, if we choose < μ < Amu (dkw) Aua )ml(K)) Jvv.
Me then (ūw)) > 44 Dw E R \ (0), Additionally, then (ūw) is radially unbounded as 11 + 00 7
Vu(ūw) + cc.
Proof. Step 1: It can be easily that V (ūw) in 5 smooth and continuously differentiable.
Step 2: Proof that (0) O.
(0) = 0 Σ / (Tw, + t b,) do Σ / tanh(Tw, + b:) do = O.
a
(29)
Step 3: Proof that the Lyapunov candidate is positive definite; i.e., (ūw) > 0 Virw E ER" [0].
As the gradient of the Lyapunov candidate, as defined (14). is zero for Ny
0V, (0)
taah( + b) tanh + b)
0 = -0
(30)
0iw
0
Tw = 0 critical point of (pw.) According Lemma 5 the Hessian in 113 is positive-definite [56]: Hv(ww) > 0 Vÿw E 2 With that,
5 is convex and its global minimum is at yo 0, where Vú₂ = 0. In
summary, we state Vi, > Vgw
{0).
Step Proof that Lyapunov candidate radially unbounded: i.e., |/w + ac ((( + oc. Lemma exploited for identifying a lower bound on V₁, (ūw):
Bwn
(ū((w)
T Vo
Py Sw
E/
tanh(Tw.a to + b.) do
tanh(5w.s + bi) do,
m=1
(31)
T Uw
Pv ju
1
Am (Pv) 11
2
Lemma 5 tells us Pv > and with that Am Am() > 0. Therefore, if [|@ + oo, it also follows that (ūw) 0c.
A.4 Proof of Theorem 1 Global asymptotic stability of unforced network
The Lemmas introduced below are used to prove Theorem T
Lemma 7, Suppose Bo Tw; E R" and n E N Then, the function by sb. (Tw) defined as
(Tu) = (tanh(Ew + 2n + b) tanh(ā + b)) Zw:
(32)
positive semi-defiiit.
Proof. Proving by (
(āw)
0 equivalent to proving that the scalar function (r rra E R, where we Sw.i and Sw.i b.. Now,
(tanh(r a) tan)(a))
pand the hyperbolic tangent:
Av (r) = (tanh(r + a) tanh(a)) r= ((+4))
2 -
(l(++)) +1
co +
r,
(33)
20 (62r -1)
2 0,
22r+40 + 2r+ + 20
as the denominator e2 2r+4a + 2r + +1> Vr IR and as sign (2e2n (e2- - 11) = sigg(r)). For
example, 2r -120V: 2Vr Analog, 2r 1<00
Lemma 8. Let Ma > Kw > 0, and Dw > Also, let p E R+ be chosen such that
< P
Aua (Dw.)
12 == pv- Then, the matrix Pv = L (Am((W)
μKw
1/D Dw WMm
E B is positive definite.
Am (Mdw) +
Proof. The Schur complement of Py is given by
1
Spy = Dw pM
uDD K= Dw.
(34)
17

---

The lower bound on the smallest Eigenvalue of Spy can be identified as
D
Am (Spr) 2 Am (Dw) - (M
(35)
Ama (mr(Kw)
As Kw, Dw > O, we know that 1.Dw 0. Therefore, the case μ the lower bound on Am (Sp
Amn (Dw)
Hig determines
[m(Kw]
Aui (w
4am
12x
1D
Am (Spp) > Am Am HV (Mm)
O.
(36)
4 Am (mKw))
We conclude, based on the Eigenvalue sensitivity theorem of symmetric matrices that Pv. [56].
[46] that Sp. > 0 and with
A.5 Proof of Theorem Proof of Input-to-State Stability (ISS)
Lemma 9. Let Rw Rw, b 1e". Then,
(v =
tanh( to + b.) do tanh(Tw. + bi) do S 2 [wwl- (37)
Proof. Proving un S [ww] is equivalent proving that scalar function Ay (r))))) = Sh" tanli(o + a) do So lanh((a) do 2/r Vr, R, where set Zw,s and Zw,t + b.. We perform the integration contained in (r)::(r):
(ry I." Lanli(o + a) do - ''" Lanli(a) da,
(38)
(r) log(coshi(r a)) log(cosh((a) tanh(a)
Next, demonstrate that the slope of |r| always larger than the magnitude of the slope of hy,
ttanh(r + a) tanh(a)) 2 Ör (2
(39)
Ör
Additionally, iv = | O. We conclude that (r) & 2/r| V E IR and with that, lov (Tw) 2 [iw] V Jw E R".
Lemma 10. Let Mw > 0 and K O. Then, (5) is bounded by the two scalar, class Koe functions an(7) = Am (Pv) and un (r) = AM vnn: on (11@ S Va(ūw) az((] 11)2
we already showed that (w)) is Lyapunov candidate. Now, we additionally also Proof. With Lemma 2
verify the conditions for ISS-Lyapunov candidates (4)]
Establishing bounds Vú, (y)w). We first identify lower bound ((x) by leveraging Lemma
Step 1:
Bw.s tanh( + b:) do,
Vū,( Uur Pv Tw
Σ/ mm
tanh(Tw, to +
2 Vu
-T
Pv Ju hv,))(
(40)
T Pv Ju 2 Am (Pv) 11 a1(1)=1l2).
1
2 V
2
Similarly, we derive an upper bound for (ūw) exploiting Lemma
Bw.s tanh(Tms + bi) do
(in)
T 2 1e
Py Vn
tanh(T to b.) do
1
1 Pv) |1 11 + || 5 3 XM((Pv) 11 / 12lll
(41)
<
2
Pv) 11 2 11 an((lūm1122)
Step 2: Proof that (r) (r) belong to class Kize-
According Lemma 5 Pv > 0 and with that Am. (Pv) O. First, analyze the behavior of (r): as it is belongs to class K. Furthermore, Therefore, an K... 40 (rz(r) is also strictly increasing for E [0, so), az()) 0, and it radially
strictly increasing 01 (0)
evaluate lim (1( co.
unbounded lim
For that reason, a2(r)) E Noo well.
L
18

---

A.5.1 Proof of Theorem 2
Theorem restated. Suppose Ma, Kw, Dw, O, < 0 1, and that we choose 0 < / < min(uy,
Then, is globally Input-to-Stat Stable (ISS) such that the solution Bw (1) verifies
ll lla S (110m(()))22: to) +
sup [lo(u(T)]la). MOSTEL
V to
(42)
where B(r,t) E K.C. 7(r) (1 (m() E K. 84
Proof. Step 1: Bounds on ISS-Lyapunov candidate.
Lemma 10 provides the Xo functions on (r) = 1 Am Am (Pv) and az(r) = + Vur such that
or @1(1
12) < (ūw))) I 02 (2(1@..1))
Step 2: Minimum energy dissipation.
Let < f minffv (annfxee)))) [4:]
proof Theorem We compute the input-dependent time-derivative repeat the derivations already made as part 7 (e.g.. exploiting
ISS Lyapunov candidate.
Lemmas 7 Jand S).
Un P& DN (tanh(æw Tw 6) tanh (#w + +))² 2w + FV [ug(u(t)] [(((tt))
Va u(t))
[ug(u(t))
S - Am (Pv) 11 +
Borr
[((((t))
(43)
((((1)))]
S Am (Pv) 11 100ll2
g(u(t))
Am (Py) || 1 /I+ u² (1@
where leveraged Hölder's inequality. We choose 0 such that < As consequence,
#2
(ww u(()) < -(1 0) Am (Pv) |1 11.
A llūwl12 2
[9(u(())l2 0.
(44)
0 Am (Py)
We define
2 0 Xm (Pv)
(3(r) (1 0) Am (Py) 22 and p(r)
(45)
Lemma 8 shows that Am (Py) Therefore, ag(r) is continuous positive function. Furthermore, μ
p(r) strictly increasing for r E co). Additionally with p(0) verified, can be stated that p(r) belongs to class [40]. We conclude that
Valu,, u(t)) -03 (11=11), V1ū 2 (Ilo(a(t)) 0.
(46)
Step 3: Conclusions.
As a result of Steps and 2, the system input-to-state stable, and with that, the solution je satisfies 40
1/w S 3 (1) (10)112, - 10) + [ls(u(())lla)
sup 105
(47)
with
(1 + n2)) AM (Pv) 40 n 1 n² Am (Py)r
7(r) = aj
02 p(r)
(48)
0 Am (Pv) X2. (Pv)
Indeed, based on Theorem [and and the associated proof. can easily verify that ((r) strietly increasing for [0, oo) and that 70)) As consequence, 7(r) E K. [40]
19

---

B Appendix on an approximate closed-fore solution for coupled oscillator networks
Kvi
Ft
dj
mi
Di
(a) Parameters of harmonic oscillator
1.0
AA
0.5
0.3
0.1
00.0.0)
(4:4)
0.000)
(A: (1.0.)
10
(0.0.0)
(R.
0..01))
-10
(b) Time evolution of a harmonic oscillator
(c) Damping regimes of a harmonic oscillator
Figure 7: Panel 7(a) Parameters of a forced and damped harmonic oscillator: my E RR+ denotes the mass,
RR the stiffness, and k, IRt the damping coefficient. The position and velocity of the oscillator measured (t) E and & (t) E IR, respectively. The oscillator be excited by the (potentially time-varying) external forcing (t) E TR. Panel 7bb) Time evolution of ID harmonic oscillator for different
values of Müs d, all in the undamped harmonic oscillator: undamped (C = 0). underdamped (0 <<< 1), critically damped (C = 1), and overdamped (zeta > 1).
underdamped regime. Panel 7(e): The four damping regimes of
B.1 Closed-form solution to a forced harmonic oscillator
As introduced in we consider the linear dynamics of a ID forced harmonic oscillator with state
[a 2.] R
dag
tt = (y F.) (i(() ne x:(1) di a:(1)
(49)
de:
dt
where 1i(() E IR is the externally applied force acting on the oscillator.
The characteristic equation for the unforced dynamics (i.e., (t( = 0) can be stated as 48
x² 2 Wn.s À wi. O, with the solutions λ1 Wh,t 1 Wh. C
(50)
where Wh,s = /R. and 4; पン are natural frequency and the damping factor oscillator, respectively. This harmonic oscillator exhibits three regimes: underdamped (C: 1), critically damped (G. 1), and overdamped regime (G, 1).
de
homogeneous
We approximate the forcing using the Heavyside function H(t): Et(t) F.(t.) H(t), where Fi(tx) is the constant external forcing as computed by (10). The solution for Çi / 1 is given by 48
(t((t)))) ys (1++1) (tu+1)) ((o...o c2,18.) cos(3, 5t) (c1.08, + (2.) sin 8t)) -- (S1)
(C1.s (o 5t) + c2 sin(8, 5t)) -a,
K
where St = tui1 tx. as Sa Wh.i. and 1i = Wn,t / C. After enforcing the initial conditions
(tx(,), (ta)), the integration constants
dx((1)) ai (Ii(le)
E Rg
C1,4 (1x))) Fi(te)
02.1 =
(52)
Na
AN
20

---

can be identified with AA; = As,2 A = -2 where ji is the imaginary value. While we could derive the solution for the critically damped case d, < 2/Kx separately, we instead approximate it in our network
12,3 j where R is small, positive
dynamics with by setting AA, = j value.
12 j 2 c
B.2 Algorithmic implementation
We can now leverage the closed form solution to evolution a single, decoupled damped harmonic oscillator (T) to solve the integral in [I]
x(1++1) ze
(v(te), u(te)),
(-1 i ((1211)
(er @ cos o1) + c2 o sin( o aa + F(te)
1
(53)
((e1
c2 (3 cos öt) 0,8 sin( Ge adi
with
d
diag( (11 Knn), d= (D11 Dan): Wa VR, 5
2/Kn" ( way wh v 11-5² 53
F((+) g(u((x)) (K- xx)(((((1) (D d) [1x) tanhi (Wx(I+) b),
F(tx)
[x(xx) F(fx)
x )) (54)
4 x(te)
C2 B ((tx) a
K
We summarize the approach of integrating/rolling out the CFA-CON dynamics in Algorithm|
Algorithm 1 Rollout of CFA-CON.
Inputs: initial state y(to), input sequence {u(tt).....uuf u((x). u(tw))
Outputs: state sequence (y(lo),..- -(())))), ((4x),..)((xx) v(N))
1: & -
2: while k S N do
3:
(I(tk) x(x)) + y(ta)
4:
F(tx) t g(u(ts) - K - K) x(ta) (D (tx)) - tanh (Wa(ta) +)
5:
Whr + VF: 20x
Compute the characteristics of the decoupled harmonic oscillators.
6:
a, 3 + -- Whi Wa c-
7:
a x((x)) F((e)
Compute integration constants using initial conditions.
8:
2 + 1/2 (i(tz) + +
(z(tx) F(L))
9:
8t = l++1 t
> Set time step.
10:
Update state with approximated closed-form solution. x(tx+1) (c) @ cos( (3 [t) + sin(3 öt)) 3 e +
11:
#(1++1) -- ((e) 02 + B) cBs(B 8t) (cr B0Be 120200 sin(3 (()) e
kt 1-k+1
Update time index.
14:
end while
B.3 Approximation bounds for CFA-CON
Lemma TT demonstrates how, for the particular of external input and linearly decoupled oscillators establish bounds the approximation
(which always free to choose), closed-form solution instead of the ground-truth coupled oscillator dynamics.
when using the
Lemma 11. Suppose that the network is unforced with g(u(t)) 0 and that K = diag(Ri,..... .Rn). D diag((d (dr..dd........ dn)) such that the oscillators are not linearly coupled. Then, given any 2 0, and the initial state (()) E R²" the ervor between the continuous dynamics (ft of 2 and the approximated dynamics &t(t) in
10 is bounded by [(((t) - a(t)]l S 2.
Proof.
10 and fa. (v()), 0) give us (fr,ia(v(t), + nld (v(1).0) -fana)(()) F);
[(t)
Kx(()) Dř Laml(( + 6) + Kw(() + Då + tanh(Wa(() + 6), ||- tan(( + b) + tanh(W=(()) b)ll < 2.
(55)
21

---

B.4 Empirical evaluation of approximation error
In Table 3 we present a comparison of CFA-CON with several other strategies for integrating nonlinear dynamics, such as CON. Following the implicit assumption made in Section [ we consider the case of g(u) = 0, K = diag(Ki,... ..... .K)) and D = diag(di...... dx) but with the hyperbolic coupling between the oscillators active (i.c., full W matrix). Integrating the dynamics at a very small time step (i.c., öt = 5e-5 s) with a high-order ODE solver would give us a very accurate solution, but this is computationally infeasible in practice. We, therefore, regard this as the upper bound on the accuracy the solution. A feasihle solution would he to implement either a high-order solver such as Tsit5 at a larger integration time-stop, c.g., = 1e-1 s) or low-order solver with slightly smaller integration time step, e.g.. 51 5e-2 2 s). Therefore, we also benchmark these options. We also benchmark an implementation specialized on the underdamped case (i.e.. & < 1): Closed-Fore Form Approximation specialized implementation allows us avoid using complex numbers in the algorithm and reduces the number of computations necessary for calculating the approximated solution. As result, we considerable increase in the sim-time to real-time factor.
the Underdamped Coupled Oscillator Network (CFA UDCON), This
B.4.1 Integration error
We perform the integration error benchmark over 100 different network configurations. all consisting of 50 oscillators (n = 50); First, we sample the natural frequency of the ith oscillator from uniform distribution as Wh.i (0.05 Hz, 0.5 Hz). then we sample x; ~ (((22 N/m, 2 2/mm) such that = diag(Ki..... Kn) >
which lets us determine each mass m, _ Next, the damping ratio is determined as G ~ (4 0.9) and
.,,
(0.11, 2.0) for the underdamped and general case, respectively. As result, D diag(d1, (d1, ..... d.)) > sample W and
with d, = 2 /InG K. given. Finally, by leveraging the Cholesky decomposition, ~ (d,-1, 1). We compute the estimation error of all integrated trajectories with respect to the high-precision solution (i.e., Tsitouras' 5/4 method (Tsit5) at ót 5e-4 s). For this, we compute the RMSE for each 60 trajectory and then take the mean and standard deviation across the 100 different network configurations.
B.4.2 Simulation-time to real-time factor
The simulation vs. real-time factor is computed as the simulated rollout duration per second of computational time, For this, we let cach method simulate 60s trajectory for 100 times and record the minimum run time on an Intel Core 17-10870H CPU (single core) over 10 trials. Because of computational constraints, we simulated with the high-precision Tsit5 solver the trajectory only 5 times.
B.4.3 Results
The results in Table show that CFA-CON is 30 % more accurate than the Euler integrator at half of the speed. Compared against the Tsit5 integrator, CFA-CON exhibits a 1.56x speed increase while heing significantly less accurate. For the underdamped case with Ç < 1, the specialized implementation CFA-UDCON is 14.8 % faster and at the same time 32 % more accurate than the Euler integrator. Furthermore, CFA-UDCON is 3.7x faster and significantly less accurate than the Tsit5 inlegrator. We can conclude that in the pure rollout setting (i.e., no backpropagation involved) for a generic CON, the CFA-CON does not show clear advantages to an appropriately tuned Fuler or Tsit5 solver. However, the specialized version CFA-UDCON demonstrates 2.4x speed-up at no reduction of accuracy vs. CFA-CON for underdamped oscillator networks.
22

---

Method
RMSE [m] 4
RMSE < [m]]
Complexity J
Sim. time Real time
CON with Tsit5 St 5e-5
n/a
lodg 7ph
(3.5el
5.68x
CON Tsit5 at 61 = 1e- a
Se-5 1e-5
8e-6 1e-3
(1.8888)
11310x
CON with Fuler 8t 5e-2
0.010 0.003
0.022 + 0.005
loga Th 61
(7.1ee7
36500x
CFA-CON (our) St = 1e-1s
0.007 0.002
0.015 0.003
loxa Vi
(3.5e7)
CFA-UDCON (our) with St = le- 1
n/a
0.015 1: 0.003
long Tn.
6t
(3.5ce7)
41900x
Table 3: Benchmarking of various methods for integrating the CON dynamics. The RMSE is computed with respect to the Tsitouras' S/4 method (Tsit5) (i.e. extremely high accuracy but also extremely high computational complexity). We denote with the number of oscillators in the network (in this case of the numerical ODE solver, and with öt the time-step. When stating the complexity. we refer to h tx to as the rollout horizon in seconds. In this casc, we report the results for a horizon of h = 60 s. The RMSE column states the RMSE of the various integration strategics with respect to the CON with Tsit5 at 8t 5e-55 5s solution, which we consider to be the ground-truth. The RMSE ( < 1 computes the same metrics, but this time for dataset that contains only underdamped oscillators. The the simulation achieved (in seconds) per second real- real-time (i.e., computational time). We report the mean and standard deviation RMSE over different network configurations.
50), with the order
column states the ratio between the duration of
23

---

C
Appendix on experimental setup and datasets
C.1
Datasets
For all datascts, we generate images of size 32 x 32px and subsequently normalize the pixels to the interval [-1,1].
C.1.1 Unactuated mechanical datasets
We consider multiple mechanical datascts based standard implementation included in the Toy Physics category of the NeurlPS 2021 Track on Datasets and Benchmarks publication by Botev et al. [26]: mass-spring with friction (M-SP+F), single pendulum with friction (S-P+F), and double pendulum with friction (D-P+F). All datasets contain 5000 system trajectories in the training set and 1000 trajectories each in the validation and lest set. Each trajectory is generated by first randomly initializing the system. then rolling it out for 3 using an Euler integrator with time step size of ms. Samples are recorded at rate of 20 Hz (i.c., a time step of 0.05 s). As a result, each trajectory contains 60 images of the system's 's state. As all of these datasets are unactuated, we can deactivate the impu-to-forcing mapping component from all models (e.g., set g(u) = 0 for the CON model).
The M-SP+F F dataset contains motion samples of damped harmonic oscillator with mass of 0.5 kg, a spring stiffness of 2 N/ m, and a damping coefficient of 0.05 Ns/ m. For each trajectory, the initial condition of the mass- spring is randomly sampled by combining a random sign(g) (q) with a uniformly sampled |4| ~ (((0.1 m, 1 m). The position of the mass is rendered with filled cirele in grayscale image.
The SP+F and DP+F datasets include the evolutions of a single-link link pendulum and double- link pendulum. respectively. with mass of 0.5 kg attached the end of each link, which has length of 1 m. The dataset considers a gravitational accelcration of 3 m/s". rotational damper with coefficient 0.05 Nms/ra provides the friction. Similarly to the M-SP+F dataset, both the sign and the absolute value of the initial configuration are randomly sampled, where (g0)] ~ (1.3 rad, 2.3 rad). The position of each mass is rendered with a filled for the double pendulum, each mass is circle. For the single-link pendulum, this is done in grayscale, rendered with a different color (i.e., blue and red).
C.1.2 Actuated continuum soft robot datasets
The shape of slender and deformable rods can be approximated by considering the deformations along the ID curve of the backbone [57]. While this curve is still infinite-dimensionall it is possible to discretize the backbone into (many) segments with piccewise constant strain (58) 57. Accordingly, we describe the kinematics of a planar continuum soft robot consisting of nis segments with the PCS model (58). We assume each segment has a length of 100 mm and a diameter of mm. The PCS model assumes each segment to have constant strain. In the planar case, this means that the shape of the ith segment can be parametrized by & [Mbe, Tsh.s axx,] E D where Khe, is the bending strain (i.e., the curvature) in the unit rad/un o/mm dah, is the shear strain (dimensionless). and aax,s is the axial elongation strain (dimensionless). The robot's configuration is then defined as q = [&]
<T, In the case of Piecewise Constant Curvature (PCC), only the hending strain is active as shear strains and axial strains are neglected, and the configuration is now q E R"b, The PCS model generates EOM in the form of (59]
B(q) q + C(y. a) 4 + G(g) + Kg9 + Dū u((),
(56)
where B(4)= > O and Cly. 4) are the inertia and Corioli matrices, respectively, G(4) collects the gravitational forces, Kg >- is the stiffness matrix, and D2 > 0 contains the damping coefficients. u(t) E Rd an extemal force acting on the generalized coordinates, and now m
7b..
We derive the corresponding dynamics for continuum soft robot of material density 600 600kg//m kg/mm, /m" elastic modulus of 20 000 Pa, shear modulus of 10 100 Pa, and damping coefficients of 0. 01 Nm for bending strains, 0.01 Ns for shear strains, and 0.01 Ns for axial strains, respectively. Gravity is pointing downwards. The implementation of the dynamics in JAX (60] is hased on the JSRM library [54 61]. and we simulate the
robot using a constant integration time step of O. mss We render grayscale images of the robot with a size of 32 x 32px at rate of 50 Hz using OpenCV (62] We generate 10000 trajectories. each of duration 2.0s and sampling time-step of 0.02 s. We use 60 % training, 20% % validation, and 20 % test split. For cach trajectory. we randomly sample a constant actuation/input u ~ (-(-mmax; (epax)) ). We choose the maximum actuation magnitude to be equal to the sum of the contribution of the potential forces (i.c., clastic and gravitational forces):
Wmax = G(mmax) + K fmax with fmax. = [5 rad/ m, 0.2, 0.2)".
We generate three datasets based this continuum soft rohot model: in the CS dataset, we consider one segment with all three planar strains active (i.e.. bending, shear, and elongation). This results in three DOF and six-state variables in the dynamical model. In the case of the PCC-NS-2 and PCC-NS-3 datasets, simulated system consisting of two planar bending segments, respectively. Each segment is parametrized using Constant Curvature (CC) [63] (4), which results in two configuration variables and a state dimension of four.
base the dataset on
24

---

C.1.3 Unactuated PDE
dataset
We consider the 1*-order Reaction-diffusing (R-D) PDE which (Champion al, 2019). 2 evaluated their SINDy Autocncoder on. The PDE of the high-dimensional lambda-omega reaction-diftusion system is defined
du (1 (u² v²))) 2 + 3 (u² + + di
(uu 8'u
0u
(57)
0u A(u" + + uu (u² + e²)) + d2
8v'
0'v
Øt
0q2
where u(t, g) R x R² + R and v(t, g) R R + R are time-dependent two vector fields defined over the
spatial domain y E R We choose the same system parameters and initial condition as Champion et al. [9]: di, d2 0.1, and and
(0 q) tanh (Voi q cos [(m iqg) / 42
(58)
(v + q sin <(q + +
V
(00 q) tanh
After discretizing the spatial domain into 32 points along each dimension, we solve the PDE with a MATLABB ODE45 solver solution u(t, q) and v(t, q) at each time step and grid point. Subsequently, the solution multiplied with Gaussian centered at origin [9]
a(t, 2) = exp(-0.01 0.01 (qf + x ū)t,
(59)
ū(t, exp(-0.01 (qi x v(t,
We integrate the system from the specified initial condition for 500 and store samples at time step of 0.05 We divide the entire sequence into 99 subsequences each containing 101 samples. We train the models to predict these subsequences that have horizon of 5.0s each.
We stack the solution of u(t, q) and v(t,g) contained in the grids ou()), (t) 32x33 respectively,
gather the images oft) E R 922232x2 containing two channels. A sample sequence of the generated images is presented in Apx. 17 We use 60 % of the subsequences (i.c., 59) as our training set, and employ 20 % (i.c., 19)
for the validation and test sets, respectively.
C.2 Autoencoder architecture
For the encoder and decoder, we rely on a vanilla CNNs implemented as a 3 VAE 5
Encoder. The encoder consists two convolutional layers with kernel size (3, 3) and stride (1, 1) mapping to 16, 32, respectively. The features flattened and then passed two linear layers with hidden dimension 256 and n.. Fach layer (except for the last) followed by layer norm (55) and LeakyRel.U nonlincarity.
The decoder first uses two linear layers to map to hidden dimensions of 256 and 32768, respectively. We then apply two 2D transposed convolutions [66] reducing the number of channels first to 16, and then to 1. Each layer (except for the last linear and last convolutional) is followed by layer norm [65 and a LeakyReLU nonlinearity, Finally, we apply sigmoid function to clip the output into the range [-1,
Decoder.
C.3 Latent dynamic models
In the following section, we provide implementation details for the latent dynamic models that we evaluated as part of this work.
C.3.1 Coupled Oscillator Network (CON)
We leverage the CON in W-coordinates given by 3 for learning latent space dynamics. Specifically, consider the input-to-force mapping g(u)) = (u) u(t), where B(u) e Rexm is parametrized hy few-layer MLP. We report results for two different sizes of the MLP: one medium-sized variant consisting five layers with a hidden dimension of 30 and a small variant with two layers and hidden dimension 12. In both cases, we use a hyperbolic tangent as a nonlinearity.
When training the model, we jointly optimize M Kw, Dw, b and g(u). However, we also need to make sure that adhere the stability constraints M Kw. Dw > O. For this, leverage the Cholvesky dccomposition [67]. Instead of directly learning the full matrix E X",, we designate the elements of an upper triangular matrix U X"= as the trainable parameters. The Cholesky decomposition demands that diag (U11.... Un > 0. Therefore, we apply the operation
Un = log (1 + Catee)) +
(60)
25

---

where Ü is the learned upper triangular matrix, and 61 1e 6 and 62 2e-6 are two small, positive values. The positive-definite matrix now given by U U
C.3.2 Neural ODEs
We consider two kinds of Neural ODEs (27]: the vanilla INODE: : E(t) x utt) + E(t) maps latent state and system actuation directly into time derivative of the latent state. In contrast, for the MECH-NODE, we enforce the latent dynamics to have mechanical structure
a(t)
{(t)
=
(ACCH-NNOD(E(t))) u(t))
(61)
We represent both SNODE and /MECH-NODS MLPs consisting of layers, a hidden dimension of 30, and hyperbolic tangent nonlinearity.
C.3.3 Autoregressive models
For the below stated autoregressive models, we divide the integration hetween two (latent) samples E(tx)) and {(tx+1) into Nint integration steps E(tx + 6t),...,((t,. + k'ōt)..... + Nmmrt) where &t is the integration step size and fa+1 = fg +N Nf. The autoregressive model now describes the transition E(fer41) far (E(t))))) u((+))) vk' L...., Nint.
RNN. We implement a standard, single-layer Elman RNN with tanh nonlinearity. The hidden state captures the latent state of the system. The latent state transition functions are given by
<((x+1)) = (aahh E(4x) + bbn + Wi u((x) + bm))
(62)
where Whn 2 bah E 2 Win E xm and by, E I22ms
GRU. We implement a standard, single-layer GRU (30) with si gmo d activation function where we interpret the latent state of the system as the hidden state the cell. The latent state transition functions are given by
r a (Wnr E(tx) bn. + Wi, u(ta) + bix)
a (Wip E((er) + bup + Wip r((()) + bip)
(63)
tanh (Win E(tar) bua)) Win u(tx) + bia)
<(te++1) (1 p) ④ n P E(lx)
where a is the sigmoid function, @ the Hadamard product, Whe, Whps Win R X2nx Wh. Wip, Win 2 xm and bue, bir, dips dm 2
CORNN. A lime-discrete coRNN is defined by the transition function
[x((e442) = 8(tw/ + St (-xz(t=)-=xtt =(((x) tanh (WE(tw) + Vu(ta) (((t))-b))) b)) [(te(+1)
z(4w) 8t z(tx)
{(te/11)=
(64)
where 7,5 Rt positive, scalar hyperparameters representing the stiffness and damping coefficients, respectively. The term tanh (WE(ft)) Vu(tz) + b) with W E R2n, x2n,, V E xm and E R" contributes nonlinear state-to-state connections. is implemented with a lincar layer operating on (E(te); ) (ttx)) followed by hyperbolic tangent nonlinearity.
CFA-CON. We adapt the Alg. foor predicting the time evolution in latent-space
<(tu+1) fCFA- (E(()))))())))))))))))))) m(te)),
(65)
where SCFA-CON describe the autoregressive state transition by the CFA-CON model as introduced in Eq 11
C.4 First-order variants of dynamical models
learning (latent) dynamics 1d (o-order systems (c.g.. the reaction reaction-diffusion dataset R- D), might beneficial also to formulate the dynamical model to of s (ooorerr. While this is straightforward for some dynamics that do not explicitly take the order account (c.g., RNN, GRU, NODE), for other models such as coRNN, CON, and CFA-CON more adjustments are necessary. Namely, we substitute the da component of the ODE with the as expression for dt- Furthermore, we remove any terms that depend on the velocity å (e.g., damping effects). Below, we report in detail the adapted, -rrrreer formulations for the coRNN, CON, and CFA-CON models.
de
CON. the or-order version, we adapt the standard, 2nd -oddr ODE of the CON network as defined in to
{t)) a(t) = M (o(u(t)) - Kwa(f) - tanh(z(t) + b))
(66)
26

---

CORNN.
In the 1ororder version, we define the transition function
{(tv+4)) =(tw41) z(twe 6t 2(tw) ) öt tanh (WE(tw)) Vu(ta) + 6).
(67)
CFA-CON. We adapt 1"-order version of Alg. 1 for predicting the time evolution in latent-space
<(4x+1) =(40++1) = 2((x))+ (
fortte
- ((w))x=((()))
(68)
F(te)= g(u(tt)) (K- xx)=(tt)) tanh(W b),
where the closed-form solution for the integral is given by
( F(tx) - xx((t) dt') F((N) (1 (1- e
((w))naz(tt))
(69)
C.5 Estimation of the initial latent velocity
For 2nd -order systems and when integrating the evolution of the latent state {(I) = [=- 5 =T())]" in time, we need have access initial latent velocity ((4) such that we roll the latent state {(t) in time. naive approach estimating such initial latent velocity would be encode multiple (at least two) images of the system the start of the trajectory into latent space and then perform numerical differentiation (c.g. finite difterences) in latent space. However, found the resulting ((t, to be relatively noisy and susceptible small encoding errors. Instead, velocity into latent space using the encoder's Jacobian. First, estimate the image-space velocity at ta using o(f+=1)-o(tt-1)) The latent velocity is then estimated ((tx) Do
propose perform numerical differentiation in image space and then map this
finite differences: ō(tx) 2 f
as (o(ti)) ô(tx).
where 02 is obtained with forward-mode automatic differentiation.
do
C.6 Training
We implement the network dynamics and the neural networks (e.g.. encoder, decoder, and MI.Ps) in JAX (0)
and Flax [68]. respectively. When training or inferring time-continoous dynamical models (e.g.. NeuralODE.
CON), we rely on Diffrax (47) for numerical integration
the ODE using the Dormand Prince's s 5/4 method (69
(Dopri5). For the numerical integration of both the coneicontiuuuus and the time-discrete models (c.g.. RNN, coRNN, CFA-CON), we use an integration time-step St of 0.025 s and 0.01 robotic datasets, respectively.
the Toy Physics (26] and soft
Because of the GPU memory constraints, we limit ourselves batch size of and 80 trajectories for the Toy Physics 266 and soft robotic datasets, respectively. We implement learning rate schedule consisting of warm-up (5 epochs) and cosine annealing (70] period (remaining epochs). We employ Adam optimizer 71 [22 with 81 0.9, 82 0.999 for updating both neural network weights (e.g.. encoder, decoder) and parameters of the dynamical model (c.g., K, Dw, Mw. etc.).
Before training, we conduct a hyperparameter selection study using Optuna (53) For this, we leverage a Tree-Structured Parzen Estimator (73) for identifying hyperparameters such as the base learning rate, the weight decay, the loss function weights, and model-specific hy hyperparameters such as the number of MI.P layers, the hidden dimension of the MLP layers, the y and values for the coRNN model etc. that minimize the RMSE of the predicted images. To reduce computational requirements, we employ the Asynchronous Successive Halving Algorithm [74] to stop unpromising trials carly.
C.7 Evaluation metrics
we state the RMSE, the Peak Signal-to-Noiss Ratio (PSNR) and the Structural Similarity Index Measure (SSIM) 78] between the ground-truth image o E R and the predicted image image ó E RR guu,
Similar to other publications in the field 75 76 77
We use the separated test set for all evaluation results.
C.7.1 Root Mean-Square Error
The RMSE between the two images is given by
(Owv oun) ho We
RMSE(0, 6)
(70)
27

---

C.7.2 Peak Signal-to-Noise Ratio
The PSNR function of the total Mean Squared Error (MSF) loss and the maximum dynamic range of the image L.
Po 6) = 20 (L) 10 log
Σ
(Ouo Ous))
(71)
ho Wo
As we work with normalized images with pixels in the interval [-1, 11. the dynamic range is L = 2.
C.7.3 Structural Similarity Index Measure
As simple pixel-by bixeeell metrics such RMSE PSNR tend to average out any encountered errors, this could lead to a situation in which significant reconstruction error in part of the image is not seen in the RMSE metric but has a huge impact the visual appearance of the reconstruction. SSIM [78) incorporates not just the absolute errors, but also the strong inter-dependencies between pixels, especially when they are spatially close. The SSIM metric between two observations and ô is given by
SSIM(o, 6) = 1 (o
(72)
where
1(o, 6) = 2/ooss + +1
+ + C2' c(o, 0) 20.00 +
s(o, Tob + C
(73)
d3 +
0.00 + C3
We use the constants C₁ L)) 2 C2 (kaL)² and C3 C2/2, where L signifies the dynamic range
previously used for the PSNR metric, and k = 0.01 and l22 = 0.03. The average and the variance x² is computed with a Caussian filter with a ID kernel of size 11 and sigma 1.5. We set the weight exponents a, B. and y for the luminance, contrast, and structure comparisons all onc. We rely on the PIX library (79 for efficiently computing the SSIM metric.
C.8 Compute resources
We trained the models on several desktop workstations for a total duration of roughly 150 h. In total, we relied on 10x RTX 3090/4090 GPUs, each with 24 GB of VRAM, training the models in parallel. Each workstation contained between 64 and 128 GB of RAM. and we used roughly 100 GB of total storage. Training each model on one random seed took between 45 min and 4 h depending on the model type, the integration time constant, and the number of trainable parameters. The hyperparameter tuning we conducted beforehand (only on one random seed) took roughly the same time and computational resources as generating the final results.
For the control experiments, we additionally used laptop with 16-core Intel Core 17-10870H CPU and 32 GB RAM. We did not need to use GPU for evaluating the model during closed-loop control.
28

---

D Appendix on learning latent dynamics
We report the full set quantitative results, including the additional evaluation metries PSNR and SSIM, in Tables to For PCC-NS-2 Tab. [8] we additionally also recorded the training steps per second on Nvidia RTX 3090 GPU with batch size of 80 (which leads to 8080 images per batch). We plot the results
sweep
across the latent dimensions for the additional evaluation metrics PSNR and SSIM Fig. Correspondingly.
visualize the number of trainable parameters of each model vs. the latent dimension in Fig. 9 In Figs. 10 16
we present sequences of stills for the rollout of the trained CON(-M) models the various datasets.
Model
RMSE
PSNR f
SSIM
Parameters
RNN
0.2739 + 0.0057
1.16 + 0.02
0.6958 + 0.0122
88
GRU 30
0.0267 + 0.0033
6.13 + 0.09
0.9861 + 0.0022
248
CORNN 35
0.0265 + 0.0002
6.13 0.01
0.9853 + 0.0006
40
NODE 27
0.0264 + 0.0010
6.14 0.03
0.9858 + 0.0009
3368
MECH-NODE
0.0328 + 0.0034
5.99 + 0.07
0.9821 + 0.0024
3211
CON (our)
0.0303 + 0.0053
6.05 + 0.13
0.9817 + 0.0027
34
CFA-CON (our)
0.0313 + 0.0026
6.02 + 0.06
0.9843 + 0.0008
34
Table 4: Benchmarking of CON and CFA-CON at learning latent dynamics on the M-SP+F (mass-spring with friction) dataset. For all models, latent dimension of 4 is chosen. As this dataset does not consider any inputs, we remove all parameters in the RNN, GRU, coRNN, CON, and CFA-CON models related to the input mapping. MECH-NODE is NODE with prior knowledge about the mechanical structure of the system (i.e.. di .). We report the mean and standard deviation over three different random seeds the number of parameters of each latent dynamics model.
de
Model
RMSE
PSNR f
SSIM
Parameters
RNN
0.2378 + 0.0352
4.31
0.7568 + 0.0350
88
GRU 30
0.1457 1 0.0078
4.78 0.05
0.9168 1 0.0093
248
cORNN 35
0.1333 + 0.0044
4.86 0.03
0.9191 11994:00 0.0055
40
NODE 27
0.1260 = 0.0013
4.91
0.9379 10879-00 0.0009
3368
MECH-NODE
0.1650 + 0.0205
4.67 0.12
0.8985 + 0.0153
3244
CON (our)
0.1303 + 0.0064
4.88 0.04
0.9175 + 0.0095
34
CFA-CON (our)
0.1352 1 0.0073
4.85 1 0.05
0.9133 1 0.0052
34
Table 5: Benchmarking of CON and CFA-CON learning latent dynamics the S-P+F (single pendulum with friction) dataset. For all models, latent dimension of n. 4 is chosen. As this dataset does not consider any inputs, remove all parameters in the RNN, GRU, coRNN, CON, and CFA-CON models related to the input mapping. MECH-NODE is NODE with prior knowledge about the mechanical structure of the system (i.c., dI 2). We report the mean and standard deviation over three different random seeds and the number of parameters of each latent dynamics model.
de
Model
RMSE 4
PSNR 1
SSIM 1
Parameters
RNN
0.1694 0.0004
4.631 = 0.002
0.7082 + 0.0032
672
GRU 30
0.1329
0.0005
4.858 + 0.003
0.8340 + 0.0021
1968
caRNN 35
0.1324 0.0016
4.862 0.012
0.8229 + 0.0039
348
NODE 27
0.1324 0.0024
4.861 1 0.016
0.8101 1 0.0024
4404
MECH-NODE
0.1710 0.0111
4.624 0.063
0.7170 +: 0.0439
4032
CON (our)
0.1323 0.0018
4.862 + 0.013
0.8067 + 0.0038
246
CFA-CON (our)
0.1307 0.0012
4.873 1 0.008
0.8147 + 0.0034
246
Table 6: Benchmarking of CON and CFA-CON learning latent dynamics on the D-P+F (double pendulum with friction) dataset. For all models, latent dimension of nı 12 is chosen. As this dataset do not consider any inputs, we remove all parameters in the RNN, GRU, coRNN, CON, and CFA-CON models related to the input mapping. (i.c., de parameters of each latent dynamics model.
ELECH-NODA is NODE with prior knowledge about the mechanical structure of the system dt æ). We report the mean and standard deviation over three different random seeds and the number of
29

---

Model
RMSE
PSNR f
SSIM
# Parameters
RNN
0.1011 1 0.0009
25.92 1 0.08
0.9777
0.0004
696
GRU 30
0.0100
24.99
0,9730 = 0.0040
2040
CORNN
35
0.2537 0.0018
0.006
0.8820 0.0024
336
NODE 27
0.2415 0.0021
18.36
0.8946 0.0023
4374
MECH-NODE
0.2494 1 0.0028
18.08 0.10
0.8898 0.0016
4002
CON-S (our)
0.1993 +: 0.0646
20.03
2.44
0.9218 0.0380
1386
CON-M (our)
0.1063 + 0.0027
25.49 0.22
0.9758 0.0011
8568
CFA-CON (our)
0.1462 1 0.0211
22.72 1.17
0.9573 0.0103
8568
Table 7: Benchmarking of CON and CFA-CON at learning latent dynamics on the CS (soft robot with one
constant strain segment) dataset. For all models, a latent dimension of Ttz = 12 is chosen. CON-S and CON-M are small and medium-sized versions of the CON model, respectively. MECH-NODE is a NODE with prior knowledge about the mechanical structure of the system (i.e., da æ). We report the mean and standard deviation over three different random seeds and the number of parameters of each latent dynamics model.
Model
RMSE 4
PSNR
SSIM
# Paramelers
Init slepe scsnd
1 Inf. time [ins] 1
RNN
0.1373 0.0185
23.27
+ 1.10
0.9613 + 0. 0.0077
320
1.87
02.6
GRU 3D
0.0951
0.0021
26.45
0.19
0.9730
0.0040
928
1.83
03.2
coRNN
0.2504
0.0899
18.05
2.66
0.9814
152
02.7
NODE 27
0.1867 1. 0.0561
20.60
1: 2.28
0.8774
1 0.0857
3856
0.79
50.2
MECH-NODE
0.1035
0.0012
25.07
0.06
0.9778
0.0004
3062
50.3
CON S (our)
0.0996
0.0012
26.05
0.11
0.9792
0.0007
676
0.78
50.2
CON-M (our)
0.1008
0.0006
25.95
0.05
0.9786
O. 0.0003
7048
0.60
60.1
(our)
0.1124
25.01
0.9734
Table Benchmarking CON and CFA-CON at learning latent dynamics on PCC-NS-2 (soft robot with constant curvature segments) dataset. For all models, a latent dimension of n. 8 chosen. CON-S and CON-M are small and coddium-sizes versions the CON model, respectively. MECH-NODE,
a NODE dar 1). We report the mcan and with prior knowledge about the mechanical structure of the system (i.c., az standard deviation over three different random seeds. Furthermore, we state the number of parameters of each latent dynamics model and the training steps per second on a Nvidia RTX 3090 GPU. Fach batch contains 80 trajectories and 8080 images of resolution 32x32px in total. Finally, we report the inference time averaged over 5000 runs for performing a rollout of 2.02 s (while encoding and decoding all images along the trajectory) on an Nvidia RTX 3090 GPU with a batch size of 1.
D.1 Results for Reaction-diffusion dataset
As all previous examples exampled ODEs, strive the proposed approach also system that governed by PDEs. Specifically, consider the To address the unactuated nature of the dataset, remove, analog the M-SP+F, M- S-P+F, and D-P+F datasets, the input-to-state mapping parameters of the dynamical models (c.g., the B(u) and E(T) MLPs for the CON models). Furthermore, the PDE describing system dynamics is ordrr.... Therefore, we leverage 1"-order versions the latent dynamics as specified in Apx. C.4
(R-D) dataset introduced in Apx. C.
1.1.
We report the metries of the test set evaluations in Tab. 10 Furthermore, also present sequence of stills of the rollout of trained latent dynamics CON model in Fig. 17 We find it impressive that CON with its strong stability guarantees can accurately model the dynamics of high-dimensiooal PDE system.
30

---

Model
RMSE
PSNR f
SSIM
Parameters
RNN
0.2232 + 0.0075
19.05 0.29
0.8955
0.0083
696
GRU
30
0.2148 + 0.0196
19.38 0.76
0.9039
0.0223
2040
CORNN
35
0.2474 0.0018
18.15 0.06
0.8877
0.0011
336
NODE
27
0.3373 + 0.0565
15.46 1.34
0.7432 0.0935
4374
MECH-NODE
0.1900 + 0.0024
20.45 +.11
0.9315
0.0011
4002
CON-S (our)
0.1792 1 0.0038
20.96
1 0.18
0.9392
0.0023
1386
CON-M (our)
0.1785 1 0.0023
20.99
0.11
0.9395
0.0018
8568
CFA-CON (our)
0.1803 + 0.0003
20.90
0.01
0.9366 = 0.0004
8568
Table Benchmarking of CON and CFA-CON learning latent dynamics the PCC-NS-3 (soft robot with three constant curvature segments) dataset. For all models, latent dimension of 12 is chosen. CON-S NODE as = á). We report the mean and
and CON-M small and medium-sized versions of the CON model, respectively. MECH-NODE with prior knowledge about the mechanical structure of the system (i.e., da standard deviation over three different random seeds and the number parameters each latent dynamics model.
Model
RMSE
PSNR +
SSIM 1
Parameters
RNN
0.3763 = 0.0374
3.82
0.4463 = 0.1358
20
GRU 30
0.3232 + 0.0368
3.9910.13
0.6798 0.0949
52
order coRNN
0.0741 + 0.0001
5.05
0.9724
0.0014
20
NODE
0.0738 = 0.0007
5.06
0.9683
0.0022
3064
CON (our)
0.1110 0.0160
5.03
0.9372 0.0109
24
CFA-CON (our)
0.1068 + 0.0059
5.05 + 0.05
0.9418 0.0026
24
Table 10: Benchmarking of CON and CFA-CON at learning latent dynamics on the R-D -defftio-diffluiio) dataset. For all models, a latent dimension of n, = 4 chosen. As this dataset does not consider inputs, we remove all parameters in the RNN, GRU, coRNN, CON, and CFA-CON models related to the input mapping. Also, as reaction-diftusion system is governed by 1o-order PDE dynamics, we use specialized, 1 (o-order version of CON, CFA-CON, coRNN dynamics. report the mean standard deviation over three different random seeds and number parameters of each latent dynamics model.
31

---

PSNR
PSNR
16
HEUN MERCH
IONS
NODES
10
10
GPA-DON
CON
10
2000
4000
600D
8000
10000
12000
75g
Model porameters
PSNR latent dimension fh
(b) PSNR model parameters
1.0
AP
SSDM
0.8
SSIM
0.8
MISCH NOTR
MIRCH NODE
S
CON-M CPA CON
COX-M
CFA CON
12000
Model parameters
(c) SSIM latent dimension
(d) SSIM model parameters
Figure Evaluation prediction performance various models vs. the dimension of their latent representation n, and the number of trainable parameters of the dynamics model, respectively, We optimize the hyperparameters for case nz = 8, and execute the tuning separately each model and dataset.
12500
RNN
MECH-NODESS
Model paranelers
10000
MECH-NODE
CON-S
CON-M
7500
CFA-CON
5000
2500
10
15
20
25
30
78
Figure 9: Plot of number of trainable parameters vs. the latent dimension of various models trained on the PCC-NS-2 dataset. As we have configured them, CON-M and CFA-CON always have the same number of parameters (i.c., overlaying lines).
32

---

(a) s00.s
(b) s5.55
(c) 11.10s
(d) t=1.65
(c) 1=2.20s
(f) 1
Figure Prediction sequence of CON model with latent dimension nt, trained damped harmonic oscillator (M-SP+F) dataset [26]. Top row: Ground-truth evolution the system. Bottom row: Predictions of the CON model.
prediction model is given three images centered around for encoding the initial latent =(0) and estimation of the initial latent velocity (0). Subsequently, roll out the autonomous network dynamics (i.c., unforced) and compare the decoded predictions with the ground-trun truth evolution the system.
(a) s0.0s
(b) s50.5
(c) 11.10
(d) 1=1.65s
(e) 1=2.20
(f) 2
Figure 11: Prediction sequence a CON model with latent dimension n. 4 trained on single pendulum with friction (S-P+F) dataset [26). Top row: Ground-truth evolution system. Bottom row: Predictions of the CON model.
The prediction model is given three images centered around estimation of the initial latent velucity 2(0). Subsequently, roll out the autonomous network dynamics (i.e., unforced) and compare the decoded predictions with the ground-truth evolution of the system.
encoding the initial latent =(0) and
33

---

(a) s0.0s
(b) s5055
(c) 11.10s
(d) t=1.65
(c) 1=2.20s
1
Figure Prediction sequence of model with latent dimension nt, trained on the double pendulum with friction (D-P+F) dataset (26]. Top row: ground-truth evolution of system. Bottom row: predictions of the CON model.
prediction model is given three images centered around for encoding the initial latent =(0) and estimation the initial latent velocity (0). Subsequently, roll out the autonomous network dynamics (i.c., unforced) and compare the decoded predictions with the ground-trutt truth evolution the system.
(a) s0.00
(b) 1=0.44s
(c) 0
(d) 1=0.12s
(e) 1=0.16s
(f) t-0.20s
Figure 13: Prediction sequence of forced CON model with latent dimension n. trained on soft robotic dataset containing trajectories of simulated constant strain robot with one segment. Top row: Ground-truth evolution of the system. Bottom row: Predictions the CON-M model. prediction model is given three images centered around
encoding the initial latent =(0) and estimation of the initial latent velucity 2(0). Subsequently, roll out the autonomous network dynamics (i.e., unforced) and compare the decoded predictions with the ground-truth evolution of the system.
34

---

(a) t00.ss
(b) s
(c) s
(d) t=0.9
(c) 1=1.2s
Figure 14: Prediction sequence of an unforced CON model with latent dimension nt 8 trained on the PCC-NS-2 dataset. Top row: Ground-truth evolution the system. Bottom row: Predictions model.
CON-M
prediction model is given three images centered around estimation the initial latent velocity 2()). Subsequently, we roll out the autonomous network dynamics (i.e.. unforced) compare the decoded predictions with the ground- truth evolution the system.
encoding the initial latent =0) and
(a) t00.0s
(b) t003s
(c) 5
(d) t=0.95
(c) t=1.2s
Figure 15: Prodiction sequence of forced CON model with latent dimension n. 8 trained on the PCC-NS-2 dataset. Top row: Ground-truth evolution of system. Bottom row: Predictions the CON-M model. prediction model is given three images centered around estimation of initial latent velocity 2(0). Subsequently, we provide same constant imput u to both simulator and the network dynamics (i.c., unforced) and compare the decoded predictions with the ground-truth evolution of the system.
encoding the initial latent =(0) and
35

---

(a) s=0.8
(b) s
(c) 100.2s
(d) 1=1.08s
(e) 1=1.44s
(I) t=1.80
12 trained on the soft robotic PCC-NS-3 dataset containing trajectories of a simulated piccewise constant curvature robot with three segments. Top row: Ground-truth evolution of the system. Bottom row: Predictions of the CON-M model. The prediction model is given three images centered around t = 0 for encoding the initial latent z00) and estimation of the initial latent velocity 3(0). Subsequently, we roll out the autonomous network dynamics (i.e. unforced) and compare the decoded predictions with the ground-truth evolution the system.
Figure 16: Prediction sequence of forced CON model with latent dimension
00.0s
(b) 1=1.0s
(c) t=2.0s
1=3.0
(c) t=4.0s
(f) t=5.0s
Figure 17: Prediction sequence of an unforced, 1 CON model with latent dimension 4 trained on the reaction-diffe cofttiondiffusion (R-D) dataset. Top row: Ground-truth evolution of the system. Bottom row: Predictions of the CON-M model. We roll out the autonomous, order network dynamics (i.c., unforced) and compare the decoded predictions with the ground -rutth evolution of the system.
36

---

E Appendix on latent-space control
E.1 Latent-space control of damped harmonic oscillator
We consider an actuated version of the M-SP+F dataset (i.e., damped harmonic oscillator) and denote it with M-SP+F+A. All system, trajectory sampling and rendering parameters remain the same, except that for each trajectory in the dataset we randomly sample a forcing v ~ ((11N, N).
We train CON model with latent dimension ni = 1 over three random seeds on the M-SP+F+A dataset. This means that the network consists of a single oscillator. From the three different random seeds, we choose the model that achieves the best validation loss, which results in an RMSE of 0.0327, a PSNR of 5.99, and SSIM of 0.9796 on the test set.
shows how the encoder learns an almost linear relationship between the actual configuration of the system and the predicted latent space representation. Furthermore, we notice that both the ground-truth and the learned potential energy are convex and exhibit global minimum at 0
Fig. [8]
We compare the performance of P-satI-D, D+FF, and P-satl-D+FF controllers based the CON model in Fig. 19 For the P-satl-D controller, we choose the control gains K, 10, K: = 10, Ka 5, v = 1. The D+FF controller uses Ka 3.5. Finally, the P-satt-DDFF FF is configured with K 2, K; 0.3, Ka = =3.5, e 1. The results show that the sutl-D+FFF controller exhibits thanks to its feedforward term no overshooting and a faster response time than the pure feedback controller P-sail-D. The high accuracy of the feedfoward term can be seen from the performance of the D+ FF controller, that only exhibits relatively small iteady-state error. Adding small proportional and integral feedback actions in the P RFFI-D+FF controller keeps the compliance high while removing the steady-state error and reducing the response time.
Finally, we visualize the behavior of the P-satl-D+ controller as a sequence of stills in Fig. 20
1.0
016
0-1
0.7
-32
-1.00
-0.75
35
0955
0.50
175
1.00
75
-0.00
2.00
0.50
0.75
[e]]
[mc]
(a) Configuration g to latent & mapping
(b) Ground-truth and learned potential energy 24
Figure 18: Panel (a): Learned mapping from configuration to latent space for the CON model with n. (i.c., consisting of single oscillator) trained on the actuated damped harmonic oscillator (M-SP+F+A) dataset. Panel (b): The blue line represents the ground-truth potential energy of the damped harmonic oscillator. The orange line represents the learned potential energy of the CON model evaluated vs. the system configuration by rendering and subsequently encoding into latent space each configuration value.
E.2
Latent-space control of two segment PCC soft robot
E.2.1
Potential energy landscape
When leveraging (learned) dynamical models for setpoint regulation, it is essential to accurately estimate the potential energy as this dictates the efficacy of the feedforward terms. Therefore, we qualitatively evaluate the potential energy landscape of the CON latent dynamic model.
In Fig. [1 we can see how CON contains a single, isolated, and globally asymptotically stable equilibrium as proven In Appendix A.2 and Section 2
respectively.
Furthermore, we want to verify that the learned potential corresponds to the actual potential energy of the simu- lated system. An autonomous continuum soft robot with the tip pointing downwards in straight configuration exhibits isolated, globally asymptotically stable equilibrium q 0 (i.c., zero strains) [42]. For this purpose, we compare the learned potential energy field in Fig. 21(b) with the ground-truth potential energy field in
Fig. (1(c) We confirm, based on Fig, [1b that, indeed, the learned potential also has its minimum close to/at
q U. Although the field is shaped slightly differently, the potential forces clearly pointing inwards towards the global attractor.
37

---

0.2
0.50
Caniggaakk (bi)
0.25
Targe P
wrishb
0.00
PoaiL-D
0.0
-0.25
Paut-DDFFF
Poat-DeFF
-0.00
-0.1
-0.2
10
20
30
30
30
Time
Time
(a) Configuration g(t) E R
(b) Latent representation z(1) E IR?
Poun--D
D&FF
P-eat-D it D4FP
Contral taput [N)
2
Po
0.00
P-aat-DDFF
-2-
0.02
0.00
15
Timn
Time
(c) Control input u(t) IR?
(d) Potential energy (t) E R
Figure 19: Latent-space control of an actuated damped harmonic oscillator (M-SP+F+A) following sequence of setpoints. We compare multiple controllers based on trained CON network with nı = 1. The CON model weights are initialized using a random seed of O. The blue line represents a pure feedback controller (P-salI-D) The orange line visualizes the behavior (D+FF), The green line shows the performance our proposed combination of feedback and feedforward terms (P-satl-D++ FF). The dotted and solid lines show the reference and actual values, respectively. For each setpoint, randomly sample desired shape q² and render the corresponding image o This image is then encoded to target latent 24 The controller then computes a latent-space torque pd which is decoded to an input u. Finally, we provide this input to simulator, which performs a roll-out of the closed-loop dynamics. Important: The robot's 'sss configure configuration (i.c., first-prit princciie,,, minimal-order state) is solely used for generating target image simulating the closed-loop system.
feedforward controller with only damping term applied in feedback
E.2.2 Model selection
the control experiments. we train instances of MECH-NODE and CON-M models with latent dimension Y& and with neural network weights initialized with three different random seeds. MECH-NODE, we choose the model with the lowest validation loss (seed
For the CON network, found that model-based control does not perform as well when the latent stiffness Tw (as visualized in Fig. (2)))) is significantly larger along one of the Eigenvectors than along the other one. Therefore, we evaluate the Eigenvalues of the learned stiffness matrix in W-coordinates after training: (T Particularly, we choose the seed that minimizes the normalized standard deviation of the Eigenvalues
(T >T()w))
(A (Tw) nA) (Ne(Tw)) nA)
(74)
2
seed arg min
E.2.3 Additional control results
Additional results for the P-satI-D feedback controller based on the MECH-NODE and CON models are provided in Fig. (3] Fig. [4 respectively and for the P-satI-D+FF controller hased on the CON model in Fig. Sequences of stills for the CON P-satI-D+FF FF controller are provided Fig. 22
38

---

(a) s=0.0s
(b) 1=0.5s
(c) t=1.0
(d) t=1.5s
(e) 1=2.0s
(l) Target 0 4.9s
(g) 155.0s
(h) 5
(i) t=
(i) t=6.5s
(k) 1=7.5s
(l) Target 5-9.99
Figure 20: Sequence of closed-looe loop control of actuated damped harmonic oscillator (M-SP+F+A) with
P-sat- Satt-D+FFF D*Foonrrrolerr based trained CON with n, 1. Columns
actual behavior of
closed-loop system. Column 5: demonstrates the target image that the control for all time instances in the row.
39

---

0.521
0.3
0.25
0.252
0.1
0.216
0.1
-0.08
02
0.072
0.036
-04
0,000
-0.15
-0.10
00.05
0.00
0.05
0.10
0.15
(a) Learned potential energy a function of
0.397
0.1016
0284
CUONTT
0.98
0.10
[alin]
D/0728
0.000
0U0864
0.06
0/0012
15
0000
(ead
(b) Learned potential energy as function g
(c) Ground-truth potential energy as function of
Figure 21: Potential energy landscapes of a CON with = 2 trained to learn the latent space dynamics of continuum soft robots (simulated with two PCC segments). Panel (a): Here, we visualize the learned potential energy of CON using the color scale function of the latent representation the gradient of the potential field d84
Zw R2.. The arrows denote 22 (i.e., the potential force), with the magnitude of the gradient expressed the length of the arrow. Panel (b): Again, display the learned potential energy of CON using the color scale, but in this case, as function of the configuration y R of the robot (that is hidden from the model), First,
we render an image o the shape the robot for each configuration & = [9] (2]
IR. Then, we encode
the image into latent space as z ((o). This allows us then to compute the potential energy U((z) of the CON latent dynamics model. Panel (c): Here, we display the potential energy and its associated potential forees the actual (i.c., simulated) system.
40

---

(a) s=0.ss
(b) 10.2s
(c) 10.4s
(d) 1=0.6s
(e) t=0.8s
(I) Target 4.9s
(g) 15.00s
(h) 1
(i) 5
5
(k) 15.88
Target 9.99
Figure 22: Sequence closed-loop control of continuum soft robot consisting of two constant curvature segments with the P-satl-D+FF based on trained CON with n. Columns 1-4: show the actual behavior closed-loop system. Column demonstrates the target image that control sees for all instances the row.
10
0.01
Culliqunibe jad/ui
0-
0.00
-5-
-0.01
-10
00,05
-15
[H]
(a) Configuration g(t) E R
(b) Latent representation z(t) IR
0,004
Contrdd Eapet [Nua]
0.002
0,000
0.002
-0.004
30 15
25
30
35
Timae (4)
(c) Control input u(t) IR2
Figure 23: Latent-space control of continuum soft robot (simulated using PCC segments) following sequence of setpoints with pure P-satI-D feedback controller operating in 2D latent space learned with the MÉCH-NODE model. The CON model weights are initialized using random seed of O. The dotted and solid lines show the reference and actual values. respectively. For each setpoint, randomly sample desired shape g and render the corresponding image o This image is then encoded to a target latent zd The controller then computes satent-space torque fa which decoded to input Finally, we provide this input to the simulator, which performs roll-out of the closed-loop dynamics. Important: The robot's configuration (i.e., the first- princiiiee,,,, minimal-ort state) is solely used for generating target image and simulating the closed-loop system.
41

---

POOORS
0.2
Cintigzzaion (cos))n))
Latent sorssile
-5-
5 St
-0.1
-10
-0.3
15
10
15
20
20
30
25
30
33
(a) Configuration q(t) E R
(b) Latent representation x(1) E uR?
0.20
0.20
0.15
0.15
Dime
0.30
0.30
0.05
0.05
Time [s]
Time
[H]
(c) Control input u(t) UR?
(d) Potential energy U(L) E
Figure 24: Latent-space control of continuum soft robot (simulated using two PCC segments) following
sequence of setpoints with pure P-satI-D feedback controller operating in 2D latent space learned with the CON model. The CON model weights are initialized using random seed of O. The dotted and solid lines show the reference and actual values, respectively. For each setpoint, we randomly sample desired shape q" and render corresponding image oº. This image is then encoded a target latent 2d. The controller then computes a latent-space torque r², which is decoded to input Finally, we provide this input to the simulator, which performs roll-out of the closed-loop dynamics. Important: The robot's configuration (i.e., the first principle, minimal-orde order state) is solely used for generating target image and simulating the closed-loop system.
42

---

20
Cinligzzaion (coo))o)
10
-10
5
-022
-20
-0.4
-0.0
10
15
20
30
30
135
(a) Configuration q(t) E R
(b) Latent representation z(1) E uR?
0.4
Dieme
0.2
0.2
Perin
0.1
0.1
Time [Hl
Time [el
(c) Control input u(t) uR?
(d) Potential energy U(L) E E
Figure 25: Latent-space control of continuum soft robot (simulated using two PCC segments) following sequence of setpoints with pure P-satI-++FF feedback & feedforward controller operating in 2D latent space learned with the CON model. The CON model weights are initialized using random seed of 0. The dotted and solid lines show the reference and actual values, respectively. For each setpoint, we randomly sample desired shape a and render corresponding image oa. This image then computes a latent-space torque 1rl, which is decoded to input Finally, we provide this input to the simulator, which performs roll-out of the closed-loop dynamics. Important: The robot's configuration (i.e., the first principle, canimal-orded state) is solely used for generating target image and simulating the closed-loop system.
encoded to target latent 2 The controller
43

---

F Extended discussion on future applications and limitations
F1 Systems for which we would expect the proposed method to work
Mechanical systems with continuous dynamics, dissipation, and a single, attractive equilibrium point. The proposed method is a very good fit for mechanical systems with continuous dynamics, dissipation, and a single, attractive equilibrium point. In this case, the real system and the latent dynamics share the energetic structure and stability guarantecs. Examples of such systems dominant elastic behavior, and other mechanical structures with elasticity.
many soft robots, deformable objects with
Local modeling of (mechanical) systems that do not meet the global assumptions. Even if the global assumptions of the proposed method are not met, the method can still be applied to model the local behavior around local asymptotic equilibrium point of the system (i.e., in the case of multi-stability) For example. the method could be used to model the behavior of a robotic leg locally in contact with the ground, a cobot's interaction with its environment, etc.
F.2 Systems for which we could envision the proposed method to work under (minor) modifications
Mechanical systems without dissipation. The proposed method would currently not work well for mechanical systems without any dissipation, as (a) the original system will likely not have glohally asymp- totically stable equilibrium point, and more importantly. (b) currently force the damping learned in latent space to be positive definite. However, these systems are not common in practice as friction and other dissi- pation mechanisms suss-spring-frrittion system). A possible remedy could be to relax the positive definiteness of the damping matrix in the latent space, allowing for zero damping. This would allow the method to work for systems without dissipation, such as conservative systems. Examples of such systems include mass-spring system without damping, the n-body problem, etc.
omnipresent, and the proposed method can learn very small damping values (c.g.. the
Systems with discontinuous dynamics. The proposed method might underperform for systems with highly discontinuous dynamics, such as systems with impacts, friction, or other discontinuities. In these cases, the latent dynamics might not capture the real system's 's behavior accurately, and the control performance of feedforward feedback will very likely be worse than pure feedback. Again, the method should be able to capture local behavior well. possible remedy for learning global dynamics could be to augment the latent dynamies with additional terms that capture the discontinulties, such as contact and friction models (e.g.. stick-slip friction).
Systems with multiple equilibrium points. The original system having multiple equilibria confficts with the stability assumptions underlying the proposed CON latent dynamics. In this case, as, for example, seen on the pendulum+ rendulum-friettoo and double pendulum + friction results, the method might work locally but will not be able to capture the global behavior of the system. A possible remedy could be to relax the global stability assumptions of the CON network. For example, the latent dynamies could be learned in the original coordinates of CON while allowing W also to be negative definite. This would allow the system to have multiple equilibria & attractors. Examples of such systems include a robotic arm under gravity, pendula under gravity, etc.
Systems with periodic behavior. The proposed method will likely not work well for systems with periodic behavior, as they do not have a single, attractive equilibrium point. Examples of such systems include mass- spring system with periodic external force, a pendulum with periodic external force, some chemical reactions, etc. Again, it is likely possible apply the presented method learning local behavior (i.c., not completing the full orbit). A possible remedy could be to augment the latent dynamics with additional terms that capture the periodic behavior, such substituting the harmonic oscillators with Van der Pol oscillators to establish limit eycle a supercritical Hopf bifurcation.
F.3 Systems for which we would not expect the proposed method to work
Nonholonomic systems. The proposed method likely would not work well for nonholonomic systems, both structure (e.g.. physical constraints) and stability characteristics would not be shared between the real system and the latent dynamics. Examples of such systems include vehicles, a ball rolling surface, and many mobile robots.
Partially observable and non-markovian systems. As the CON dynamics are evaluated based on the latent position and velocity encoded by the observation of the current time step and the observation-space velocity, we implicitly assume that the system is (a) fully observable and (b) fulfills the Markov property. This assumption might not hold for partially observable systems, such as systems with hidden states or systems with delayed observations. Examples of such cases include settings where the system is partially occluded or in situations without sufficient (camera) perspectives covering the system. Furthermore, time-dependent material properties, such as viscnelasticityy or hysteresis, that are present and significant in some soft robots and deformable objects are not captured by the method in its current formulation.
44

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims and contributions made in the abstract and introduction are all supported by theoretical and/or experimental results included in the paper.
Guidelines:
The answer NA mcans that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The known limitations of the proposed method are listed and presented multiple places in the manuscript: the last paragraph of the introduction, Section 6 and an extended version in Appendix
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results these assumptions (e-g., independence assumptions, noiseless settings, model well-specification.
violations of
asymptotic approximations only holding locally). The authors should reflect on how these assumptions might he violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
applicaalle, the authors should discuss possible limitations their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't it acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: All equations, Theorems, and Lemmas are numbered and cross-referenced. In cach Theorem/Lemma, we clearly state the assumptions under which the proof is valid (e.g., positive definite matrices for the ISS stability proof). All Theorems are included in the main paper: for the global asymptotic stability proof, we directly detail the proof in the main paper, and for the ISS proof, we provide a sketch with the full proof appearing in the Appendix. We also include auxiliary Lemmas in the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
45

---

All the theorems, formulas, and proofs in the paper should be numbered and cross roos-reffreeced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by
Tormal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: All the code for the experiments has been open-sourced on GitHub. Furthermore, we provide detailed description the implementation details in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can accomplished various ways. example, contribution novel architecture, describing the architecture fully might suffice, if the contribution is specific model and empirical evaluation, may be necessary to either make possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data often good way accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to hosted model (e.g.. in the checkpoint, other means that appropriate to the research performed.
of a large language model), releasing of model
While NeurlPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily new algorithm, the paper should make clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is way to access this model for reproducing the results way to reproduce the model (c.g.. with open-source dataset instructions for how construct the dataset).
model (e.g.. large language model), then there should either be
(d) We recognize that reproducibility may be tricky some cases, which case authors are welcome to describe the particular way they provide glosed-source models, it may be that access to the model is limited in some way (e.g.. registered users), hut it should be possible for other researchers have some path reproducing verifying the results.
reproducibility. In the case of
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instructions faithfully reproduce the main experimental results, described supplemental material?
Answer: [Yes]
The code associated with paper is available on Gut It allows the user to generate the datasets, run hyperparameters selection, train the models, and generats result plots based on training checkpoints.
Justification:
Guidelincs:
The answer NA means that paper does include experiments requiring code.
Please see the NeurlPS code and data submission guidelines httt:s::/nips.cc///
guides des /CodeSubm ississsinn oniionnoliiii for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g.. for a new open-source benchmark).
attpp:://ithub... phip wtt.///ccveiigggg. coupled-oscilleen aetworks pixels
46

---

The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https://nips. cc/
publ http///lcc/cccceuuaeeeeeee Swwmm/cccccccssscccccee ssionnollssiioooooo (omp)llcyyyy.. for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Dettils
Question: Does the paper specify training and test details (e.g.. data splits, hyperparameters.
how they were chosen, type of optimizer, etc.) necessary to understand results?
Answer: [Yes]
Justification: In Appendix C we include implementation details such the used libraries, algorithms, optimizers, and evaluation procedures. All hyperparameters (c.g., learning rate, loss weights, weight decay, etc.) can be found in the code on GitHub (specifically, in the sweep_generln e_dynami ca_autooencoerr py Python script that is placed in the e stamples/eee sweeppins/weppigg folder).
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, the appendix, as supplemental material.
Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate informa- tion about the statistical significance of the experiments?
Answer: [Yes]
Justification: We all experiments with various initializations (i.e., different random seeds), and each result table/plot is accompanied by description of how the variability of results is captured.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors). It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error hars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables or figures rates).
symmetric bars that would yield results that are of range (e.g. negative
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: In Section C.8 we report details about the necessary compute for performing the experiments reported in this paper.
Guidelines:
The answer NA means that the paper does not include experiments.
47

---

The paper should indicate the type of compute workers CPU or GPU, internal cluster. or cloud provider, including relevant memory and storage.
The paper should provide amount compute required for each the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn' It make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethies htttps://neurrt http:///............aaaaaaaaaaaa
Answer: [Yes]
Justification: The authors have reviewed the NeurIPS Code of Ethics, and this research conforms, in every respect, with this code.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
the authors answer No, they should explain special circumstances that require deviation from Code of Ethics.
The authors should make sure preserve anonymity (e.g., if there is special consideration due to laws regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [NA]
Justification: This paper primarily involves fundamental research, and the presented application of predicting and controlling the future evolutions of dynamical systems does not directly have any societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed. the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g.- disinformation, generating fake profiles, surveillance). fairness considerations (e.g., deploy- ment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications. the authors should point it out. For example, is legitimate to point out that an improvement the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it not needed point out that generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is heing used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitor- ing misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place responsible release data or models that have a high risk for misuse (e.g. pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This work by itself does not pose any risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
48

---

We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
The code accompanying this submission is original. While we do rely on common open-source 3rd-party packages (e.g.. JAX, flax, diffrrx, etc.), we clearly document these dependencies in the requi rements..... txt file the accompanying code archive. Furthermore, we leverage the datasets that are part of the DeepMind Humiltonian Dynamics Suite 26 and have been open-sourced with an Apache 2.0 license.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from particular source (e.g., website), the copyright terms of service of that source should provided.
If assets are released, the license, copyright information, and terms of use the package should be provided. For popular datasets, datasets.ppape |papersw thooee. com/ cat aset has curated licenses for some datasets. Their licensing guide help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
this information not available online, the authors are encouraged reach out to asset's creators.
13. New Assets
Question: Are alongside the assets?
introduced in the paper well documented and is the documentation provided
Answer: [Yes]
Justification: We release the necessary code to generate the datasets that we used in this paper alongside the submission.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate details of the dataset/coulemmodel as part of their sub- missions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, rememher to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This research did not involve any crowdsourcing experiments or trials with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
attps: http://ttttttccteee wtt.://..com....eeeeeeaee Stmi www.mmmmminni... wtw.miitiiiiiddddiiieeeee..
49

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This research did not involve any trials with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country which research conducted, approval (or equivalent) may be required for any human subjects research. you obtained IRB approval. you should clearly state this the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
initial submissions, not include any information that would break anonymity applica- ble), such institution conducting the review.
50


Reviews:
[{"review_id": "zL1z9j5n3H", "reviewer": "Reviewer_W9L3", "summary": "- This paper proposes Coupled Oscillator Networks (CONs); Networks consisting of coupled one-dimensional dampened harmonic oscillators, coupled through a neuron-like connection.\n- It is shown that under some constraints, the unforced coupled oscillators have a single, globally stable equilibrium and the forced system is globally ISS stable.\n- As the CON has no closed-form solution, the authors propose to split the solution in one decoupled linear, analytically solvable term and a residual non-linear term which is. For this system, a closed-form approximation is derived (CON-CFA).\n- High dimensional observations are encoded via a VAE into a compressed latent space that is trained to reconstruct future observations that are predicted via a CON latent dynamics model. This system is trained jointly.\n- “CONs are an ideal fit for learning latent dynamics as they guarantee that the latent states stay bounded.”\n- The approaches (CON with two different latent dimensions,S and M, and CON-CFA) are evaluated on a simulated soft-robots dataset and benchmarked against other latent-space dynamics methods based on neural ODEs ore autoregressive models.\n- When predicting system state, CON-M performs on par with sota methods, CON-S and CON-CFA perform almost as good, while cheaper computationally.\n- For control, potential-shaping strategies are combined with PID, resulting in a latent-space control law. A forcing decoder is then trained to to predict control inputs. This approach strongly outperforms a simple latent space PID controller.", "questions": "- Have you tried the method on other (non-soft) robots? Do you have any intuition on how your method would perform on environments with many contact forces?\n- Have you considered non-physical systems?", "limitations": "Many limitations of the method are clearly stated.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- Very interesting and (to my knowledge) novel approach of controlling robots in a latent dynamics space, modeled through a network of coupled oscillators with provable stability guarantees.\n- Extensive proofs and detailed descriptions of implementation.\n- Promising experimental evaluation.\n- Very concise and complete presentation of their approach.", "weaknesses": "- Limited experimental evaluation on a single simulated soft robot.\n- It is unclear how the obtained results generalize to other robots, in particular with less smooth contact dynamics.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer W9L3 (R2)\n\nWe thank the Reviewer for the careful reading and the encouraging comments.\nIn the following, we will respond to the questions raised by the reviewer, which also relate to the weaknesses mentioned by the reviewer.\n\n## Application of CON to non-soft robots\n\n> R2: Have you tried the method on other (non-soft) robots?\n\nWe thank the Reviewer for inquiring about results on non-soft robots. We provide additional results on various mechanical/robotic systems in the global rebuttal.\n\n> R2: Do you have any intuition on how your method would perform on environments with many contact forces?\n\nFirst, we want to stress that we have not conducted any experiments with contact-rich systems yet and reserve this interesting challenge for future work.\nWe hypothesize that the highly discontinuous dynamics of contact-rich systems could present a challenge for the method in its present form.\nStill, we envision that the method could be augmented to handle such systems. For example, we could add stick-slip friction or similar mechanisms to the CON dynamics to increase their expressiveness while maintaining the physical structure.\nFurthermore, contact-rich systems often have multiple equilibria, which would be a violation of the stability conditions introduced in this work. As touched on in Section 6 of the paper, these strong global stability guarantees would probably need to be relaxed.\n\n## Application of CON to non-physical systems\n\n> R2: Have you considered non-physical systems?\n\nWe thank the Reviewer for their interest in this topic. While we do think that the CON / CFA-CON models could be potentially useful in other applications where strong stability guarantees are required, we strive to focus in this paper on learning the dynamics of physical/mechanical systems as we can here leverage (a) the shared stability characteristics between the original and the latent space systems, and (b) exploit the mechanical structure of the CON model for control.\nTo make this focus clear, we will, in the final paper, modify the sentence starting on line 69 to (with changes marked in **bold**):\n\n_We resolve all the above-mentioned challenges by proposing Coupled Oscillator Networks (CONs), a new formulation of a coupled oscillator network that is inherently Input-to-State Stability (ISS) stable, **for learning the dynamics of physical systems,** and subsequently exploiting its structure for model-based control in latent space._\n\nThis work was specifically focused on how to impose structural biases that enforce second order Lagrangian structure and stability properties. However, we believe that the proposed model can be used more broadly, beyond mechanical systems, whenever strong stability guarantees are needed. This will be the focus of future research.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "2XKnZNzCqI", "reviewer": "Reviewer_Ys99", "summary": "In this paper, the authors introduce Coupled Oscillator Networks (CONs), extending the idea of coRNNs but with the difference that the generalized force is symmetric, enabling the definition of a potential energy expression which can be exploited for energy shaping based control methods. It is proven, that this model is asymptotically stable and, given an input, input-to-state-stable. It is shown that CONs provide similar performance compared to other state of the art methods. Then, the proposed network is used together with a variational autoencoder to learn a low-dimensional representation of the dynamics of a soft robot. Finally, a PID controller with an energy shaping part is designed to steer the soft robot to a given desired state. The main contribution is the new network structure which allows us to exploit the potential energy expression.", "questions": "- Can you add the inference time for the different methods in Table 1? If that’s not possible, could you give some general comments on the inference time of CONs?\n- Figure 3 visualizes that the controller with feed-forward part leads to heavy oscillations in the systems. Is that due to a poorly tuned controller or do you see the reason in the CON model?\n\nUPDATE: New score after discussion", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- It was a pleasure to read the article as it is nicely written, and all steps are explained in detail and clarity.\n- Computing the solution of the CON by splitting the dynamics into a linear part (where the solution can be computed in closed form) and a nonlinear part is an original idea.\n- The new structure of the proposed network that allows us to exploit energy-based control methods open doors for many interesting applications and might be a significant contribution.", "weaknesses": "- Theorem 1 seems to be trivial as it elaborates on the (nonlinear) coupling of passive systems (damped mass-spring systems) which is always passive since it does not “produce” any energy and thus, is inherently stable. (Same for Theorem 2). In this sense, the theoretical contribution seems limited. Is there anything I missed?\n- The statement in Line 266 “PID controller has several well-known drawbacks, such […] steady-state errors (in case the integral gain is chosen to be zero)” is misleading as in this case we would say it’s a PD controller.\n- The main contribution seems to be the new network structure that allows to exploit the potential energy for energy shaping methods. The VAE and controller are existing methods. Thus, a more detailed elaboration on the performance and limitation of CONs would be beneficial as it is evaluated for soft robotics data sets only.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer Ys99 (R1)\n\nWe thank the Reviewer very much for the kind words, for their interest in our research activities, and for the very insightful comments that they have provided.\nBecause of space constraints, we respond to the request for inference times and questions about the oscillations of the controller in the global rebuttal.\n\n## Is the stability of CON trivial?\n\n> R1: Theorem 1 seems to be trivial as it elaborates on the (nonlinear) coupling of passive systems (damped mass-spring systems), which is always passive since it does not “produce” any energy and thus, is inherently stable. (Same for Theorem 2). In this sense, the theoretical contribution seems limited. Is there anything I missed?\n\nWe thank the Reviewer for their question and for allowing us to elaborate further on this topic. In short, we differ with the Reviewer's technical arguments, and, consequently, with their conclusions on the triviality of Theorem 1. \nIn essence, we recognize that the Reviewer's arguments are well funded in the context of port-Hamiltonian theory. However, substantial extra properties/assumptions are needed to yield to the Reviewer's conclusion. Thus, the PH way of proving Theorem 1 would be interesting and proper, but it would definitely not result in the trivial proof the Reviewer is hinting at.\nMore precisely, we point out two main points of technical disagreement, which we list below.\n\n### (A): Passivity does not imply Stability\n\nPassivity, even in the scalar case, does not imply (global) stability. Additional assumptions on the (global) convexity of the energy landscape are needed, which would be hard to impose. As a trivial counterexample, think of a damped mass placed atop an infinite hill. \nTo further illustrate this point, we propose two simple passive-but-not-globally-stable CONs that deviate slightly from the assumptions in Theorem 1 and 2.\n\n1. Take a passive CON as stated in Eq. (2) in the scalar case with $K=1$ and $D=0.4$. Conflicting with the assumptions in Theorem 1, we select $W=-5$. Now, depending on the choice of $b$, this scalar system either has 1, 2, or 3 equilibria (see Fig. R1 in global response). All global stability guarantees are lost if the system has multiple attractors. \n2. We study a scalar CON with $K=-1$, positive damping $D=0.4$, and $W = 0$, $b=0$ resulting in the EOM $\\ddot{x} + K x + D \\dot{x} = \\tau$. We take the passive output $o = \\dot{x}$, and we prove passivity according to Def. 6.3 (Khalil, 2001) using the storage function $V(x) = K x^2 + \\dot{x}^2$: $\\dot{V}(x, \\tau) = \\tau \\dot{x} - D \\dot{x}^2 \\leq \\tau o. $\nHowever, this system is unstable, as it can be easily assessed by looking at the linearization at the equilibrium. Fig. R2 reports the globally repulsive vector field, which exhibits a single unstable equilibrium at the origin.\n\nTherefore, we conclude that passivity is not sufficient to prove global asymptotic stability.\n\n### (B): Stable harmonic oscillators do not imply Stable Networks\n\nEven if the individual systems (i.e., the harmonic oscillators) are stable, care needs to be taken when coupling them not to create an unstable network.\n\nWe illustrate this with the following example: consider a CON of dimensionality two with $K = [[1.0, -1.4],[-1.4, 1.0]]$, $D = \\mathrm{diag}(0.4, 0.4)$ $W = \\mathrm{diag}(3, 3)$, \nIt can be easily shown that the oscillators individually with the EoM $\\ddot{x}_i + k x_i + d \\dot{x}_i + \\tanh(w x_i)= 0,$ where $k=1$, $d=0.4$, and $w=3$ are globally asymptotically stable. \n\nHowever, the linear stiffness matrix is negative definite, and therefore, the system is not globally asymptotically stable. This is illustrated in Fig. R3 of the global response.\n\n## PID vs. PD\n\n> R1: The statement in Line 266 “PID controller has several well-known drawbacks, such […] steady-state errors (in case the integral gain is chosen to be zero)” is misleading as in this case we would say it’s a PD controller.\n\nWe thank the Reviewer for their comment and for pointing out the mistake. We agree that this sentence is indeed badly written and confusing. In the final version of the paper, we will remove the subsentence \"_steady-state errors (in case the integral gain is chosen to be zero)_\".\n\n## Performance and Limitations of CON\n\n> R1: The main contribution seems to be the new network structure that allows to exploit the potential energy for energy shaping methods. The VAE and controller are existing methods.\n\nWe thank the Reviewer for their comment. We want to stress that the proposed network imposes a beneficial inductive bias for conserving global stability and ISS. Therefore, we consider the two proofs (i.e., GAS and ISS) important contributions of the paper.\nFurthermore, we also regard the closed-form approximation of the CON dynamics as an important tool for deploying oscillator networks in practice.\n\n> R1: Thus, a more detailed elaboration on the performance and limitation of CONs would be beneficial as it is evaluated for soft robotics data sets only.\n\nFor this rebuttal, we have performed additional experiments involving non-soft-robotic datasets. They demonstrate that the CON network can also learn the latent dynamics of other mechanical systems, such as a mass spring, a pendulum, and a double pendulum with friction, effectively and with SOA performance. We refer to the global rebuttal for more details.\n\nIn totality, the results provided in the paper and in the rebuttal show that the performance of the CON model is on par with other SOA methods while adding physical structure and stability guarantees.\nAs detailed in Section 6.2, the strong assumptions needed to provide global stability guarantees are the primary limitation of the presented CON network, making it unsuitable for applications where complex attractor dynamics (e.g., multiple attractors, strange attractors, etc.) are required.\nRelaxing the stability assumptions could make the CON network also suitable for these applications.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
# Global Rebuttal

## Performance of CON on non-soft-robotic datasets (R1 & R2)

> R1: Thus, a more detailed elaboration on the performance and limitations of CONs would be beneficial as it is evaluated for soft robotics data sets only.

> R2: Have you tried the method on other (non-soft) robots?

We thank both reviewers for their interest in the performance of CON on non-soft-robotic datasets.
For this rebuttal, we compared the performance of CON against the baseline method on three additional mechanical, non-soft-robotic datasets: a mass-spring with friction (_M-SP+F_) (i.e., a damped harmonic oscillator), a single pendulum with friction (_S-P+F_), and a double pendulum with friction (_D-P+F_). These datasets are based on an interesting publication by Botev et al. (2021) [25], which appeared in the _NeurIPS 2021 Track on Datasets and Benchmarks_ and benchmarks various models for learning latent space dynamics.

The results, which we will refer to as Table R1 of the global response PDF, show that the NODE model slightly outperforms the CON network on the _M-SP+F_ and _S-P+F_. However, as the datasets do not consider system inputs, we can remove the input mapping from all models (e.g., RNN, GRU, coRNN, CON, and CFA-CON). With that adjustment, the CON network has the fewest parameters among all models and particularly two orders of magnitude less than the NODE model. Therefore, we find it very impressive that the CON network is roughly on par with the NODE model.
For the _D-P+F_ dataset, we can conclude that the CFA-CON model offers the best performance across all methods.
Finally, most of the time, the CON & CFA-CON networks outperform the other baseline methods that have more trainable parameters.

## Inference time of CON (R1)

> R1: Can you add the inference time for the different methods in Table 1? If that’s not possible, could you give some general comments on the inference time of CONs?

We thank the reviewer for their question about the inference time of the various methods. We note that the number of training steps per second of all methods included in Table 1 was already reported in the original submission in Table 4 of Appendix D.

For this rebuttal, we performed additional evaluations of the inference time (i.e., without computation of loss function and gradient descent) of the various models and report the results in Table R2 of the global response PDF.

## Oscillations of the controller with FF term. (R1)

> R1: Figure 3 visualizes that the controller with a feed-forward part leads to heavy oscillations in the systems. Is that due to a poorly tuned controller or do you see the reason in the CON model?

We thank the Reviewer for their question and for raising the topic. When tuning the gains PID-like controllers, a trade-off naturally exists between transient behavior (e.g., oscillations and overshooting) and response time.
In this case, we chose gains that minimized the response time but allowed for stable behavior. The oscillations are caused by a combination of (a) the underdamped nature of the system and (b) the magnitude of the proportional term.
Importantly, to have a fair comparison, we kept the gains of the feedback controller the same for both the _P-satI-D_ and _P-satI-D + FF_ cases. A higher proportional term is beneficial for the response time (and the performance) of the _P-satI-D_,
while it leads to overshooting and oscillations in the _P-satI-D + FF_ case. We stress that this is not an inherent problem of the feedback controller but can be mitigated by tuning the feedback gains differently.
For this rebuttal, we tuned a controller with reduced proportional and increased damping term and the results, included as Fig. R4 in the global response PDF, show that the oscillations and overshooting are both significantly reduced.

Meta-review: The paper presents a very interesting combination of traditional control methods and machine learning, applied to robotic control (non-soft robot results were added at rebuttal time). Unfortunately there were only 2 reviews but both are very enthusiastic and I concur with their assessment. The authors should make sure that the final version reflects the clarifications provided in the rebuttal.
Decision: Accept (spotlight)