Subgradient algorithm for computing contraction metrics for equilibria

Peter Giesl; Sigurdur Hafstein; Magnea Haraldsdottir; David Thorsteinsson; Christoph Kawan

doi:10.3934/jcd.2022030

Article Contents

2023, Volume 10, Issue 2: 281-303. Doi: 10.3934/jcd.2022030

This issue Previous Article Neural dynamic mode decomposition for end-to-end modeling of nonlinear dynamics Next Article On the preservation of second integrals by Runge-Kutta methods

Subgradient algorithm for computing contraction metrics for equilibria

1.
Department of Mathematics, University of Sussex, Falmer, BN1 9QH, UK
2.
The Science Institute, University of Iceland, Dunhagi 5,107 Reykjavik, Iceland
3.
Institute of Informatics, LMU Munich, Oettingenstraße 67, 80538 Munich, Germany

^*Corresponding author: Sigurdur Hafstein
^*Corresponding author: Sigurdur Hafstein

Received: March 2022

Revised: November 2022

Early access: January 9, 2023

Published online: January 9, 2023

The second author is supported in part by the Icelandic Research Fund under Grant 228725-051.

Abstract / Introduction Full Text(HTML) Figure(9) / Table(8) Related Papers Cited by

Abstract

We propose a subgradient algorithm for the computation of contraction metrics for systems with an exponentially stable equilibrium. We show that for sufficiently smooth systems our method is always able to compute a contraction metric on any forward-invariant compact neighbourhood of the equilibrium, which is a subset its basin of attraction. We demonstrate the applicability of our method by constructing contraction metrics for three planar and one three-dimensional systems.

Contraction metrics; stability analysis; contraction metrics; subgradient method

Keywords:

Mathematics Subject Classification: Primary: 37C75, 93D05, 37M22; Secondary: 39A30, 34D20.

Citation:

Full Text(HTML)

Figure 1. Visualization of the optimal metric, see Table 2, obtained for the time-reversed van der Pol system (5.1) on the area $ K $ using a 4th degree polynomial and 10,000 iterations. The ellipse around a point $ x_0\in {\mathbb {R}}^2 $ denotes the points of equal distance to the point with respect to the metric, i.e. $ \{x\in {\mathbb {R}}^2\mid (x-x_0) ^{ \top} M(x_0)(x-x_0) ^{ \top} = c^2\} $ with fixed $ c>0 $

Download: Full-size image PowerPoint slide

Figure 2. Visualization of the polynomial over the area $ K $ using a 4th degree polynomial and 10,000 iterations. Since $ M(x_0) = \mathrm{e}^{p(x_0)}P $, the shape of the ellipse in Figure 1 is the same in the entire area due to $ P $, and we hence show $ p_a(x,y) $ in this figure

Download: Full-size image PowerPoint slide

Figure 3. The black curve is the forward-invariant set computed in [19] for system (5.2) and the red quadrilateral is the area, on which we compute metrics for the equilibrium at $ (0,0) $. The black curve is a level set of the approximation $ v $ of the solution of $ V'(x) = -\sqrt{10^{-8}+\|f(x)\|^2} $, while the blue circles denote the areas where $ v'(x)\ge 0 $. Since the black curve does not touch any blue circles, it is a forward-invariant set

Download: Full-size image PowerPoint slide

Figure 4. Graph of the 6th degree polynomial $ p_a(x,y) $ in the metric $ M(x,y) = \mathrm{e}^{p_a(x,y)}P $ computed with the algorithm for the Speed Control system (5.2) at the equilibrium at $ (0,0) $. Note that it is not a Lyapunov function, in contrast to the construction in the proof of Theorem 3.1

Download: Full-size image PowerPoint slide

