Robustness of Numerically Computed Contraction Metrics

Abstract

Contraction metrics are an important tool to show the existence of exponentially stable equilibria or periodic orbits, and to determine a subset of their basin of attraction. One of the main advantages is that contraction metrics are robust with respect to perturbations of the system, i.e. a contraction metric for one particular system is also a contraction metric for a perturbed system. In this paper, we discuss numerical methods to compute contraction metrics for dynamical systems, with exponentially stable equilibria or periodic orbits, and perform perturbation analysis. In particular, we prove the robustness of such metrics to perturbations of the system and give concrete bounds. The results imply that a contraction metric, computed for a particular system, remains a contraction metric for the perturbed system. We illustrate our results by computing contraction metrics for systems from the literature, both with exponentially stable equilibria and exponentially stable periodic orbits, and then investigate the validity of the metrics for perturbed systems. Parts of the results are published in Giesl et al. (Proceedings of the 18th International Conference on Informatics in Control, Automation and Robotics (ICINCO), 2021).

Appendix A Description of the Algorithm

In this section, we will provide a detailed version of the algorithm both for the periodic orbit case, and the equilibrium point case, presented next to each other for the convenience of the reader. Given is a system ${{\dot{\varvec{x}}}}=\varvec{f}(\varvec{x})$, with $\varvec{f}\in C^{\sigma + 1}(\mathbb {R}^n;\mathbb {R}^n)$, where $\sigma =\lceil k+\frac{n+1}{2} \rceil$ and $k\ge 2$ if n is odd and $k\ge 3$ if n is even, so that the minimum smoothness needed for the contraction metric M and its optimal recovery S is guaranteed.

The idea is to increase the number of collocation points and vertices gradually so that we obtain a small enough fill distance and fine enough triangulation. The steps of the algorithm are as follows:

• STEP 0. setting the constants and parameters.

  Fix $d \ge 2\sqrt{n}$, $k\ge 2$ if n is odd and $k\ge 3$ if n is even, $c>0$, and the Wendland function $\psi _0(r):=\psi _{l,k}(cr)$ with $l=\lfloor \frac{n}{2}\rfloor +k+1$. Denote $\psi _{q+1}(r)=\frac{1}{r}\frac{d\psi _q}{dr}(r)$ for $q=0,1$.

  Further, fix a compact set ${\mathcal C}\subset \mathbb {R}^n$ that we triangulate and where we want to compute a contraction metric for the system, and an open set ${\mathcal O}\supset {\mathcal C}$. Start with $h_\text {collo}>0$, $h_\text {triang}>0$.

