Learning Deep Input-Output Stable Dynamics

Supposing $f_{\bf n}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, $G_{\bf n}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times m}$, and $h_{\bf n}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{l}$, a triplet map $(f_{\bf n},G_{\bf n},h_{\bf n})$ denote as nominal dynamics. Introducing a distance on the triplet maps, the nearest triplet satisfying the Hamilton-Jacobi inequality from the nominal dynamics $(f_{\bf n},G_{\bf n},h_{\bf n})$ is called modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$. The purpose of this section is to describe the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$ associated with the nominal dynamics $(f_{\bf n},G_{\bf n},h_{\bf n})$ by analytically deriving a projection onto the space satisfying the Hamilton-Jacobi inequality.

The problem of finding the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$ is written as a quadratic constraint quadratic programming (QCQP) problem for the nominal dynamics $(f_{\bf n},G_{\bf n},h_{\bf n})$. 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, we define the ramp and the clamp functions as

and define the Hamilton-Jacobi function as

where $V$ is 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, $x=0$ can be regarded as the positive invariant set and the increasing function $V$ can be designed as $V(x)=\tfrac{1}{2}x^{2}$.

The objective function of (5a) is a new distance between the nominal dynamics $(f_{\bf n},G_{\bf n},h_{\bf n})$ and the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$. 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 $k$ represents the ratio of the distance scale of $f$ to $G$ and $h$. When $k=0$, the result of this problem (5) are consistent with the projection of the conventional method that guarantee internal stability [9]. As the constant $k$ converges $1$, the modification method of Theorem 1 approaches $G_{\bf m}$ and $h_{\bf m}$ to $G_{\bf n}$ and $h_{\bf n}$, respectively. In this case, the objective function (5a) becomes the distance between $f_{\bf m}$ to $f_{\bf n}$. Therefore, the following corollary is satisfied.

Corollary 1 derives a solution of a linear programming problem rather than QCQP problems. Replacing the Hamilton-Jacobi function $\mathrm{HJ}(f_{\bf m},G_{\bf n},h_{\bf n})$ with the time derivative of a positive definite function $\nabla V^{\mathrm{T}}f$, this corollary matches the result of the conventional study [9].

When the map $h_{\bf m}$ is fixed as $h_{\bf n}$, a similar solution as Theorem 1 is derived. Although Corollary 1 is proved by changing $k$, the modified dynamics with the fixed $h_{\bf n}$ are not derived. We reprove this modified dynamics in a similar way to Theorem 1.

We represent $(f_{\bf n},G_{\bf n},h_{\bf n})$ as neural networks and denote $(f_{\bf m},G_{\bf m},h_{\bf m})$ as the $\mathcal{L}_{2}$ stable dynamics modified by Theorem 1, Corollary 1, or 2. Note that the modification depends on the candidate of $\mathcal{L}_{2}$ gain $\gamma$. Figure 2 (A) shows the sketch of this modification, where the blue region satisfies the Hamilton-Jacobi inequality.

Since the nonlinear system of the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$ is represented as ordinary differential equations (ODEs) consisting of the differentiable functions $f_{\bf n}$,$G_{\bf n}$, and $h_{\bf n}$, the techniques of training neural ODEs can be applied [14]. Once a loss function is designed, the parameters of neural networks in the modified system can be learned from given data.

where $\lambda$ and $\alpha$ are non-negative coefficients.

The first term shows the prediction error of the output signal $y$ ( Figure 2 (A)). A dataset $\cal{D}$ consists of tuples $(x_{0},u,y)$ where the initial value $x_{0}$, the input signal $u$, and the output signal $y$. The predicted output $\hat{y}$ is calculated from $x_{0}$, $u$, and the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$.

The second term aims to improve the nominal dynamics $(f_{\bf n},G_{\bf n},h_{\bf n})$ closer to the modified dynamics $(f_{\bf m},G_{\bf m},h_{\bf m})$ and is defined as

where $\varepsilon$ is a positive constant ( Figure 2 (B)). Since this term is a form of the hinge loss, $\varepsilon$ represents the magnitude of the penalties for the Hamilton-Jacobi inequality. To evaluate this inequality for any $x\in D$, we introduce a distribution of $x$ over the domain $D$. In our experiments, this distribution is decided as a Gaussian distribution ${\cal N}(\mu,\sigma^{2})$ where the mean $\mu$ is placed at the asymptotically stable point and the variance $\sigma^{2}$ is an experimental parameter. Without this second term of the loss function, there are 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 $\gamma$ can be manually designed for the application or automatically trained from data by introducing the last term. This training explores smaller $\gamma$ while keeping the same level of the prediction error (Figure 2 (C)).

