Neural Contraction Metrics for Robust Estimation and Control: A Convex Optimization Approach

Consider the following perturbed nonlinear system:

where $t\in\mathbb{R}_{\geq 0}$, $x:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$, $f:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$, and $d:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ with $\overline{d}=\sup_{t\geq 0}\|d(t)\|<+\infty$.

An LSTM-RNN is a neural network designed for processing sequential data with inputs $\{x_{i}\}_{i=0}^{N}$ and outputs $\{y_{i}\}_{i=0}^{N}$, and defined as $y_{i}=W_{hy}h_{i}+b_{y}$ where $h_{i}=\phi(x_{i},h_{i-1},c_{i-1})$. The activation function $\phi$ is given by the following relations: $h_{i}=o_{i}\tanh(c_{i})$, $o_{i}=\sigma(W_{xo}x_{i}+W_{ho}h_{i-1}+W_{co}c_{i}+b_{o})$, $c_{i}=f_{i}c_{i-1}+\iota_{i}\tanh(W_{xc}x_{i}+W_{hc}h_{i-1}+b_{c})$, $f_{i}=\sigma(W_{xf}x_{i}+W_{hf}h_{i-1}+W_{cf}c_{i-1}+b_{f})$, and $\iota_{i}=\sigma(W_{xi}x_{i}+W_{hi}h_{i-1}+W_{ci}c_{i-1}+b_{i})$, where $\sigma$ is the logistic sigmoid function and $W$ and $b$ terms represent weight matrices and bias vectors to be optimized, respectively. The deep LSTM-RNN can be constructed by stacking multiple of these layers [1, 2].

Since contraction analysis for discrete-time systems leads to similar results [5, 11], we define the inputs $x_{i}$ as discretized states $\{x_{i}=x(t_{i})\}_{t=0}^{N}$, and the outputs $y_{i}$ as non-zero components of the unique Cholesky decomposition of the optimal contraction metric, as will be discussed in Sec. III.

We derive one approach to sample contraction metrics by using the ConVex optimization-based Steady-state Tracking Error Minimization (CV-STEM) method [10, 11], which could handle the control design of Itô stochastic nonlinear systems. In this section, we propose its simpler formulation for nonlinear systems with bounded disturbances in order to be of practical use in engineering applications.

By Theorem 1, the problem to minimize an upper bound of the steady-state Euclidean distance between the trajectory $x_{1}(t)$ of the unperturbed system and $x_{2}(t)$ of the perturbed system (1) ((4) as $t\to\infty$) can be formulated as follows:

where $W(x,t)=M(x,t)^{-1}$ is used as a decision variable instead of $M(x,t)$. We assume that the contraction rate $\alpha$ and disturbance bound $\overline{d}$ are given in (5) (see Remark 1 on how to select $\alpha$). We need the following lemma to convexify this nonlinear optimization problem.

We are now ready to state and prove our main result on the convex optimization-based sampling.

Instead of directly using sequential data of optimal contraction metrics $\{M(x(t_{i}),t_{i})\}^{N}_{i=0}$ for neural network training, the positive definiteness of $M(x,t)$ is utilized to reduce the dimension of the target output $\{y_{i}\}_{i=0}^{N}$ defined in Sec. II.

As a result of Lemma 2, we can select $\Theta(x,t)$ defined in Theorem 1 as the unique Cholesky decomposition of $M(x,t)$ and train the deep LSTM-RNN using only the non-zero entries of the unique upper triangular matrices $\{\Theta(x(t_{i}),t_{i})\}^{N}_{i=0}$. We denote these nonzero entries as $\theta(x,t)\in\mathbb{R}^{\frac{1}{2}n(n+1)}$. As a result, the dimension of the target data $\theta(x,t)$ is reduced by $n(n-1)/2$ without losing any information on $M(x,t)$.

The pseudocode to obtain an NCM depicted in Fig. 1 is presented in Algorithm 1. The deep LSTM-RNN in Sec. II is trained with the sequential state data $\{x(t_{i})\}_{i=0}^{N}$ and the target data $\{\theta(x(t_{i}),t_{i})\}^{N}_{i=0}$ using Stochastic Gradient Descent (SGD). We note these pairs will be sampled for multiple trajectories to increase sample size and to avoid overfitting.