Figure 5. The black curve is the forward-invariant set computed in [19] for system (5.2) and the red quadrilateral is the area, on which we compute metrics for the equilibrium at $ (-0.7887,0) $. The black curve is a level set of the approximation $ v $ of the solution of $ V'(x) = -\sqrt{10^{-8}+\|f(x)\|^2} $, while the blue circles denote the areas where $ v'(x)\ge 0 $. Since the black curve does not touch any blue circles, it is a forward-invariant set

Download: Full-size image PowerPoint slide

Figure 6. Graph of the polynomial $ p(x,y) $ in the metric $ M(x,y) = \mathrm{e}^{p_a(x,y)}P $ for system (5.2) at the equilibrium at $ (-0.7887,0) $, obtained with the algorithm and a polynomial of degree 10

Download: Full-size image PowerPoint slide

Figure 7. The black curve shows the forward-invariant set we computed with the method from [19] for system (5.4) and the red quadrilateral is the area, on which we compute metrics for the equilibrium at the origin. The black curve is a level set of the approximation $ v $ of the solution of $ V'(x) = -\sqrt{10^{-8}+\|f(x)\|^2} $, while the blue circles denote the areas where $ v'(x)\ge 0 $. Since the black curve does not touch any blue circles, it is a forward-invariant set

Download: Full-size image PowerPoint slide

Figure 8. Graph of the polynomial $ p(x,y) $ in the metric $ M(x,y) = \mathrm{e}^{p_a(x,y)}P $ for system (5.4) at the equilibrium at the origin, obtained with the algorithm and a polynomial of degree 10

Download: Full-size image PowerPoint slide

Figure 9. The red area in the figure is the forward-invariant set computed in [20] for system (5.5)

Download: Full-size image PowerPoint slide

Table 1. Results for the time-reversed van der Pol system (5.1). The degree of the polynomial in the metric is in the first column and the second column contains the first iteration where $ -\nu<0 $, i.e. exponential stability is asserted ($ \infty $ when $ -\nu \ge 0 $ for all iterations). The third column reports the lowest value for $ -\nu $ in $ 10,000 $ iterations and the last column reports the time needed for the computation of the $ 10,000 $ iterations on AMD ThreadRipper 3990X (64 cores)

Time-reversed van der Pol (5.1)
degree	first neg. itr.	lowest $ -\nu $	time
0	$ \infty $	1.447348563020335	445 s
1	$ \infty $	1.448404890308665	468 s
2	2	-0.3133495012442786	545 s
3	7	-0.2952854122048584	639 s
4	3	-0.3692555275877225	885 s
5	7	-0.357751157505084	1134 s
6	3	-0.3639232968985002	1439 s

| Show Table

DownLoad: CSV

Table 2. The contraction metric with the best bound on the exponential rate of attraction for (5.1)

$ p_a(x,y) $
term	coefficient	term	coefficient
$ x $	0.005874011679514269	$ xy^2 $	-0.023594561234174
$ y $	0.002544185767556669	$ y^3 $	0.2025391618590853
$ x^2 $	3.133966383013717	$ x^4 $	1.61103471945693
$ xy $	-0.6857166309584187	$ x^3y $	-0.2487310019195522
$ y^2 $	0.3002580326389737	$ x^2y^2 $	0.3543359078177429
$ x^3 $	-0.01456990449537407	$ xy^3 $	0.2471473588049839
$ x^2y $	-0.05483029835239338	$ y^4 $	0.9675867347507244

| Show Table

DownLoad: CSV

Table 3. Results for the Speed Control system (5.2) at equilibrium $ (0,0) $. The degree of the polynomial in the metric is in the first column and the second column contains the first iteration where $ -\nu<0 $, i.e. exponential stability is asserted ($ \infty $ when $ -\nu \ge 0 $ for all iterations). The third column reports the lowest value for $ -\nu $ in $ 10,000 $ iterations and the last column reports the time needed for the computation of the $ 10,000 $ iterations on AMD ThreadRipper 3990X (64 cores)

Speed Control system (5.2) at equilibrium (0, 0)
degree	first neg. itr.	lowest $ -\nu $	time
0	$ \infty $	0.02658980657159118	431 s
1	1628	-0.003808211430038129	462 s
2	1085	-0.004862705555080483	546 s
3	1081	-0.004874280676476273	688 s
4	1081	-0.004874531739145288	860 s
5	1081	-0.004874532732511677	1097 s
6	1081	-0.004874532751723261	1419 s

| Show Table

DownLoad: CSV

Table 4. The contraction metric with the best bound on the exponential rate of attraction for (5.2)

$ p_a(x,y) $
term	coefficient	term	coefficient
$ x $	-0.6406990025124931	$ x^5 $	-0.0001312413688617601
$ y $	0.5224323624975147	$ x^4y $	0.000103857484910813
$ x^2 $	0.1189296701550616	$ x^3y^2 $	-8.838844320265229$ \cdot10^{-5} $
$ xy $	-0.07400507177495169	$ x^2y^3 $	8.619846913687634$ \cdot10^{-5} $
$ y^2 $	0.06034845248479573	$ xy^4 $	-5.08894875827791$ \cdot10^{-5} $
$ x^3 $	-0.01228598722340182	$ y^5 $	7.926099979155098$ \cdot10^{-5} $
$ x^2y $	0.009450037931185229	$ x^6 $	1.458034136517088$ \cdot10^{-5} $
$ xy^2 $	-0.006313699379560556	$ x^5y $	-1.311184472170316$ \cdot10^{-5} $
$ y^3 $	0.008513732143946087	$ x^4y^2 $	1.072007441348303$ \cdot10^{-5} $
$ x^4 $	0.001521991939021997	$ x^3y^3 $	-8.196794919518857$ \cdot10^{-6} $
$ x^3y $	-0.001284245468051136	$ x^2y^4 $	7.786137296316773$ \cdot10^{-6} $
$ x^2y^2 $	0.001027217693957349	$ xy^5 $	-3.724746494883292$ \cdot10^{-6} $
$ xy^3 $	-0.000539000261930083	$ y^6 $	5.699738404557221$ \cdot10^{-6} $
$ y^4 $	0.0006766508585524121

| Show Table

DownLoad: CSV

Table 5. Results for the Speed Control system (5.2) at equilibrium $ (-0.7887,0) $. The degree of the polynomial in the metric is in the first column and the second column contains the first iteration where $ -\nu<0 $, i.e. exponential stability is asserted ($ \infty $ when $ -\nu \ge 0 $ for all iterations). The third column reports the lowest value for $ -\nu $ in $ 10,000 $ iterations and the last column reports the time needed for the computation of the $ 10,000 $ iterations on AMD ThreadRipper 3990X (64 cores)

Speed Control system (5.2) at equilibrium (−0.7887, 0)
degree	first neg. itr.	lowest $ -\nu $	time
0	$ \infty $	0.1222798935657396	457 s
1	173	-0.05644550658291747	481 s
2	109	-0.01663507890691926	551 s
3	885	-0.01316624664094682	697 s
4	182	-0.04806445535693586	873 s
5	90	-0.1022075700247213	1129 s
6	64	-0.1330993315433471	1413 s
7	54	-0.1486330006840473	1785 s
8	55	-0.1623406730007267	2160 s
9	54	-0.1785195291961935	2607 s
10	56	-0.189294001167907	3092 s

| Show Table

DownLoad: CSV

Table 6. Results for the Moore–Greitzer jet engine model (5.4). The degree of the polynomial in the metric is in the first column and the second column contains the first iteration where $ -\nu<0 $ and exponential stability is asserted ($ \infty $ when $ -\nu \ge 0 $ for all iterations). The third column reports the lowest value for $ -\nu $ in $ 10,000 $ iterations and the last column reports the time needed for the computation of the iterations on AMD ThreadRipper 3990X (64 cores)

Moore–Greitzer jet engine model (5.4)
degree	first neg. itr.	lowest $ -\nu $	time
0	$ \infty $	0.206199297541573	462 s
1	$ \infty $	0.2052385876784619	513 s
2	878	-0.01738829435500433	577 s
3	625	-0.01920738923339776	709 s
4	44	-0.06668089661574375	885 s
5	36	-0.06946535756779648	1148 s
6	28	-0.07271224180134611	1426 s
7	28	-0.07323485077040126	1794 s
8	28	-0.07361891063939163	2186 s
9	28	-0.07366357746042262	2621 s
10	28	-0.0737091329581808	3115 s

| Show Table

DownLoad: CSV

Table 7. Results for the three-dimensional system (5.4). The degree of the polynomial in the metric is in the first column and the second column contains the first iteration where $ -\nu<0 $ and exponential stability is asserted ($ \infty $ when $ -\nu \ge 0 $ for all iterations). The third column reports the lowest value for $ -\nu $ in $ 10,000 $ iterations and the last column reports the time needed for the computation of the iterations on Intel i9900K (8 cores, much slower than the processor used in the other computations)

Three dimensional example (5.4)
degree	first neg. itr.	lowest $ -\nu $	time
0	$ \infty $	0.8099724835857676	7,734 s
1	$ \infty $	0.8530578751350688	8,486 s
2	32	-0.1896101138448129	10,889 s
3	48	-0.1712718552180064	17,334 s
4	26	-0.221581678986785	28,712 s
5	30	-0.2160216478649938	45,265 s
6	26	-0.2229981124760225	71,993 s

| Show Table

DownLoad: CSV

Table 8. The contraction metric with the best bound on the exponential rate of attraction for (5.4) with a polynomial of degree $ 2 $

$ p_a(x,y,z) $
term	coefficient	term	coefficient
$ x $	2.886921642216002$ \cdot10^{-5} $	$ xz $	-0.01496681273321166
$ y $	3.103209165308687$ \cdot10^{-6} $	$ y^2 $	1.188121045477437
$ z $	2.350525324502533$ \cdot10^{-7} $	$ yz $	0.02688196544970737
$ x^2 $	1.153938280913646	$ z^2 $	1.070635223238271
$ xy $	-0.01809571166344076

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	N. Aghannan and P. Rouchon, An intrinsic observer for a class of Lagrangian systems, IEEE Trans. Automat. Control, 48 (2003), 936-944. doi: 10.1109/TAC.2003.812778.
[2]	W. Alt, Numerische Verfahren der Konvexen, Nichtglatten Optimierung, Teubner, 2004. doi: 10.1007/978-3-322-80083-1.
[3]	P. N. Anh and L. T. H. An, New subgradient extragradient methods for solving monotone bilevel equilibrium problems, Optimization, 68 (2019), 2097-2122. doi: 10.1080/02331934.2019.1656204.
[4]	P. N. Anh and Q. H. Ansari, Auxiliary principle technique for hierarchical equilibrium problems, J. Optimiz. Theory App., 188 (2021), 882-912. doi: 10.1007/s10957-021-01814-1.
[5]	P. N. Anh, Q. H. Ansari and H. P. Tu, DC auxiliary principle methods for solving lexicographic equilibrium problems, J. Global. Optim., 2022. doi: 10.1007/s10898-022-01200-9.
[6]	E. Aylward, P. Parrilo and J.-J. Slotine, Stability and robustness analysis of nonlinear systems via contraction metrics and SOS programming, LIDS Technical Report 2691, (2008), 1-25. arXiv: math/0603313v1.
[7]	E. M. Aylward, P. A. Parrilo and J.-J. Slotine, Stability and robustness analysis of nonlinear systems via contraction metrics and SOS programming, Automatica, 44 (2008), 2163-2170. doi: 10.1016/j.automatica.2007.12.012.
[8]	J. Björnsson and S. Hafstein, Efficient Lyapunov function computation for systems with multiple exponentially stable equilibria, Procedia Computer Science, 108 (2017), 655-664, Proceedings of the International Conference on Computational Science (ICCS), Zurich, Switzerland, 2017.
[9]	H.-D. Chiang, M. W. Hirsch and F. F. Wu, Stability regions of nonlinear autonomous dynamical systems, IEEE Trans. Automat. Control, 33 (1988), 16-27. doi: 10.1109/9.357.
[10]	M. Dellnitz and O. Junge, Set oriented numerical methods for dynamical systems, in Handbook of Dynamical Systems, Vol. 2, North-Holland, Amsterdam, 2002,221-264. doi: 10.1016/S1874-575X(02)80026-1.
[11]	F. Fallside and M. R. Patel, Step-response behaviour of a speed-control system with a back-e.m.f. nonlinearity, Proceedings of the Institution of Electrical Engineers, 112 (1965), 1979-1984. doi: 10.1049/piee.1965.0327.
[12]	O. P. Ferreira, M. S. Louzeiro and L. F. Prudente, Iteration-complexity of the subgradient method on Riemannian manifolds with lower bounded curvature, Optimization, 68 (2019), 713-729. doi: 10.1080/02331934.2018.1542532.
[13]	O. P. Ferreira and P. R. Oliveira, Subgradient algorithm on Riemannian manifolds, J. Optim. Theory Appl., 97 (1998), 93-104. doi: 10.1023/A:1022675100677.
[14]	P. Giesl, Construction of Global Lyapunov Functions Using Radial Basis Functions, Lecture Notes in Math. 1904, Springer, 2007.
[15]	P. Giesl, Converse theorems on contraction metrics for an equilibrium, J. Math. Anal. Appl., 424 (2015), 1380-1403. doi: 10.1016/j.jmaa.2014.12.010.
[16]	P. Giesl and S. Hafstein, Construction of a CPA contraction metric for periodic orbits using semidefinite optimization, Nonlinear Anal., 86 (2013), 114-134. doi: 10.1016/j.na.2013.03.012.
[17]	P. Giesl and S. Hafstein, Review of computational methods for Lyapunov functions, Discrete Contin. Dyn. Syst. Ser. B, 20 (2015), 2291-2331. doi: 10.3934/dcdsb.2015.20.2291.
[18]	P. Giesl, S. Hafstein and C. Kawan, Review on contraction analysis and computation of contraction metrics, J. Comp. Dyn., 10 (2023), 1-47. doi: 10.3934/jcd.2022018.
[19]	P. Giesl, S. Hafstein and I. Mehrabinezhad, Computation and verification of contraction metrics for exponentially stable equilibria, J. Comput. Appl. Math., 390 (2021), Paper No. 113332, 21 pp. doi: 10.1016/j.cam.2020.113332.
[20]	P. Giesl, S. Hafstein and I. Mehrabinezhad, Computing contraction metrics for three-dimensional systems, IFAC PapersOnLine, 54 (2021), 297-303. doi: 10.1016/j.ifacol.2021.06.151.
[21]	P. Giesl and H. Wendland, Construction of a contraction metric by meshless collocation, Discrete Cont. Dyn. Syst. Ser. B, 24 (2019), 3843-3863. doi: 10.3934/dcdsb.2018333.
[22]	S. F. Hafstein, A constructive converse Lyapunov theorem on exponential stability, Discrete Contin. Dyn. Syst., 10 (2004), 657-678. doi: 10.3934/dcds.2004.10.657.
[23]	W. Hahn, Stability of Motion, Springer, Berlin, 1967.
[24]	C. S. Hsu, Global analysis by cell mapping, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 2 (1992), 727-771. doi: 10.1142/S0218127492000422.
[25]	L. Jocić, Planar regions of attraction, IEEE Trans. Automat. Control, 27 (1982), 708-710. doi: 10.1109/TAC.1982.1102991.
[26]	T. A. Johansen, Computation of Lyapunov functions for smooth, nonlinear systems using convex optimization, Automatica, 36 (2000), 1617-1626. doi: 10.1016/S0005-1098(00)00088-1.
[27]	C. Kawan, S. F. Hafstein and P. Giesl, ResEntSG: Restoration entropy estimation for dynamical systems via Riemannian metric optimization, SofwareX, 15 (2021), 100743. doi: 10.1016/j.softx.2021.100743.
[28]	C. Kawan, S. Hafstein and P. Giesl, A subgradient algorithm for data-rate optimization in the remote state estimation problem, SIAM J. Appl. Dyn. Syst., 20 (2021), 2142-2173. doi: 10.1137/21M1406453.
[29]	H. Khalil, Nonlinear Systems, 3rd edition, Pearson, 2002.
[30]	N. Krasovskiĭ, Problems of the Theory of Stability of Motion, Mir, Moskow, 1959. English translation by Stanford University Press, 1963.
[31]	B. Krauskopf and H. M. Osinga, Computing Invariant Manifolds via the Continuation of Orbit Segments, 117-154, Springer Netherlands, Dordrecht, 2007. doi: 10.1007/978-1-4020-6356-5_4.
[32]	M. Krstic, I. Kanellakopoulos and P. Kokotovic, Nonlinear and Adaptive Control Design, Wiley, 1995.
[33]	W. Lohmiller and J.-J. E. Slotine, On contraction analysis for non-linear systems, Automatica, 34 (1998), 683-696. doi: 10.1016/S0005-1098(98)00019-3.
[34]	A. M. Lyapunov, The general problem of the stability of motion, Internat. J. Control, 55 (1992), 521-790. Translated by A. T. Fuller from Édouard Davaux's French translation (1907) of the 1892 Russian original, With an editorial (historical introduction) by Fuller, a biography of Lyapunov by V. I. Smirnov, and the bibliography of Lyapunov's works collected by J. F. Barrett, Lyapunov centenary issue. doi: 10.1080/00207179208934253.
[35]	G. Osipenko, Dynamical Systems, Graphs, and Algorithms, Lecture Notes in Math. 1889, Springer, 2007.
[36]	P. Parrilo, Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimiziation, PhD thesis: California Institute of Technology Pasadena, California, 2000.
[37]	M. M. Peet, Exponentially stable nonlinear systems have polynomial Lyapunov functions on bounded regions, IEEE Trans. Automat. Control, 54 (2009), 979-987. doi: 10.1109/TAC.2009.2017116.
[38]	M. M. Peet and A. Papachristodoulou, A converse sum of squares Lyapunov result with a degree bound, IEEE Trans. Automat. Control, 57 (2012), 2281-2293. doi: 10.1109/TAC.2012.2190163.
[39]	C. Udriste, Convex Functions and Optimization Methods on Riemannian Manifolds, vol. 297, Springer Science & Business Media, 2013.
[40]	M. Vidyasagar, Nonlinear System Analysis, 2nd edition, Classics in applied mathematics, SIAM, 2002. doi: 10.1137/1.9780898719185.
[41]	W. Walter, Ordinary Differential Equation, Springer, 1998. doi: 10.1007/978-1-4612-0601-9.
[42]	X. Wang, Subgradient algorithms on Riemannian manifolds of lower bounded curvatures, Optimization, 67 (2018), 179-194. doi: 10.1080/02331934.2017.1387548.
[43]	H. Zhang and S. Sra, First-order methods for geodesically convex optimization, JMLR: Workshop and Conference Proceedings, 49 (2016), 1-21, https://arXiv.org/abs/1602.06053.
[44]	V. Zubov, Methods of A. M. Lyapunov and their Application, Translation prepared under the auspices of the United States Atomic Energy Commission; edited by Leo F. Boron, P. Noordhoff Ltd, Groningen, 1964.