Before describing the results of the experiments, this section explains the evaluation metrics and our experimental setting.

For evaluation metrics, we define the root mean square error (RMSE) and the average $\mathcal{L}_{2}$ gain (GainIO) for the given input and output signals as follows:

where $N$ is the number of signals in the dataset. $u_{i}(\cdot)$ and $y_{i}(\cdot)$ are the input and output signal at the $i$-th index, respectively. The prediction signal $\hat{y}_{i}(\cdot)$ is computed from $u_{i}(\cdot)$, the trained dynamics, and the initial state. Note that the integral contained in the $\mathcal{L}_{2}$ 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 a metric of the property of the $\mathcal{L}_{2}$ stability. Whether the target system satisfies or does not satisfy the $\mathcal{L}_{2}$ stability, the GainIO with a given finite-size dataset can be calculated. The GainIO error is defined by the absolute error between the GainIO of the test dataset and that of the prediction.

In our experiments, 90$\%$ of the dataset is 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 $\Delta t$ for the output $y$ is set as constant and the Euler method is used to solve ODEs. $x_{0}$ 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 dynamics, the clipping operation is used so that the absolute values of the states are less than ten.

For comparative methods, we use vanilla neural networks, ARX, ORT [18], MOESP [17], and piece-wise ARX[21]. In the method of vanilla neural networks, the maps $(f,G,h)$ in the nonlinear system (1) is represented by using three neural networks, i.e., this method is consistent with a method used in Figure 1. To determine the hyperparameters of comparative methods except for neural networks, the grid search is used. Note that these comparative methods only consider the prediction errors.

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.

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 is defined as

The internal system of this model has two asymptotically stable equilibrium points $x\equiv 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 at random. The input and output signals on the period $[0,10]$ are sampled with an interval $\Delta t=0.1$. In this benchmark, we set the number of dimensions of the internal system as one and use a fixed function $V(x)=\min((x-1)^{2},(x+1)^{2})$, a mixture of the two positive definite functions.

In the result of our experiments, we name the proposed methods modified by Theorem 1, Corollary 1, and 2 as DIOS-fgh, DIOS-f, and DIOS-fg, respectively. In these methods, the parameters of our loss function are set as $\lambda=0$ and $\alpha=0.01$. Also, DIOS-fgh+ uses $\lambda=0.01$ and $\alpha=0.01$ under the same conditions as DIOS-fgh. For this example (1000 input signal), it took about 1 hour using 1 GPU to learn one model training.

The results of the RMSE and the GainIO error in this experiment are shown in Figure 3. Figures 3 (A) and (B) demonstrate that a piece-wise ARX model (PWARX) gives very low RMSE but its GainIO 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 dataset size was varied to 100, 1000, and 3000, the result of the larger dataset provided the smaller RMSE for all methods. The vanilla neural networks do not consider the $\mathcal{L}_{2}$ gain, so the GainIO error was large in the case of the low RMSE.

Figure 4 shows the relationship between $x$ and $f(x)$ in the trained system in (1) to compare the vanilla neural networks, DIOS-fgh and DIOS-fgh+. DIOS-fgh and DIOS-fgh+ successfully find two stable points, i.e., the trained function $f(x)$ had roots of two stable points $x=\pm 1$ and an unstable point $x=0$. 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.

This section addresses an example of the identification of biological systems. We consider the task of learning the glucose-insulin system using a simulator [32] to construct responses for various inputs and evaluate the robustness of the proposed method for unexpected inputs. This simulator outputs the concentrations of plasma glucose $y_{1}$ and insulin $y_{2}$ 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 $\Delta t=1$ and 1000 steps for each sequence. In this benchmark, we set the number of dimensions of the internal system as six, fix a positive definite function $V(x)=x^{2}$, and use our loss function with $\lambda=0.001$ and $\alpha=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 are 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 5 shows input and output signals in the test dataset and the predicted output by vanilla neural networks and our method (DIOS-fgh+). The $\mathcal{L}_{2}$ stability of the system using the vanilla neural networks is not guaranteed. Since the proposed method guarantees the $\mathcal{L}_{2}$ stability, the output signals of DIOS-fgh+ are bounded even if the input signals are unexpectedly large.

To demonstrate this, we conducted an additional experiment using the trained system. Figure 6 shows the transition of output behavior caused by the magnitude of the input signal changing from 2 to 10. Note that the maximum magnitude in the training dataset is one. In this experiment, $\Delta t$ 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 is always bounded. Therefore, we actually confirmed that our proposed method satisfies the $\mathcal{L}_{2}$ stability.