• STEP I. RBF subroutine

  Choose a set of pairwise distinct points $X=\{x_1,\ldots ,x_N\}$ in ${\mathcal O}$ as collocation points with fill distance $h_{X,{\mathcal O}} \le h_\text {collo}$. To obtain a solution of the optimal recovery problem based on RBF approximation we follow these steps:

  (a):

  Compute the coefficients $b_{(\ell ,i,j),(k,\mu ,\nu )}$ defined as

  $$\begin{aligned}&b_{(\ell ,i,j),(k,\mu ,\nu )}\nonumber \\&\quad = \psi _0(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\bigg [\sum _{p=1}^n Df_{pi}(\varvec{x}_\ell )Df_{p\mu }(\varvec{x}_k)\delta _{\nu j}+Df_{\mu i}(\varvec{x}_\ell ) Df_{j\nu }(\varvec{x}_k)\nonumber \\&\qquad +Df_{i\mu }(\varvec{x}_k) Df_{\nu j}(\varvec{x}_\ell )+\delta _{i\mu } \sum _{p=1}^n Df_{p\nu }(\varvec{x}_k)Df_{pj}(\varvec{x}_\ell )\bigg ]\nonumber \\&\qquad +\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_k-\varvec{x}_\ell ,\varvec{f}(\varvec{x}_k)\rangle \left[ Df_{\mu i}(\varvec{x}_\ell )\delta _{\nu j}+\delta _{i\mu } Df_{\nu j}(\varvec{x}_\ell )\right] \nonumber \\&\qquad +\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_\ell -\varvec{x}_k,\varvec{f}(\varvec{x}_\ell )\rangle \left[ Df_{i\mu }(\varvec{x}_k)\delta _{\nu j}+\delta _{i\mu }Df_{j\nu }(\varvec{x}_k)\right] \nonumber \\&\qquad -\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{f}(\varvec{x}_\ell ),\varvec{f}(\varvec{x}_k)\rangle \delta _{i\mu }\delta _{j\nu }\nonumber \\&\qquad +\psi _2(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_k-\varvec{x}_\ell ,\varvec{f}(\varvec{x}_k)\rangle \langle \varvec{x}_\ell -\varvec{x}_k,\varvec{f}(\varvec{x}_\ell )\rangle \delta _{i\mu }\delta _{j\nu }. \end{aligned}$$

  (A1)

  for $1 \le k, \ell \le N$, and $1 \le i, j, \mu , \nu \le n$. Here and below, $\langle \varvec{x},\varvec{y}\rangle =\sum _{k=1}^nx_ky_k$ denotes the Euclidean scalar product.

  ($a'$):

  For the periodic orbit case we need to fix a $c_0 \in \mathbb {R}^+$, and $\varvec{y}_0 \in {\mathcal A}(\Omega )$, then calculate the coefficients $b_{(\ell ,i,j),(k,\mu ,\nu )}$, $b_{0,(k,\mu ,\nu )}$, $b_{(\ell ,i,j),0}$ and $b_{0,0}$ for $1\le k,\ell \le N$, $1\le i,j,\mu ,\nu \le n$ using the following formulas

  $$\begin{aligned}&b_{0,(k,\mu ,\nu )}\nonumber \\&\quad = \psi _0(\Vert \varvec{x}_k-\varvec{y}_0\Vert _2)\bigg [\sum _{p=1}^n V_{p\mu }(\varvec{x}_k)f_p(\varvec{y}_0)f_\nu (\varvec{y}_0)+\sum _{p=1}^n V_{p\nu }(\varvec{x}_k)f_p(\varvec{y}_0)f_\mu (\varvec{y}_0)\bigg ]\nonumber \\&\qquad +\psi _1(\Vert \varvec{x}_k-\varvec{y}_0\Vert _2)\langle \varvec{x}_k-\varvec{y}_0,\varvec{f}(\varvec{x}_k)\rangle f_\mu (\varvec{y}_0)f_\nu (\varvec{y}_0) \end{aligned}$$

  (A2)

  $$\begin{aligned}&b_{0,0}\nonumber \\&\quad =\psi _0(0)\Vert \varvec{f}(\varvec{y}_0)\Vert _2^4. \end{aligned}$$

  (A3)

  $$\begin{aligned}&b_{(\ell ,i,j),(k,\mu ,\nu )}\nonumber \\&\quad = \psi _0(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\bigg [\sum _{p=1}^n V_{pi}(\varvec{x}_\ell )V_{p\mu }(\varvec{x}_k)\delta _{\nu j}+V_{\mu i}(\varvec{x}_\ell ) V_{j\nu }(\varvec{x}_k)\nonumber \\&\qquad +V_{i\mu }(\varvec{x}_k) V_{\nu j}(\varvec{x}_\ell )+\delta _{i\mu } \sum _{p=1}^n V_{p\nu }(\varvec{x}_k)V_{pj}(\varvec{x}_\ell )\bigg ]\nonumber \\&\qquad +\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_k-\varvec{x}_\ell ,\varvec{f}(\varvec{x}_k)\rangle \left[ V_{\mu i}(\varvec{x}_\ell )\delta _{\nu j}+\delta _{i\mu } V_{\nu j}(\varvec{x}_\ell )\right] \nonumber \\&\qquad +\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_\ell -\varvec{x}_k,\varvec{f}(\varvec{x}_\ell )\rangle \left[ V_{i\mu }(\varvec{x}_k)\delta _{\nu j}+\delta _{i\mu }V_{j\nu }(\varvec{x}_k)\right] \nonumber \\&\qquad -\psi _1(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{f}(\varvec{x}_\ell ),\varvec{f}(\varvec{x}_k)\rangle \delta _{i\mu }\delta _{j\nu }\nonumber \\&\qquad +\psi _2(\Vert \varvec{x}_k-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_k-\varvec{x}_\ell ,\varvec{f}(\varvec{x}_k)\rangle \langle \varvec{x}_\ell -\varvec{x}_k,\varvec{f}(\varvec{x}_\ell )\rangle \delta _{i\mu }\delta _{j\nu } \end{aligned}$$

  (A4)

  $$\begin{aligned}&\text {and }b_{(\ell ,i,j),0}\nonumber \\&\quad = \psi _0(\Vert \varvec{y}_0-\varvec{x}_\ell \Vert _2)\bigg [\sum _{p=1}^n V_{pi}(\varvec{x}_\ell )f_p(\varvec{y}_0)f_j(\varvec{y}_0)+\sum _{p=1}^n V_{pj}(\varvec{x}_\ell )f_p(\varvec{y}_0)f_i(\varvec{y}_0)]\nonumber \\&\qquad +\psi _1(\Vert \varvec{y}_0-\varvec{x}_\ell \Vert _2)\langle \varvec{x}_\ell -\varvec{y}_0,\varvec{f}(\varvec{x}_\ell )\rangle f_i(\varvec{y}_0)f_j(\varvec{y}_0). \end{aligned}$$

  (A5)

  (b):

  Calculate the coefficients $c_{(\ell ,i,j),(k,\mu ,\nu )}$ defined as

  $$\begin{aligned} c_{(\ell ,i,j),(k,\mu ,\nu )} = \frac{1}{4}\left( b_{(\ell ,i,j),(k,\mu ,\nu )}+b_{(\ell ,j,i),(k,\nu ,\mu )}+b_{(\ell ,i,j),(k,\nu ,\mu )} +b_{(\ell ,j,i),(k,\mu ,\nu )} \right) , \end{aligned}$$

  (A6)

  where we assume $\mu \le \nu$ and $i\le j$. The coefficients $c_{\cdot ,\cdot }$ form a symmetric matrix of size $N\frac{n(n+1)}{2}$.

  ($b'$):

  For the periodic orbit case, in addition to the coefficients formula in (b), we need one more row and column in the matrix so for all $(\ell ,i,j)$ with $1\le \ell \le N$, $1\le i\le j\le n$ we have

  $$\begin{aligned} c_{0,0}\,=\, & {} b_{0,0}, \end{aligned}$$

  (A7)

  $$\begin{aligned} c_{(\ell ,i,j),0}\,=\, & {} b_{(\ell ,i,j),0}, \end{aligned}$$

  (A8)

  $$\begin{aligned} c_{0,(k,\mu ,\nu )}= & {} \frac{1}{2}\left( b_{0,(k,\mu ,\nu )}+b_{0,(k,\nu ,\mu )}\right) \ = \ b_{0,(k,\mu ,\nu )}, \end{aligned}$$

  (A9)

  where we assume $\mu \le \nu$. The coefficients $c_{\cdot ,\cdot }$ form a symmetric matrix of size $N\frac{n(n+1)}{2}+1$.

  (c):

  Determine $\gamma _k^{(\mu ,\nu )}$, by solving the linear system

  $$\begin{aligned}&\sum _{k=1}^N\sum _{1\le \mu \le \nu \le n} c_{(\ell ,i,j),(k,\mu ,\nu )}\gamma _k^{(\mu ,\nu )}\nonumber \\&\quad =(F (S)(\varvec{x}_\ell ))_{i,j} =\lambda _\ell ^{(i,j)}(S) = -C_{ij} \end{aligned}$$

  (A10)

  for $1 \le \ell \le N$, and $1 \le i \le j \le n$. Note that (A10) is a system of $N\frac{n(n+1)}{2}$ equations in $N\frac{n(n+1)}{2}$ unknowns and we usually choose $C_{ij}$ as constants, although the method also works for $C_{ij}(\varvec{x}_\ell )$.

  ($c'$):

  For the periodic orbit case the $\gamma _k^{(\mu ,\nu )}$ solve the modified linear system

  $$\begin{aligned} \sum _{k=1}^N\sum _{1\le \mu \le \nu \le n} c_{(\ell ,i,j),(k,\mu ,\nu )}\gamma _k^{(\mu ,\nu )}+ c_{(\ell ,i,j),0}\gamma _0= & {} -C_{ij}(\varvec{x}_\ell ) \end{aligned}$$

  (A11)

  $$\begin{aligned} \sum _{k=1}^N\sum _{1\le \mu \le \nu \le n} c_{0,(k,\mu ,\nu )}\gamma _k^{(\mu ,\nu )}+ c_{0,0}\gamma _0= & {} c_0 \Vert \varvec{f}(\varvec{y}_0)\Vert _2^4 \end{aligned}$$

  (A12)

  for $1\le \ell \le N$, $1\le i\le j\le n$. Note that this is a system of $N\frac{n(n+1)}{2}+1$ equations in $N\frac{n(n+1)}{2}+1$ unknowns.

  (d):

  Compute $\beta _k\in \mathbb {S}^{n\times n}$ from $\gamma _k$; recalling that

  $$\begin{aligned} \beta _k^{(j,i)}= & {} \beta _k^{(i,j)}=\frac{1}{2}\gamma _k^{(i,j)} \text{ if } i\not = j, \\ \beta _k^{(i,i)}= & {} \gamma _k^{(i,i)}. \end{aligned}$$

  ($d'$):

  For the periodic orbit case, we need relations in (d) and $\beta _0=\gamma _0$.

  (e):

  We now have a formula for the optimal recovery

  $$\begin{aligned} S(\varvec{x})= & {} \sum _{k=1}^N \bigg [ \psi _0(\Vert \varvec{x}_k-\varvec{x}\Vert _2)\left[ D\varvec{f}(\varvec{x}_k)\beta _k+ \beta _k D\varvec{f}(\varvec{x}_k)^T \right] \nonumber \\&~~~~~ + \, \psi _1(\Vert \varvec{x}_k-\varvec{x}\Vert _2)\langle \varvec{x}_k-\varvec{x}, \varvec{f}(\varvec{x}_k)\rangle \beta _k\bigg ]. \end{aligned}$$

  (A 13)

  ($e'$):

  For the periodic orbit case, we can compute $S(\varvec{x})$ with

  $$\begin{aligned} S(\varvec{x})= & {} \sum _{k=1}^N\bigg [\psi _0(\Vert \varvec{x}_k-\varvec{x}\Vert _2)\left[ V(\varvec{x}_k)\beta _k+ \beta _k V (\varvec{x}_k)^T \right] \nonumber \\&+\psi _1(\Vert \varvec{x}_k-\varvec{x}\Vert _2)\langle \varvec{x}_k-\varvec{x},\varvec{f}(\varvec{x}_k)\rangle \beta _k\bigg ]\nonumber \\&+\beta _0 \psi _0(\Vert \varvec{y}_0-\varvec{x}\Vert _2)\varvec{f}(\varvec{y}_0)\varvec{f}(\varvec{y}_0)^T\,. \end{aligned}$$

  (A 14)

• STEP II. CPA subroutine

  Fix an (h, d)-bounded triangulation ${\mathcal T}$ with $h \le h_\text {triang}$ and ${\mathcal D}_{\mathcal T}={\mathcal C}$. Note that d denotes a bound on the degeneracy of the simplices, for details see [39].

  (a):

  Compute the values $S(\varvec{y})$ at the vertices of the triangulation $\varvec{y}\in {\mathcal V}_{\mathcal T}$ and check if they are positive definite. If not, decrease $h_\text {collo}$ by a factor, e.g. $h_\text {collo} \leftarrow h_\text {collo} /2$, and go back to STEP I.

  (b):

  Fix constants $B_{2,\nu }$, and $B_{3,\nu }$ as upper bounds on the second- and third-order derivatives of the components $f_k$ of $\varvec{f}$ on each simplex ${\mathfrak S}_\nu$. Then compute the error term $E^e_{\nu } := n^{2} (1+4\sqrt{n}) B_{2,\nu } \Vert \nabla P_{ij}^\nu \Vert _1 +2 \, n^{3} B_{3,\nu } P_\nu ,$ where

  $$\begin{aligned} P_\nu :=\max _{\varvec{x}\in {\mathfrak S}_\nu }\Vert P(\varvec{x})\Vert _2=\max _{k=0,1,\ldots ,n}\Vert P(\varvec{x}_k)\Vert _2. \end{aligned}$$

  ($b'$):

  For the periodic orbit case we need $B_{i,\nu }$ with $i = 0, \cdots , 3$; together with $B_{V_{1,\nu }}$, and $B_{V_{2,\nu }}$ as the upper bounds on the first- and second-order derivatives of the components $V_{lj}$ of V. Then compute

  $$\begin{aligned}&E^p_{\nu } := n^2 \cdot \\&\quad \left( (4 \sqrt{n} B_{V_{1,\nu }} + B_{2,\nu }) \Vert \nabla P_{ij}^\nu \Vert _1 + 2 n B_{V_{2,\nu }} P_\nu + 2 \,\kappa _{\nu }^*\, B_{0,\nu } B_{2,\nu } + 2 \,\kappa _{\nu }^*\,B_{1,\nu }^2 \right) . \end{aligned}$$

  (c):

  Check whether Constraints (VP2) of Verification Problem 2 (or 1) are fulfilled. If not, then reduce $h_\text {triang}$ by a factor, e.g. $h_\text {triang} \leftarrow h_\text {triang}/8$, and repeat STEP II from the beginning. If the conditions still do not hold, decrease $h_\text {collo}$ by a factor, e.g. $h_\text {collo} \leftarrow h_\text {collo} /2$, and go back to STEP I.

• STEP III. creating the contraction metric

  Build the metric $P:{\mathcal C}\rightarrow \mathbb {S}^{n\times n}$ as the CPA interpolation of the values $P(\varvec{y})$, $\varvec{y}\in {\mathcal V}_{\mathcal T}$. P is a contraction metric for the system on any compact $\widetilde{K}\subset {\mathcal C}^\circ$.

Remark 1

Note that in most applications it is more practical to use a relaxed version of the procedure above to compute a contraction metric. If the matrices $P(\varvec{y})$, $\varvec{y}\in {\mathcal V}_{\mathcal T}$, in STEP II (a) are positive definite in a reasonably large part of ${\mathcal C}$, then one can proceed to the next sub-steps. Further, if additionally Constraints (VP2) of Verification Problem 2 (or 1) are fulfilled in a reasonably large part of ${\mathcal C}$ in STEP II (c), then one can proceed to STEP III.

The CPA interpolation P will then be a contraction metric on any compact subset $\widetilde{K}$ of the interior of the area where P is both positive definite, asserted by Constraints (VP1) of Verification Problem 1, and fulfills Constraints (VP2) of Verification Problem 2 (or 1). We use this relaxed procedure in our examples.