We apply an NCM to the state estimation problem for the following nonlinear system with bounded disturbances:

where $d_{1}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{k_{1}}$, $B:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{n\times k_{1}}$, $y:\mathbb{R}_{\geq 0}\to\mathbb{R}^{m}$, $d_{2}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{k_{2}}$, $h:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{m}$, and $G:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{m\times k_{2}}$ with $\overline{d}_{1}=\sup_{t}\|d_{1}(t)\|<+\infty$ and $\overline{d}_{2}=\sup_{t}\|d_{2}(t)\|<+\infty$. Let $W=M(\hat{x},t)^{-1}\succ 0$, $A(x,t)=(\partial f/\partial x)$, and $C(x,t)=(\partial h/\partial x)$. We design an estimator as

where $\alpha>0$. The virtual system of (9) and (10) is given as

where $d_{e}(q,t)$ is defined as $d_{e}(x,t)=B(x,t)d_{1}(t)$ and $d_{e}(\hat{x},t)=M(\hat{x},t)C(\hat{x},t)^{T}G(x,t)d_{2}(t)$. Note that (12) has $q=x$ and $q=\hat{x}$ as its particular solutions. The differential dynamics of (12) with $d_{e}=0$ is given as

We have the following lemma for deriving a condition to guarantee the local contraction of (12) in Theorem 3.

The following theorem along with this lemma guarantees the exponential stability of the estimator (10).

We have the following proposition to sample optimal contraction metrics for the NCM-based state estimation.

We have an analogous result for state feedback control.

Algorithm 1 along with Proposition 1 and Corollary 1 returns NCMs to compute $\hat{x}(t)$ of (10) and $u(t)$ of (17) for state estimation and control in real-time. They also provide the bounded error tube (see Theorem 1, [15, 16]) for robust motion planning problems as will be seen in Sec. V.

The similarity of Corollary 1 to Proposition 1 stems from the estimation and control duality due to the differential nature of contraction analysis as is evident from (13) and (20). Analogously to the discussion of the Kalman filter and LQR duality in LTV systems, this leads to two different interpretations on the weight of $\nu$ ($\overline{d}_{2}\bar{c}\bar{g}/\gamma$ in (16) and $\lambda$ in (19)). As discussed in Remark 1, one way is to see it as a penalty on the induced 2-norm of feedback gains. Since $\overline{d}_{2}=0$ in (16) means no noise acts on $y(t)$, it can also be viewed as an indicator of how much we trust the measurement $y(t)$: the larger the weight of $\nu$, the less confident we are in $y(t)$. These agree with our intuition as smaller feedback gains are suitable for systems with larger measurement uncertainty.

We first consider state estimation of the Lorentz oscillator with bounded disturbances described as $\dot{x}=f(x)+d_{1}(t)$ and $y=Cx+d_{2}(t)$, where $f(x)=[\sigma(x_{2}-x_{1}),x_{1}(\rho-x_{3})-x_{2},x_{1}x_{2}-\beta x_{3}]^{T}$, $x=[x_{1},x_{2},x_{3}]^{T}$, $\sigma=10$, $\beta=8/3$, $\rho=28$, $C=[1~{}0~{}0]$, $\sup_{t}\|d_{1}(t)\|=\sqrt{3}$, and $\sup_{t}\|d_{2}(t)\|=1$. We use $dt=0.1$ for integration, with one measurement $y$ per $dt$.

We consider an optimal motion planning problem of the planar spacecraft dynamical system, given as $\dot{x}=Ax+B(x)u+d(t)$, where $u\in\mathbb{R}^{8}$, $\sup_{t}\|d(t)\|=0.15$, and $x=[p_{x},p_{y},\phi,\dot{p}_{x},\dot{p}_{y},\dot{\phi}]^{T}$ with $p_{x}$, $p_{y}$, and $\phi$ being the horizontal coordinate, vertical coordinate, and yaw angle of the spacecraft, respectively. The constant matrix $A$ and the state-dependent actuation matrix $B(x)$ are defined in [30]. All the parameters of the spacecraft are normalized to $1$.