Learning Deep Input-Output Stable Dynamics
Authors: Ryosuke Kojima
Yuji Okamoto
Venue: NeurIPS 2022 Accept (2022)
Keywords: Dynamical systems
Input-output stability
Nonlinear system identification

Abstract: Learning stable dynamics from observed time-series data is an essential problem in robotics, physical modeling, and systems biology. Many of these dynamics are represented as an inputs-output system to communicate with the external environment. In this study, we focus on input-output stable systems, exhibiting robustness against unexpected stimuli and noise. We propose a method to learn nonlinear systems guaranteeing the input-output stability. Our proposed method utilizes the differentiable projection onto the space satisfying the Hamilton-Jacobi inequality to realize the input-output stability. The problem of finding this projection can be formulated as a quadratic constraint quadratic programming problem, and we derive the particular solution analytically. Also, we apply our method to a toy bistable model and the task of training a benchmark generated from a glucose-insulin simulator. The results show that the nonlinear system with neural networks by our method achieves the input-output stability, unlike naive neural networks. Our code is available at https://github.com/clinfo/DeepIOStability .


Learning Deep Input-Output Stable Dynamics
Ryosuke Kojima'
Yuji Okamoto" Graduate School Medicine
Graduate School Medicine
Kyoto University
Kyoto University Kyoto, 606-8501
Kyoto, 606-8501
koj iroouukk Be@kyoto-u.acc ac..
okamoto wwt.mooo.yu..222ttytoou.........
Abstract
Learning stable dynamics from observed time-series data is an essential problem in robotics, physical modeling, and systems biology. Many of these dynamics are represented as an inputs-output system to communicate with the external environ- ment. In this study, focus on onput- stable systems, exhibiting robustness against unexpected stimuli and noise. We propose a method to learn nonlinear systems guaranteeing the input-output stability. Our proposed method utilizes the differentiable projection onto the space satisfying the Hamilton-Jacobi inequality realize the input-output stability. The problem of finding this projection be formulated as a quadratic constraint quadratic programming problem, and we derive the particular solution analytically. Also, we apply our method to a toy bistable model and the task of training benchmark generated from glucose- insulin simulator. The results show that the nonlinear system with neural networks by our method achieves the input-output stability. unlike naive neural networks. Our code is available at https:: httpps://githb... www. httm///lnnft/DeppSSaai1111
Introduction
Learning dynamics from time-series data has many applications such as industrial robot systems [. physical systems and biological systems 3 Many of these real-world systems equipped with inputs and outputs to connect for each other, which are called input-output systems [5]]. For example, biological systems sustain life by obtaining energy from the external environment through their inputs. Such real-world systems have various properties such stability, controllability, and observability, which provide clues analyze the complex systems.
Our purpose is to learn a complex system with "desired properties" from dataset consisting of pairs of input and output signals. To represent the target system, this paper considers the following nonlincar dynamics:
4 = f(x) + G(x)u y M(x).
x00)= = To
where the inner state æ, the input u, and the output y belong to signal space that maps time interval [0, oc) to the Euclidean space. We denote the dimension of x, u, and as n, m, and 1, respectively.
Recently, with the development of deep learning, many methods to learn systems from time-series data using neural networks have been proposed [ By representing the maps (f, G, h) Eq. ① neural networks, complex systems
be modeled and trained from given dataset. However, guaranteeing that a trained system has the desired properties is challenging.
A naively trained system fits the input signals contained in the training dataset, but does not always fit for new input signals. For example, Figure I shows our preliminary experiments where we
Equul contribution.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

1.0
data
1.0
y: ouuubut
predict
0.5
0.5
0.0
0.0
1.0
1.0
U: Input
iniut
0.0
0.0
t:tme
t: time
(A)
(B)
Figure 1: (A) The model prediction of neural networks (B) the reaction of trained model. These results are min-max normalized.
naively learned neural networks G, in Eq. ! from imput and output signals. The trained neural networks provided small predictive errors for input signal in the training dataset (Figure (A) ) Contrastingly, the unbounded output signals were computed by the trained system with new step input signals (Figure T (B) The reason be expected that the magnitude (integral value) of this input signal was larger than that of the input signals in the training dataset.
The internal stability is known as one of the attractive properties that should be often satisfied in real-world dynamical systems. The conventional methods to train internal stable systems consisting of neural networks adopt the syapunov-based approaches [-II These methods focus on the internal with the
system, à = (ar in the input-output system I Thus, how to learn the entire system |
desired property related to the influence of the input signal is still challenging.
We propose a novel method learn dynamical system consisting neural networks considering the input-output stability. The notion of the input-output stability is often used together with the Hamilton- Jacobi inequality for controller design of the target system the field of control theory. The Mamilton-Jacohe inequality is one of the sufficient conditions for the input-output stability. The feature of this condition that the variable for an input signal u does not appear in the expression, i.e, we do not need to evaluate the condition for unknown inputs [I22. To the best of our knowledge, this is the first work that establishes a learning method for the dynamical systems consisting of neural networks using the Hamilton-Jacole inequality.
The contributions of this paper are as follows:
This paper derives differentiable projections to the space satisfying the Hamilton-Jacobi inequality.
This paper presents the learning method for the input-output system proven to always satisfy the Hamilton-Jacobi inequality.
This paper also provides loss function derived from Hamilton-Jacobi inequality. By combining this loss function with the projection described above, efficient learning can be expected.
This paper presents experiments using two types of benchmarks to evaluate our method.
2 Background
This section describes the Lg stability, standard definition of the input-output stability, and the Hamilton-Jacohi inequality.
First, we define the L2 stability of the nonlinear dynamical system I If the norm ratio of the output signal to the input signal bounded, the system is La stable and this norm ratio called the Ca gain. The L2 norm on the input and output signal space is used for the definition of the L2 stability. To deal with the L2 stability on the nonlinear system ( we assume that the origin x = 0 of the internal system å f(x) is asymptotically stable equilibrium point and (())) 0. Since the translation of

---

any asymptotically stable equilibrium points is possible, this assumption can be satisfied without loss of generality.
Definition The C2 norm is defined as || = So [|x()) If there exists y Z 0 and
function 3 on domain D C IR" such that
llullc, S rllullc, + 3(zo),
(2)
then the system 77 is Lg stable, where the function (3-) non non-negative and (((0) Furthermore, the minimum r satisfying is called the L2 gain of the nonlinear system 1
Next, we describe sufficient condition of the L2 stability using the Lyapunov function : R. where the function V is positive semi semi-definite, i.e., V(x) 2 0 and V(0) = The Hamilton-Jacok
inequality is input-independent sufficient condition.
Proposition (5 Theorem 5.5]). Let f be locally Lipschitz and let G and be continuous. If there exist constant 7 > O and a continuously differentiable positive semi-definite function : C R" + R such that
(VT(x)f(x) + 2-72 2 IG" 110=0)11== 00, Vx E D
(3)
then the system 7 is L2 stable and the La2 gain is less than equal to 7, This condition is called the Hamilton- Jacobi inequality.
The above proposition can be generalized to allow the more complicated situations such as limit-cycle and bistable cases, where the domain D contains multiple asymptotically stable equilibrium points. The equilibrium point assumed in this proposition can be replaced with positive invariant sets by extending the L2 norm [I3]. Furthermore, by mixing multiple Lyapunov functions, this proposition can be generalized around multiple isolated equilibrium points.
3 Method
The goal of this paper is to learn the L2 stable system represented by using neural networks G,,, The Hamiltonian-Jacob inequality, which implies the L2 stability, is expressed by (f, h). We present method project (f, h) onto the space where the Hamilton-Jacobi inequality holds.
3.1 Modification of nonlinear systems
Supposing fn R R., G. R.= Raxm and hn R. -> RR. a triplet map (Jns Gn hn) denote as nominal dynamics. Introducing a distance on the triplet maps. the nearest triplet satisfying the Hamilton-Jacobi inequality from the nominal dynamics (Jn- Gn hn) is called modified dynamics (fm. Gm, hm). The purpose of this section is to describe the modified dynamics (fm, Gms hm) associated with the nominal dynamics (fn, Gn, Ann) by analytically deriving a projection onto the space satisfying the Hamilton-Jacobi inequality.
The problem of finding the modified dynamics (fm. Gm, hm) is written as a quadratic constraint quadratic programming (QCQP) problem for the nominal dynamics (fn, Gn. hn). Since there is generally no analytical solution for QCQP problems, we aim to find the particular solution by adjusting the distance on the triplets.
To prepare for the following theorem, define the ramp and the clamp functions as
a, x S a
0,
x 0 C(x; 0, b)
R(x)
v,
a < x b.
I,
x > 0
b
x> > b
and define the Hamilton-Jacobi function
HJ G. VVT
G VVI + 1141
(4)
2-
where V a given positive definite function. A way to design V is to determine the desired positive invariant set and design the increasing function around this set. For example, if the target system has a unique stable point, æ can be regarded the positive invariant set and the increasing function V can be designed V(x) 1x2
3

---

Theorem 1. Consider that the following optimal problem:
(1 (1- k)
k
minimize Jm Gm Pem
IV 11 -Sall + 27 |I
Gnll2 1
|VV]2 Mm
hnll
(5a)
subject to
(J Gm: Am) S 50,
(5b)
where & E [0, 1]. The solution of 5 is given by
fm = fn |AV|2
R (V, + VV,
Gm Gn
(o(-v%_.2..)) PrCn.
Vora
hm = =Vo((roi V :k?, hin
where
Vj 2 ST
Vah - I
2/2 27
Pv
VVVVT |
Proof: See Appendix A
The objective function of 5a) is new distance between the nominal dynamics (Sn, Gn, hn) and the modified dynamics (Sms Gm, h. This new distance allows the derivation of analytical solutions by combining three distances of f. G, and h.
Focusing on the objective function (5a). the constant & represents the ratio the distance scale of f to G and h. When A: 0, the result of this problem (5) are consistent with the projection of the conventional method that guarantee internal stability 9 As the constant & converges 1, the modification method of Theorem approaches G. and hm to Gn and hn. respectively. In this case, the objective function (5a) becomes the distance between fm to fn- Therefore, the following corollary is satisfied.
Corollary 1. The solution of
minimize fm subject to
[|5m - Snll
(6a)
HJ(mm, Gn, l(n) S 0,
(6b)
given by
fm fn
| R (HJ(fn, Gn hn)) VV.
Proof: This solution is easily derived from Theorem
Corollary I derives solution of linear programming problem rather than QCQP problems. Re-
placing the Hamilton-Jacobl function HJ Gn, /t)) with the time derivative positive definite function VVTF, this corollary matches the result of the conventional study [9]
When the map hm is fixed as hn, a similar solution as Theorem I is derived. Although Corollary is proved by changing k. the modified dynamics with the fixed hn are not derived. We reprove this modified dynamics in a similar way to Theorem T
Corollary 2. Consider the following problem:
(1 k)
k Gm Gnll2 [VV| 1m -all + 2y
minimize Sm.Gum subject to where & E [0, 1). The solution of
(7a)
HJ(Sm, km)) 50,
(7b)
is given by
fm Sn
1
R (Vrn + V2 VV,
Gm Gn 1- vo(
Va
Vyn
PvG (8)
1

---

Contour line of Ellyy 9119
HJ(J,G,) 50
Reduced region
Gu hud
Go Am
Modifidation
Lan
[mm Gar Ami
Gam Fms
mm
VEII
(A)
(B)
(C)
Figure 2: Sketches of our method (A) minimizing the prediction error (the first term of the loss function ( in the blue region, (B) moving the nominal dynamics the blue region (the second term), and (C) reducing the blue region while keeping the same level of the prediction error (the last term).
where
Vn _ VVT fn + 11al11
Va - 27
Proof: See Appendix
3.2 Loss function
We represent (fn, hn) hn)) neural networks and denote (fm. Gm. fmm the C2 stable dynamics modified by Theorem Corollary [ or 2 Note that the modification depends on the candidate of C2 gain n. Figure (A) shows the sketch of this modification, where the blue region satisfies the Hamilton-Jacoh inequality.
Since the nonlinear system of the modified dynamics (fm- Gms hm) is represented as ordinary differential equations (ODEs) consisting of the differentiable functions fn,G,, and li, the techniques of training neural ODEs can be applied [4] Once loss function is designed, the parameters of neural networks in the modified system can be learned from given data.
Loss F. lll 0112,11 + ALns + ay?
where À and are non-negative coefficients.
The first term shows the prediction error of the output signal y Figure (A)). dataset D consists of tuples (Zo. u. y) where the initial value Tu, the input signal u, and the output signal The predicted output ý is calculated from 2g, u, and the modified dynamics (fm- Gm, hmn)
The second term aims to improve the nominal dynamics (Jn, Gn, hm) closer to the modified dynamics (Jm, Gmy hm) and is defined as
LH E. [g[RH(H (n) + ++)],
where is a positive constant Figure (B)). Since this term a form of the hinge loss, e represents the magnitude of the penalties for the Hamilton-Jacobi inequality. To evaluate this incquality for any x E D, we introduce a distribution of x over the domain D. In our experiments, this distribution is decided Gaussian distribution ((u,a²) where the mean μ is placed at the asymptotically stable point and the variance o² is an experimental parameter. Without this second term of the loss function, there degrees of freedom in the nominal dynamics, i.e., multiple nominal dynamics give the same loss by the projection, which negatively affects the parameter learning.
The modifying parameter can be manually designed for the application automatically trained from data by introducing the last term. This training explores smaller 7 while keeping the same level of the prediction error (Figure (C)).
5

---

4 Related work
Estimating parameters of a given system is traditionally studied as system identification. In the field of system identification, much research on the identification of linear systems have been done, where the maps (ff, G, h) in Eq. ① assumes to be linear [15]. Linear state-space models include identification methods by impulse response like the Eigensystem Realization Algorithm (ERA) [T by the state-space representation like Multivariable Output Error State sPace (MOESP) [7 and ORThogonal decomposition method (ORT) [8]. System identification methods for non-state-space models, unlike our target system, contain Dynamic Mode Decomposition with control (DMDec 19 and its nonlinear version, Sparse Identification of Noulinear DYnamics with control (SINDYe) [20 Also, the nonlinear models such the nonlinear ARX model 21 and the Hammerstein Wiener model 2e have been developed and often trained by error minimization. For example, the system identification method using piece-wise ARX model allows more complicated functions by relaxing the assumption of linearity [21]
These traditional methods often concern gray-box systems, where the maps (f. G. h) in Eq. !I are partially known [23]. This paper deals with a case of black-box systems, where f. G, and in the system (1 are represented by neural networks. Our method can be used regardless of whether all f, G, and h functions are parameterized using neural networks. So, our method can be easily applied to piriccatiooneeeddeent gray-box systems when the functions f. G, and h are differentiable with respect to the parameters.
Because the input-output system T can be regarded as a differential equation, our study is closely related the method of combining neural networks and ODEs 7 14] These techniques have been improved in recent years, including discretization errors and computational complexity. Although we used an Euler method for simplicity, we can expect that learning efficiency would be further improved by using these techniques.
The internal stability is fundamental property in ODEs, and learning system with this property using neural networks plays an important role in the identification of real-world systems (24.. In particular, the first method to guarantee the internal stability of the trained system has been proposed in [2] Furthermore, another method 2 extends this method to apply positive invariant sets, e.g. limit cycles and line attractors. Encouraged by these methods based on the Lyapunov function, our method further generalizes these methods using the Hamilton-Jacobi inequality to guarantee the onput-output stability.
In this paper, the Lyapunov function V is considered to be given, but this function can be learned from data. Since pair of dynamics (1] and V has redundant degrees of freedom, additional assumptions are required to determine V uniquely. In [ it realized by limiting the dynamics to the internal system and restricting V to convex function 26
Lyapunov functions are also used to design controllers, where the whole system should satisfy the stability condition. A method for learning such controller using neural nets has been proposed [ This method deals with optimization problems over the space that satisfies the stability condition, which is similar to method. Whereas solve the QCQP problem to derive the projection onto the space, this method uses a Satisfiability Modulo Theories (SMT) solver to satisfy this condition. The method has also been extended to apply unknown systems
28
Although this paper deals with deterministic systems, neural networks for stochastic dynamics with the variational inference is well studied 6 31) Guaranteeing the internal stability of trained
stochastic dynamics is important for noise filtering and robust controller design ITC
5 Experiments
We conduct two experiments to evaluate our proposed method. The first experiment uses a benchmark dataset generated from nonlinear model with multiple asymptotically equilibrium points. In the next experiment, we applied our method to biological system using simulator of the glucose-insulin system.

---

5.1 Experimental sefting
Before describing the results of the experiments, this section explains the evaluation metries and our experimental setting.
For evaluation metrics, define the root mean square enror (RMSE) and the average C2 gain (GainIO) for the given input and output signals as follows:
N
1 N
N
1 24 C2
RMSE A
1u 2=1
ds vill2a
GainlO
A
N
where N is the number of signals in the dataset. uz(+) and (-) are the input and output signal the i-th index, respectively. The prediction signal (:(-) is computed from w.(-), the trained dynamics, and the initial state. Note that the integral contained in the Le norm is approximated by a finite summation. The RMSE is a metric of the prediction errors related to the output signal, and the GainIO is metric of the property of the La stability. Whether the target system satisfies does not satisfy the L2 stability, the GainlO with a given finite-size dataset can be calculated. The GainIO error is defined by the absolute error between the GainlO of the test dataset and that of the prediction.
In our experiments, 90% of the dataset used for training and the remaining 10% is used for testing. We retry five times for all experiments and show the mean and standard deviations of the metrics.
For simplicity in our experiments, the sampling step AL for the output y is set as constant and the Euler method is used to solve ODEs. To is put at an asymptotically stable point for each benchmark and is known. In this experiment, to prevent the state from diverging during learning of dynamies. the clipping operation is used so that the absolute values of the states are less than ten.
For comparative methods, piece-wise ARX 2 In the method of vanilla neural networks, the maps (f, G, h) in the nonlinear system ① is represented by using three neural networks, i.e., this method is consistent with method used in Figure | To determine the hyperparameters of comparative methods except for neural networks, the grid search used. Note that these comparative methods only consider the prediction errors.
vanilla neural networks, ARX, ORT 18. MOESP
and
For training each method with neural networks, an NVIDIA Tesla T4 GPU was used. Our experiments are totally run on 20 GPUs over about three days.
5.2 Bistable model benchmark
The first experiment is carried out using a bistable model, which is known as a bounded system with multiple asymptotically equilibrium points. This bistable model defined as
i = ==1(1 - x²) +
z00)- 11 x.
The internal system of this model has two asymptotically stable equilibrium points 1, 1..
We generate 1000 input and output signals for this experiment. To construct this dataset, we prepare input signals using positive and negative pulse wave signals whose pulse width is changed random. The input and output signals on the period [o, 10] are sampled with an interval At = 0.1. In this benchmark, set the number of dimensions of the internal system as and use fixed function (xx) = min((x 1)2, (x + 1)²), a mixture of the two positive definite functions.
In the result of our experiments, we name the proposed methods modified by Theorem | Corollary and as DIOS-fgh, DIOS-f, and DIOS- -g, respectively. In these methods, the parameters of our loss function set as À and 0.01. Also, DIOS-fgh+ uses À 0.01 and a 0.01 under the same conditions as DIOS -igh. For this example (1000 input signal), it took about hour using I GPU to learn one model training.
The results of the RMSE and the GainlO error this experiment are shown in Figure Figures (A) and (B) demonstrate that a piece-wise ARX model (PWARX) gives very low RMSE but its Gainlo error is high. Our proposed methods achieve a small GainIO error while keeping the RMSE. Note that linear models such as MOESP and ORT only approach one point although the bistable model has two asymptotically stable equilibrium points. Figures 3 (C) and (D) display the effect of dataset sizes. When the datasct size was varied to 100, 1000, and 3000, the result of the larger dataset provided
7

---

GainIo error
RMSE
Vanilla
DIOS-f
DIOS-fg
DIOS-fgh
DIOS-fgh+
ARX
MOESP
ORT
PWARK
Vanilla
DIOS-f
D10S-fg
DIOS-fgh
DiOS-fgh+
ARX
MOESP
ORT
PWARX
(B)
2.0
0.8
Gainlo error
1.5
RMSE
0.4
1.0
0.5
0.0
0.0
Vanilla
DIOS-fgh fgh
DIOS-fgh+
Vanilla
DIOS-fgh
DIOS-fgh+
Vanilla
DIOS-fgh fih
DIOS-fgh+
Vanilla
DiOS-(gh fggh
DIOS-fght figh+
Vanilla
DIOS- fgh
pIOs-fgh+
Vanilla
DIOS-(gh
DIOS-fgh+ fgh++
Data size
Data size
Data size
Data size
Data size
Data size
100
1000
3000
100
1000
3000
(C)
(D)
Figure 3: Results the bistable model benchmark. The upper part shows (A) the RMSE and (B) the GainIO error of the vanilla neural networks (gray), our proposed methods (red), and the conventional methods (blue). The lower part shows (C) the RMSE and (D) the GainIO error for the different sizes of the datasets.
Vanilla
DIOS-fgh
DIOS-fgh+
0.23
0.15
0.15
0.15
0.2
0.02
-3.05
00.0
0.05
-325
-0.2
True
True
Data
100
1000
Data size: 3000
Figure 4: The sketch (x) displaying the internal dynamics trained from the bistable datasets with the different sizes. The sign the bistable model is shown the bottom of each figure.
the smaller RMSE for all methods. The vanilla neural networks do not consider the La gain, so the GainIO error large in the case of the low RMSE.
Figure shows the relationship between x and fz) in the trained system in ① compare the vanilla neural networks, DIOS-fgh and DIOS-fght. DIOS-fgh and DIOS-fgh+ successfully find two stable points, i.e., the trained function (x) had roots of two stable points x = +11 and an unstable point x O. Especially, these stable points were robustly estimated in the trained system using our presented loss function (DIOS-fgh+) even when the dataset size was small. The vanilla neural networks failed to obtain the two stable points in all cases.

---

Vanilla
Our method (DIOS-fgh+)
True
0.15
0.6
Y1: Giucoso
M
M
Input
0.05
5 02
0.00
0.0
500
1000
500
1000
500
1000
D: Time
Figure 5: The input and output signals of the glucose-insulin simulator and the predicted output.
Clucoss
2
0
1.5
Ys: Insulin
1.0
1.0
0.5
0.0
0.0
10
Inpuut
5
10
30
40
50
10
20
30
40
50
10
20
40
50
Time
(A) Vanilla
(B) Our method (DIOS-fght fgh+)
(B) True
Figure 6: The step reaction of the trained systems. The color of each line indicates the magnitude of the step input.
5.3 Glucose-insulin benchmark
This section addresses an example of the identification biological systems. We consider the task of learning the glucose-insulin system using simulator 32 construct responses for various inputs and evaluate the robustness of the proposed method for unexpected inputs. This simulator outputs the concentrations of plasma glucose 1/1 and insulin y2 for the appearance of plasma glucose per minute u. To determine the realistic input u, we adopt another model, an oral glucose absorption model 33
Using this simulator, 1000 input and output signals are synthesized for this experiment. The input and output signals are sampled with a sampling interval At = and 1000 steps for each sequence. In this benchmark,
set the number of dimensions of the internal system as six, fix positive definite function V(x) = x?2 and use our loss function with À = 0.001 and a = 0.001. Training one model for this examples (1000 input signal) tooks about 7.5hours using 1GPU. The hyperparameters including the number of layers in the neural networks, the learning rate, optimizer, and the weighted decay determined using the tree-structured Parzen estimator (TPE) implemented in Optuna (34)
The RMSE of vanilla and our proposed method are 0.0103 and 0.0050, respectively. So, from the perspective of RMSE, these methods achieved almost the same performance. Figure shswws input and output signals in the test dataset and the predicted output by vanilla neural networks and method (DIOS-fgh++) The L2 stability of the system using the vanilla neural networks is not guaranteed. Since the proposed method guarantees the L2 stability, the output signals of DIOS-fght bounded even if the input signals unexpectedly large.
To demonstrate this, we conducted an additional experiment using trained system. Figure 5 shows the transition of output behavior caused by the magnitude of the input signal changing from 2 to 10. Note that the maximum magnitude the training dataset one. In this experiment, Al was changed to 0.01 and the clipping operation was removed to deal with the large values of the state.

---

This result shows the output of the vanilla neural networks quickly diverged with an unexpectedly large input. Contrastingly, the output behavior of our proposed method always bounded. Therefore, we actually confirmed that our proposed method satisfies the La stability.
6 Conclusion
This paper proposed a learning method for nonlinear dynamical system guaranteeing the Le stability. By theoretically deriving the projection of triplet (f, G. h) to the space satisfying the Hamilton- Jacobi inequality, our proposed method realized the La stability of trained systems. Also, we introduced a loss function to empirically achieve smaller L2 gain while reducing prediction errors. We conducted two experiments to learn dynamical systems consisting of neural networks. The first experiment used nonlinear model with multiple asymptotically equilibrium points. The result of this experiment showed that our proposed method can robustly estimate system with multiple stable points. In the next experiment, we applied our method to a biological system using a simulator of the glucose-insulin system. It was confirmed that the proposed method can successfully learn proper system that works well under unexpectedly large inputs due to the Le stability. There is limitation that our method cannot apply the system without the L2 stability. Future work will expand the Le stability- based method to more generalized learning method by dissipativity and apply our approach of this study to stochastic systems.
Acknowledgements
This research was supported by JST Moonshot R&D Grant Number JPMJMS2021 and JPMJMS2024. This work was also supported by JSPS KAKENHI Grant H0.205 and CREST Grant Number JPMICR22D3,, Japan. This paper was also based on part of results obtained from project commissioned by the New Energy and Industrial Technology Development Organization (NEDO).
References
[1] Jan Swevers, Walter Verdonck, and Joris De Schutter. Dynamic model identification for industrial robots. IEEE control systems magazine, 27(5):58-1) 71, 2007.
[2] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 2016.
[3] Timothy S. Gardner, Diego Di Bernardo, David Lorenz, and James J. Collins. Inferring genetic networks and identifying compound mode of action via expression profiling. Science, 301 2003.
[4] Geoffrey Roeder, Paul K Grant, Andrew Phillips, Neil Dalchau, and Edwards Meeds. Efficient amor- tised bayesian inference for hicrarchical and nonlincar dynamical systems. Proceeding of International Conference on Machine Learning (ICML), 2019.
[5] Hassan K. Khalil, Nonlinear systems. Prentice-Hall, 3rd edition, 2002.
[6] Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured inference networks for nonlinear state space models. In Proceedings the AAAI Conference on Artificial Intelligence, page 2101-2109, 2017.
[7] Ricky T.Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems (NeurIPS), volume 31, 2018
18] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. In Proceeding of International Conference on Learning Representations, 2019.
[9] Gaurav Manek and J. Zico Kolter. Learning stable deep dynamics models. In Advances in Neural Information Processing Systems (NeurIPS), volume 32, 2019.
[10] Nathan I Lawrence, Philip Loewen, Michael Forbes, Johan Backstrom, and Bhushan Gopaluni. Almost surely stable deep dynamics. In Advances in Neural Information Processing Systems (NeuriPS), volume 33, 2020.
[11] Andreas Schlaginhaufen. Philippe Wenk, Andreas Krause, and Florian Dörfter. Learning stable deep dynamics models for partially observed or delayed dynamical systems. In Advances in Neural Information Processing Systems (NeurIPS), volume 34, 2021.
10

---

[12] Arjan van der Schaft. C2-Guin and Passivity Techniques in Nonlinear Control. Springer Publishing Company, Incorporated, 3rd edition, 2016.
[13] W.M. Haddad and V.S. Chellaboina. Nonlinear Dynamical Systems and Control: Lyapunov- Based Approach. Princeton University Press, 2008.
[14] Xinshi Chen. Review: Ordinary differential equations for deep learning. CoRR, 191009901.002 2019.
[15] Tohru Katayama. Subspace Method for System Identification. Springer, 2005.
[16] Jer-Nan Juang and Richard Pappa. An eigensystem realization algorithu for modal parameter identifica- tion and model reduction. Journal of guidance, control, and dynamics, 8(5):620-227, 1985.
[17] Michel Verhaegen and Patrick Dewilde. Subspace model identification part 1. the output- error state-space model identification class of algorithms, International journal of control. 56(5): 187-1210. 1992.
[18] Tohru Katayama, Hideyuki Tanaka, and Takeya Enomoto. A simple subspace identification method of closed-loop systems using orthogonal docomposition. IFAC Proceedings Volumes, 38(1)5112-517, 2005.
[19] Joshua L Proctor, Steven Brunton, and J Nathan Kutz. Dynamic mode decomposition with control. SIAM Journal on Applied Dynamical Systems, 15(1):142-161, 2016.
[20] Steven I. Brunton, Joshua I, Proctor, and J Nathan Kutz. Sparse identification of nonlinear dynamics with control (SINDYc). PPAC-PapersOnnline, 49(18):710-715 2016.
[21] Andrea Garulli, Simone Paoletti, and Antonio Vicino. survey on switched and piecewise affine system identification. In IFAC Proceedings Volumes, volume 45, pages 344-355. IFAC, 2012.
[22] Adrian Wills, Thomas B. Schon, Lennart Ljung, and Brett Ninness, Identification of hammerstein-winn models. Automatica, 49(1):70 81, 2013.
[23] Lennart jung. System identification. In Signal analysis and prediction, pages 163-173. 16rieer,, 1998.
[24] Martin Braun and Martin Golubitsky. Differential equations and their applications. Springer, 1983.
[25] Naoya Takeishi and Yoshinobu Kawahara. Learning dynamics models with stable invariant sets. In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pages 9782-9790, 2021.
[26] Brandon Amos, Lei Xu, and Zico Kolter. Input convex neural networks. In International Conference Machine Learning, pages 146-155 155, 2017.
[27] Ya-Chien Chang, Nima Roohi, and Sicun Gao. Neural lyapunov control. Advances in Neural Information Processing Systems (NeurIPS), volume 32, 2019.
[28] Vrushabh Zinage and Efstathios Bakolas. Neural koopman Iyapunov control. arXiv:2201.0500 050098 2022.
[29] Daniel Gedon, Niklas Wahlström, Thomas Schön, and Lennart Ljung. Deep state space models for nonlinear system identification. In Proceedings of the 19th IFAC Symposium on System Identification 2021
[30] Junyoung Chung, Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron C. Courville, and Yoshua Bengio. A recurrent latent variable model for sequential data. In Advonces in Neural Information Processing Systems (NeurlPS), volume 28, 2015.
[31] Justin Bayer and Christian Osendorfer. Learning stochastic recurrent networks. ArXiv, abs/1411.76100, 2014.
[32] A. De Gactano and O. Arino. A statistical approach the determination of stability for dynamical systems modelling physiological processes. Mathematical and Computer Modelling, 31(4):41-51, 2000.
[33] Chiara Dalla Man, Robert Rizza, and Claudio Cobelli. Meal simulation model the glucose-insulin system. IEEE Transactions on Biomedical Engineering, 54(10): 17499 2007.
[34] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A next-generation hyperparameter optimization framework. In Proceedings the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.

---

Checklist
The checklist follows the references. Please read the checklist guidelines carefully for information how to answer these questions. For each question, change the default [TODO] to [Yes], [No], , or [N/A] You strongly encouraged to include justification to your answer, either by referencing the appropriate section of your paper or providing a brief inline description. For example:
Did you include the license to the code and datasets? [Yes] See Section 2.
Did you include the license to the code and datasets? [No] The code and the data are proprietary.
Did you include the license to the code and datasets? [N/A]
Please do not modify the questions and only use the provided macros for your answers. Note that the Checklist section does not count towards the page limit. In your paper, please delete this instructions block and only keep the Checklist section heading above along with the questions/answers below.
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes] See Section 1
(b) Did you describe the limitations of your work? [Yes] See Section 6
(c) Did you discuss any potential negative societal impacts of your work? [No]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] The paper conforms with the provided guidelines.
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section
(b) Did you include complete proofs of all theoretical results? [Yes] See the supplemental material (Appendix.pdf)
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed reproduce the main ex- perimental results (either in the supplemental material or as a URL)? [Yes]] See the supplemental material. Our codes contains README.md reproduce the our experi- mental results.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 5
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes] See Section 5
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Section 5.1
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [N/A]
(b) Did you mention the license of the assets? [Ys]] See the supplemental materials. Our codes contains LICENSE.txt.
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] See the supplemental materials containing our souree codes
(d) Did you discuss whether and how consent was obtained from people whose data you'l 're using/curating? [N/A] Our dataset was reproduced by ourselves based on the literature (c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [Yes]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
12


Reviews:
[{"review_id": "kxbQxMFQ5Lq", "reviewer": "Reviewer_shMz", "summary": "This is a system identification method for nonlinear dynamical systems represented as ODEs based on input and output signals. Authors use neural network for such a system identification. The novelty of the paper is that the identification also comes with the guarantee on the stability. The enabling technique is to enforce/projection onto Hamiltonian-Jacobi inequality manifold for L-2 stability. ", "questions": "Can you clarify on the application of this method to chaotic ODEs? Is it a limitation? What challenges do you foresee? \n\nIs it fair to say, this work is an extension of Neural ODE to enforce stability? If not, can authors briefly comment on such a possibility? Seems like this method, since it uses Euler method, can be faster than Neural ODE. What about generalization? Also what are the impact of $\\Delta T$ on the convergence and numerical stability of learning. \n\nIt's a bit confusing that \\hat{y} is generated by modified dynamics (f,G,h); while the nominal dynamics is denoted by (\\hat{f},\\hat{G},\\hat{h}). Is there a reason authors used such notation?\n\nMaybe a schematic of the method or an algorithm can help to better clarify the method. It is still unclear how do authors in general prefer to represent V. \n\nIs it possible that the original system itself is not stable and then this system stabilize the neural network representation of that?", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\nSystem identification with neural networks has the potential to better generalize than existing methods. However, there is no guarantee that such methods are stable. This paper brings ideas from control and dynamical systems community and apply them to the learning-based system identification to enforce L2 stability. The resulting networks can be used for prediction of dynamical systems in the long-term. \n\nWeaknesses:\nThe method can be applied to nonlinear dynamical systems with fixed points or limit cycles. It is not clear if the method is applicable to more generic dynamical systems with strange attractors.\nThe method described is still a bit unclear. \nAlthough authors compare with several DNN architectures, more comparison should be made with traditional system identification methods rather than deep-learning based methods (e.g. ERA, SINDY, DMD, etc.). ", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "9Yq45UBub8c", "reviewer": "Reviewer_p8Kb", "summary": "Noticing the input-output stability in dynamical systems, the paper proposed a projection method to enforce the Hamilton-Jacobi inequality in learning such dynamics with neural network. The projection problem was converted into a quadratic constraint quadratic programming problem and solved analytically. The paper also provided experimental results to show the effectiveness of such methods.", "questions": "1. I’m not convinced by the introductory examples in figure 1. The author claims the reason is “magnitude of this input signal was larger than that of the input signals in the training dataset”. However, it looks like the output signal is unbounded (thus not satisfying I/O stability under finite step activation) regardless of input signal magnitude.\n2. As for notation, the author denotes $(\\hat{f},\\hat{G},\\hat{h})$ as nominal dynamics and $(f,G,h)$ as the modified dynamics in 3.1. It is confusing since $(f,G,h)$ denotes the true dynamics in equation 1. I suggest using other styles.\n3. In the loss function, why the 2nd term is needed here? Since $\\hat{y}$ is calculated by modified $(f,G,h)$, then such trajectory must satisfy the Hamilton Jacobian inequality. nominal $(\\hat{f},\\hat{G},\\hat{h})$ is not directly used in prediction, why do we need to penalize nominal $(\\hat{f},\\hat{G},\\hat{h})$?\n4. Regarding the computational time, I’m curious what is bottleneck for such long hours (60 GPU day in total). The toy example is simply single variable ODE, the Glucose-insulin examples only have 1000 input and output signals and seems to be low dimensional as well.", "limitations": "There is no potential negative societal impact for this work.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n1. The paper structure is easy to follow and the hierarchy of math derivation is clear. The experiments are well organized and helpful for theorem illustration.\n2. Despite many recent papers to enforce certain properties in learning dynamics, the input-output stability idea is novel.\n\n\nWeakness\n\n1. The problem formulation (eqn 1) with $ \\dot{x} $ linearly dependent on u limits the application to certain systems. But this is common treatment in dynamical system literature and is not a big concern.\n\n2. The projection idea can be seen as a continuum of [A],[B], with modification on projection methods to enforce Lyapunov like properties. The setups in the above two papers are more flexible, but lack the input u and h() output function. The projection method in this paper essentially downgrades to the setup in [A] with y=x and u=0. I’m not confident but the similar model performance might also be achieved from previous works. Please correct me if I'm wrong or provide comparison experiments.\n\n[A]: Gaurav Manek and J. Zico Kolter. Learning stable deep dynamics models. In Advances in Neural 298 Information Processing Systems (NeurIPS), volume 32, 2019.\n[B]Nathan Lawrence, Philip Loewen, Michael Forbes, Johan Backstrom, and Bhushan Gopaluni. Almost 300 surely stable deep dynamics. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, 301 2020.\n\n3. In the main experimental results(figure 1 and 6), I’m concerned if the bad performance of the step input response is mainly due to model never been unexposed to such input during the training process. The results' improvement from the proposed method could be trivial if adding similar examples into training data could address such issue.\n\n4. the experiments could be stronger to illustrate the application of the method on more complex and diverse examples. The toy example is too simple, and the general audience are not familiar with the latter glucose-insulin system in terms of complexity, modeling difficulty, system stability etc. The current examples are not representative. This is my major concern.\n\n5. More experimental details are needed, such as neural network structure, total parameter number in both NN and other systemID methods, training hyperparameters etc.", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "0ttrKWU2Qo", "reviewer": "Reviewer_adB5", "summary": "They propose an approach that uses the Hamiltonian-Jacobi (HJ) inequality to guarantee that the learned nonlinear dynamics have input-output stability.\n\nThis problem is fundamental and important in the sense that it bridges the gap between theoretical control theory and data-driven control for unknown nonlinear systems using well-established tools from control theory like the HJ inequality.\n", "questions": "1) What are the assumptions for (f,G,h) after Eqn(1)? The HJ inequality does not guarantee input-output stability if (f, G,h) are not continuous or differentiable. Please recheck.\n\n2) Please cite the relevant reference for Proposition 1\n\n3) As I am aware of the HJ inequality, apart from V being a positive definite function (which is mentioned), V must also be continuously differentiable and the value of V(0) must be zero which I think is not mentioned in your manuscript.\n\n4) Guaranteeing that the HJ inequality is satisfied over dataset D (by minimizing the loss function given in (9)) does not ensure that the HJ inequality is satisfied for all points in the feasible domain of (x0,u,y). This is a major challenge that I observed in many recent relevant papers. However, some papers effectively address this issue by using SMT solvers where the SMT solver generates counter-examples that do not satisfy the inequality and are therefore fed back to the neural network to train. Please see relevant refs titled “Neural Lyapunov Control”, “Neural Koopman Lyapunov Control”, and “Neural Lyapunov Redesign” and the relevant references therein.\n\n5) In eqn after (9), what does R(HJ(..)) denote? It is not clear to me from the write-up.\n\n6) Please give some insight on whether one can compute a valid V that satisfies the HJ inequality. I think that is not mentioned in the manuscript except in the part where you write “The function V can be designed as a trainable neural network [14] or as a fixed function that achieves the desired properties of the nonlinear dynamics [15]”. However, this is not very clear to me.\n\n7) Please give a brief remark on how you compute yhat from the nominal dynamics (f, G,h) given (x0,u,y) in the simulations.\nI am wondering why is “Relevant work” in Section 4 and not after Section 1?\n", "limitations": "The main limitation of the approach is point 4 under “Questions” and it is also not clear to me on know do they compute a valid V.\n", "rating": 7, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The merits of the paper are as follows:\n1) In contrast to the most recent approaches which simultaneously design a controller and learn the dynamics, this paper goes a step ahead to also ensure input-output stability for the learned dynamics using the well-established HJ inequality.\n\n2) The results are extensive in the sense that they clearly validate their approach and compare it with other state-of-the-art and recent approaches.\n\nThe weakness of the paper is as follows:\n1) Ensuring that the HJ inequality is satisfied over a finite dataset D does not guarantee the same is true over the entire feasible time\nAssumptions of (f, G, h) are not clear\n\nMore questions are in the next section “Questions”\n", "was_revised": false, "final_justification": "", "rebuttal": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept