\begin{document}
\newcommand{\spc}{\ensuremath{\,\,}}
\def\red{\color{red}}
\def\bl{\color{blue}}
\def\ora{\color{orange}}
\def\green{\color{green}}
\def\br{\color{brown}}
\newcommand{\clist}{\texttt{c\_list}}
\newcommand{\calpha}{\texttt{c\_alpha}}
\newcommand{\slist}{\texttt{s\_list}}
\def\la{\langle}
\def\ra{\rangle}
\def\e{{\rm e}}
\def\x{\mathbf{x}}
\def\by{\mathbf{y}}
\def\bz{\mathbf{z}}
\def\F{\mathcal{F}}
\def\R{\mathbb{R}}
\def\T{\mathbf{T}}
\def\N{\mathbb{N}}
\def\K{\mathbf{K}}
\def\bK{\overline{\mathbf{K}}}
\def\Q{\mathbf{Q}}
\def\M{\mathbf{M}}
\def\O{\mathbf{O}}
\def\C{\mathbf{C}}
\def\P{\mathbf{P}}
\def\Z{\mathbb{Z}}
\def\H{\mathcal{H}}
\def\A{\mathbf{A}}
\def\V{\mathbf{V}}
\def\AA{\overline{\mathbf{A}}}
\def\B{\mathbf{B}}
\def\L{\mathbf{L}}
\def\bS{\mathbf{S}}
\def\H{\mathcal{H}}
\def\I{\mathbf{I}}
\def\Y{\mathbf{Y}}
\def\X{\mathbf{X}}
\def\G{\mathbf{G}}
\def\L{\mathbf{L}}
\def\D{\mathbf{D}}
\def\B{\mathbf{B}}
\def\f{\mathbf{f}}
\def\z{\mathbf{z}}
\def\y{\mathbf{y}}
\def\d{\hat{d}}
\def\bx{\mathbf{x}}
\def\y{\mathbf{y}}
\def\v{\mathbf{v}}
\def\g{\mathbf{g}}
\def\w{\mathbf{w}}
\def\b{\mathcal{B}}
\def\a{\mathbf{a}}
\def\q{\mathbf{q}}
\def\u{\mathbf{u}}
\def\s{\mathcal{S}}
\def\cc{\mathcal{C}}
\def\co{{\rm co}\,}
\def\cp{{\rm CP}}

\def\tg{\tilde{f}}
\def\tx{\tilde{\x}}
\def\supmu{{\rm supp}\,\mu}
\def\supnu{{\rm supp}\,\nu}
\def\m{\mathcal{M}}
\def\bR{\mathbf{R}}
\def\om{\mathbf{\Omega}}
\def\s{\mathcal{S}}
\def\k{\mathcal{K}}
\def\la{\langle}
\def\ra{\rangle}
\def\blambda{{\boldmath{\lambda}}}
\def\bsmlambda{\boldmath{\lambda}}
\def\ov{\overline{o}}
\def\und{\underline{o}}
\def\minlog{\texttt{Minlog}}
\newcommand{\comment}[3]{
\ifcomment
	{\color{#1}\bfseries\sffamily(#3)
	}
	\marginpar{\textcolor{#1}{\hspace{3em}\bfseries\sffamily #2}}
	\else
	\fi
}
\newcommand{\victor}[1]{
	\comment{blue}{V}{#1}
}
\newif\ifcomment
\commentfalse
\commenttrue
\title[Computer-assisted proofs for Lyapunov stability via SOS certificates]{Computer-assisted proofs for Lyapunov stability via Sums of Squares certificates and Constructive Analysis}

\author[G. Devadze \and V. Magron  \and S. Streif]{Grigory Devadze \and Victor Magron \and Stefan Streif}

\address{Grigory Devadze: Technische Universit\"at Chemnitz\\
Automatic Control and System Dynamics Lab\\
09107 Chemnitz, Germany}
\email{grigory.devadze@etit.tu-chemnitz.de}
\address{Victor Magron: LAAS, 7 avenue du Colonel Roche\\
31077 Toulouse C\'edex 4,France}
\email{vmagron@laas.fr}
\address{Stefan Streif: Technische Universit\"at Chemnitz\\
Automatic Control and System Dynamics Lab\\
09107 Chemnitz, Germany}
\email{stefan.streif@etit.tu-chemnitz.de}
\date{}

\begin{abstract}
We provide a computer-assisted approach to ensure that a given continuous or discrete-time polynomial system is (asymptotically) stable.
Our framework relies on constructive analysis together with formally certified sums of squares Lyapunov functions. The crucial steps are formalized within of the proof assistant $\minlog$. 
We illustrate our approach with various examples issued from the control system literature. 
\end{abstract}

\keywords{Lyapunov stability;sums of squares;  semidefinite programming relaxations; constructive analysis; proof assistants}

\subjclass{90C22}

\maketitle


\section{Introduction}
Computer-assisted theorem proving -- performed by special software called {\em proof assistants} -- have been successfully used for important problems, including the formalization of  Kepler's conjecture  \cite{pi17} and a proof of the odd order theorem in group theory \cite{gonthier2013machine}. One of the main advantages of using the proof assistants is their ability of incorporating the formal correctness of the logical statements programmatically. 

Lyapunov stability plays one of the most crucial roles in nonlinear control theory and the {\em second method of Lyapunov}, also called the {\em Lyapunov function method} is a standard tool for the controller design.
In this work we present a computer-assisted verification framework for the Lyapunov stability analysis of polynomial dynamical. 
 In the presented approach, first, Lyapunov functions are constructed via sums of squares techniques; second, the formal correctness is certified within of the framework of \texttt{Minlog}. The certificate guarantees that the constructed Lyapunov function meets its specification and the combination of both approaches allows for automation of reasoning on dynamical systems. The stability analysis requires dealing with nonlinear inequalities and in contrast to purely numerical methods \cite{Giesl2015-review}, we are interested in reasoning with the so-called exact real arithmetic. For this purpose the standard library of \texttt{Minlog} is extended in this work by several necessary definitions and theorems in the setting of constructive reals. 
 The presented concepts may also be viewed as a preliminary work for the development of formalized nonlinear control theory, hence offering the possibility of computer-assisted controller design and program extraction for controller implementation.
\subsection*{Related works}


In the context of dynamical systems,  \cite{boldo2013wave} formally proves the correctness of a program implementing the numerical resolution of a wave equation, with Frama-C \cite{cuoq2012frama}.  
In \cite{immler2019flow}, the authors formalize  variational equations related to the solutions of ordinary differential equations in Isabelle/HOL \cite{Nipkow02IsabelleHol}. 
Beyond formalization of numerical analysis procedures, several frameworks are currently available to perform computer-assisted proofs of mathematical analysis. 
The work in \cite{affeldt2020formalizing} proposes formalized mathematical structures in Coq \cite{CoqProofAssistant} for function  spaces. 
Earlier developments have been pursued towards the formalization of real analysis in several proof assistants, including Coq \cite{mayero2001formalisation}, HOL-Light,  \cite{harrison2012theorem}, Isabelle/HOL \cite{fleuriot2000mechanization} or  ACL2 \cite{gamboa2001nonstandard}, as well as constructive real analysis \cite{geuvers2000constructive,Krebbers2013-realanalysis}.

Proof assistants also play an important role in formal verification of control systems. 
Zou \cite{Zou2013-form-ver-Chinese-train} used Isabelle for implementation of their control system formal verification based on Hybrid Hoare logic. 
The Why3 platform \cite{bobot2011why3}  coupled with MATLAB/Simulink was used in \cite{Araiza-Illan2014-thm-prover-sys, Araiza2015-thm-proving-Simulink} to perform simple stability checks of linear discrete systems with quadratic Lyapunov functions. 
Gao \cite{Gao2014-descr-ctrl-thr} mentioned several proof assistant software tools to be considered for implementation of formal verification of control systems, including Coq \cite{CoqProofAssistant}, HOL-Light \cite{harrison1996hol}, Isabelle/HOL and Lean \cite{de2015lean}. The most closely related work is due to \cite{Anand2015-roscoq} where the authors present the framework for the certified programs for robots based on constructive real numbers \cite{Oconnor2008-certified}. Their work is based on the Logic of Events for distributed agent systems and generation of verified software controllers. However they do not address the fundamental principles of control theory and system dynamics like stability, robustness or optimality. In \cite{Rouhling2018-formal} the author formally proves the soundness of a control function for the inverted pendulum within the classical logic. The author in \cite{Rouhling2018-formal} also mentions the impossibility of generating and running programs with that framework.

In the present work the interactive proof system $\minlog$ \cite{schwichtenberg2011minlog,Berger2011-minlog} is used for verification. The development of $\minlog$ adresses proof-theoretical methods based on realizability and their applications to the verification and especially to the program extraction. The implementation of $\minlog$ is done in \texttt{Scheme}. The main reason which justifies the usage of $\minlog$ is its special treatment of the program extraction, i.e., in addition to the program extraction, $\minlog$ provides the proof of correctness of the extracted program. 
The machinery behind $\minlog$ is supported by the soundness theorem and the theory of computable functionals \cite{Berger2011-minlog,Berger2012-proofs}. 
The conversion of proofs to instances handled by the soundness theorem can be done automatically \cite{Berger2011-minlog,Miyamoto2013-program,Benl1998-minlog}. This kind of approach is not used in, e.g., Coq \cite{CoqProofAssistant}, where the program extraction is guided by the use of external tools. Due to the so-called normalization-and-evaluation technique, $\minlog$ achieves often more simpler and efficient terms, which matters in practice \cite{Berger2011-minlog}.
  
 Consider the autonomous continuous-time dynamical system 
 \begin{align} \label{def:dyn-sys}
 	\dot{\x} = f(\x)
 \end{align}
 where $f:D \rightarrow \mathbb{R}^n$ is a locally Lipschitz function with equilibrium point at origin and $D \subset \mathbb{R}^n$.
A related result~\cite{khalil2002nonlinear} states that the equilibrium of the system (\ref{def:dyn-sys}) is asymptotically stable if there is some positive definite homogeneous function $V(\x)$ such that the derivative of $V$ along the trajectories of the system is negative definite, or stated equivalently, that the opposite of this derivative is positive definite.
If we replace the later condition by the weaker condition that the derivative is nonnegative, then the system is stable.
In both cases, such a function $V$ is called a {\em Lyapunov function}.
For polynomial systems, i.e., when $f$ is a polynomial, these two sufficient positivity/nonnegativity conditions can be strengthened by computing a positive definite homogeneous function $V$ which is a {\em sum of squares} (SOS) of polynomials and such that $- \grad V \cdot f$ is SOS~\cite{phdParrilo}.
In this case, we refer to $V$ as an {\em SOS Lyapunov function}.

Proving stability with SOS polynomials lies in the research track of characterizing nonlinear systems through linear programs (LP), whose unknown variables are measures with support on the system constraints. 
This methodology was first introduced for static polynomial optimization problems~\cite{lasserre2001global,phdParrilo}. 
In~\cite{lasserre2001global}, Lasserre provides a  hierarchy of relaxations yielding a converging sequence of lower bounds for the minimum of a polynomial over a set of constraints also defined with polynomials. 
Each bound is obtained by solving a semidefinite program (SDP) \cite{vandenberghe1996semidefinite}, i.e., by computing the infimum of a linear function under the constraint that a symmetric matrix has only nonnegative eigenvalues. 

SDP has by its own many applications, for instance in combinatorial optimization~\cite{Gvozdenovic09}, control theory~\cite{BEFB94} or parameter estimation \cite{Streif2013-estimation}. 
Available software libraries such as SeDuMi~\cite{Sturm98usingsedumi}, SDPA~\cite{Yamashita10SDPA} or Mosek~\cite{moseksoft} can solve SDP with up to a few thousand variables.

The general idea consists of modeling a polynomial optimization problem as an infinite-dimensional LP over Borel measures. 
In practice, one approximates the solution of this LP with a converging hierarchy of semidefinite programming (SDP) primal-dual programs, called {\em moment-SOS} or Lasserre's  hierarchy.
The primal is a \emph{moment} problem, that is an optimization problem where variables are the moments of a Borel measure. The first moment is related to the mean values, the second moment is related to the variances, etc.  

The book~\cite{lasserre2009moments} provides a non-exhaustive list of various problems arising in applied mathematics, statistics and probability, economics, engineering which can be cast as instances of optimization problems over  moments. 

The dual is an SOS problem, where the variables are the coefficients of SOS polynomials. An example of SOS polynomial is $(\frac{1}{3} x - y)^2 + (x+y)^2$. 
Even though there exist positive polynomials which cannot be written with SOS decompositions, it turns out that when the set of constraints satisfies certain assumptions (slightly stronger than compactness) then one can represent positive polynomials with weighted SOS decompositions. 
These weights happen to be the polynomials defining the set of constraints. Among several similar powerful results for representation of positive polynomials, the one related to the dual of Lasserre's hierarchy is due to Putinar~\cite{Putinar1993positive}.
For dynamical systems, this hierarchy have been extended multiple times, for example to obtain converging approximations for optimal control problems~\cite{lasserre2008nonlinear}, reachable sets~\cite{magron2019semidefinite} or invariant measures~\cite{invsdp}.

However, such methods rely solely on numerical SDP solvers, implemented in double precision, thus output approximate nonnegativity certificates. 
In order to circumvent the numerical errors, one can use several algorithmic frameworks either in the univariate case~\cite{Chevillard11,univsos} or in the multivariate case~\cite{PaPe08,KLYZ08,magron2018exact}.
In particular, 
the rounding-projection algorithm from \cite{PaPe08} outputs a weighted rational SOS decompositon of polynomials lying in the interior of the SOS cone.
More recently, \cite{magron2018exact} has provided an alternative framework to
compute such decompositions. 
In a nutshell, this framework first outputs an
approximate SOS decomposition for a perturbation of the input polynomial with an
arbitrary-precision SDP solver. 
Then an exact decomposition is provided thanks to
the perturbation terms and a compensation phenomenon of the numerical errors.
The underlying algorithm is implemented in the RealCertify Maple package \cite{magron2018realcertify}.

Beyond the stability proofs of critical control systems, computing exact certificates of  nonnegativity is mandatory in many application fields, such as certified evaluation in computer arithmetics~\cite{Chevillard11} or 
formal verification of real inequalities~\cite{magron14} within proof 
assistants such as Coq~\cite{CoqProofAssistant} or HOL-Light~\cite{harrison1996hol}.


\subsection*{Contributions}

The aim of the present work is to provide  computer-assisted proofs for the stability of a given continuous/discrete-time polynomial system.  
For such a system, our framework consists of proving formally that a given function is an SOS Lyapunov function.
Our formal proofs of polynomial nonnegativity are handled with weighted SOS certificates. 
Since proof assistants such as \texttt{MinLog} have computational limitation, we rely on external tools that produce certificates, whose checking is reasonably easier from a computational point of view. 
These certificates are obtained with the RealCertify \cite{magron2018realcertify} tool available outside of the proof assistant \texttt{MinLog} and  verified inside. 
Our first contribution, presented in Section~\ref{sec:prelim_exact}, is to provide an algorithm, called \texttt{exactLyapunov} to compute exact weighted rational SOS Lyapunov functions for a given continuous-time polynomial system (\ref{def:dyn-sys}).
For the sake of clarity, we provide detailed explanation for the case of continuous systems only. Analogously, 	the results can be applied to the discrete-time systems if we require that the difference in the Lyapunov function along the trajectories is negative definite. 

Section \ref{sec:minlog} is dedicated to the presentation of the $\minlog$ formalization framework and its extension of the standard library. 

Our further contribution, provided in Section \ref{sec:formal}, is to derive a full  framework to formally prove the stability of such  systems, by verifying the output of  \texttt{exactLyapunov} within the $\minlog$ system.
We illustrate our formal framework with several  benchmarks in Section \ref{sec:benchs}.


\section{Preliminary Background and Results}
\label{sec:prelim}

\subsection{Positive polynomials, sums of squares}
\label{sec:prelim_sos}

Let $\mathbb{Q}$ (resp.~$\mathbb{R}$) be the field of rational (resp.~real) numbers.
With $d,n \in \N$, let us denote by $\R[\x]$ (resp. $\R[\x]_{d}$) the space of real-valued $n$-variate polynomials (resp. of degree at most $d$) in the variable $\x=(x_1,\ldots,x_n) \in \R^n$.
Given a real symmetric matrix $\G$, the notation $\G \succeq 0$ (resp. $\G \succ 0$) means that $\G$ is \emph{positive semidefinite} (resp. positive definite), that is $\G$ has only nonnegative (resp. positive) eigenvalues.

Let us note $\Sigma[\x]$ the set of sums of squares (SOS) of polynomials and $\Sigma[\x]_{2 d} := \Sigma[\x] \cap \R[\x]_{2 d}$ the set of SOS polynomials of degree at most $2 d$.
An $n$-variate polynomial $g$ belongs to $\Sigma[\x]_{2 d}$ if there exist a nonzero $r \in \N$ and $g_1,\dots,g_r \in \R[\x]_{d}$ such that $f = g_1^2 + \dots + g_r^2$.
An example of SOS polynomial is $g = 4 x_1^4 + 4
x_1^3 x_2 - 7 x_1^2 x_2^2 - 2 x_1 x_2^3 + 10 x_2^4 = (2 x_1 x_2 + x_2^2)^2 + (2
x_1^2 + x_1 x_2 - 3 x_2^2)^2$.
For each multi-index $\alpha = (\alpha_1,\dots,\alpha_n) \in \N^n$, we use the notation $\x^{\alpha} = x_1^{\alpha_1} \cdots x_1^{\alpha_n}$, and we write $g(\x) = \sum_{\alpha} g_{\alpha} \x^{\alpha}$.
The support of a polynomial $g$ written as above is the set $\{\alpha \in \N^n: g_\alpha \neq 0 \}$. 
For our previous example, the support is $\{(4,0), (3,1), (2,2) ,(1,3), (0,4) \}$.
We say that $g \in \R[\x]_{d}$ has {\em full support} if its support is $\N^n_d := \{ \alpha \in \N^n : \alpha_1 + \dots + \alpha_n \leq d \}$.

It is trivial that the condition of $g$ to be SOS ensures that $g$ is nonnegative over $\R^n$. 

The converse statement establishing the equivalence between nonnegative polynomials and SOS polynomials, is true only in the cases of univariate polynomials of even degree, quadratic polynomials and quartic polynomials in two variables. 

This equivalence was proven by Hilbert and we refer the interested reader to \cite{reznick2000some} and the references therein for more details.
The condition that $g \in \Sigma[\x]_{2 d}$ is equivalent to the existence of a symmetric matrix $\G \succeq 0$ such that $g(\x) = \v_d(\x)^T \G  \v_d(\x)$, where $\v_d(\x) = (1,x_1,\dots,x_n,x_1^2,x_1 x_2,\dots, x_n^2, x_1^d,\dots,x_n^d)$ is the vector of all monomials of degree at most $d$. 
The length of this vector is equal to $\binom{n+d}{n}$.

The matrix $\G$ is called a {\em Gram matrix} \cite{choi1995sums} associated to $g$ and  
computing $\G$ boils down to solving a {\em semidefinite program} (SDP for short) \cite{vandenberghe1996semidefinite}.
We refer to this technique as the {\em Gram matrix method}.
The SOS decomposition $g = g_1^2 + \dots + g_r^2$ is a by-product of a decomposition of the form $g = \v_d(\x)^T \L^T \D \L \v_d(\x)$, with $\L$ being a lower
triangular matrix with nonnegative real entries on the diagonal and
$\D$ being a diagonal matrix with nonnegative real entries, such that $\L^T \D \L = \G $.

\subsection{SOS Lyapunov functions}
\label{sec:prelim_lyapunov}
A positive definite polynomial (resp. positive definite form) is a polynomial (resp. homogeneous polynomial) which takes only positive values, except at the origin.

Given a vector of polynomials $f = (f_1,\dots,f_n)$, let us now consider a continuous polynomial system (\ref{def:dyn-sys}).
Let $d$ be the maximal degree of the polynomials $f_1,\dots,f_n$.
The Gram matrix method can be used to certify that such a system with equilibrium at the origin is stable by computing a {\em Lyapunov} polynomial function $V$. 
More precisely, one would like to obtain $V \in \R[\x]_{2 k}$ such that (1) $V$ is a positive definite form and (2)  $- \grad V \cdot f \in \R[\x]_{2 k-1+d}$ is nonnegative.

Regarding (1), since we expect $V$ to be a homogeneous polynomial, we look for a decomposition $V(\x) = \w_{k}(\x)^T \G_1 \w_{k}(\x)$, where $\w_k$ is the vector of all monomials of degree equal to $k$. 
Regarding (2), we assume that $d$ is odd, so that the degree of $- \grad V \cdot f$ is even, i.e., that $\ell := (2k-1+d)/2$ is a natural integer (otherwise $- \grad V \cdot f$ has odd degree and cannot take only nonnegative values).
To ensure that these two conditions (1) and (2) are fulfilled, we consider the following semidefinite feasibility problem of computing $V$, $\G_1$ and $\G_2$ such that:

\begin{equation}
\label{eq:lyapunovSOS}
\tag{SDP$_k$}
\begin{cases}
V(\x) \qquad  \qquad \quad \ \ = \w_{k}(\x)^T \G_1 \w_{k}(\x) \,, \\
- \grad V(\x) \cdot f(\x) = \v_{\ell}(\x)^T \G_2 \v_{\ell}(\x) \,,\\
V \in \R[\x]_{2 k} \,, \quad \G_1, \G_2 \succ 0  \,.
\end{cases}
\end{equation}

The optimization variables of the above  feasibility program~\eqref{eq:lyapunovSOS} are the coefficients of $V$ and the entries of the matrices $\G_1$ and $\G_2$.
Any feasible solution $(V,\G_1,\G_2)$ of this SDP ensures that $V$ is a positive definite form and $- \grad V \cdot f$ is nonnegative, yielding a stability proof for the initial system (\ref{def:dyn-sys}).

To solve~\eqref{eq:lyapunovSOS}, one relies on numerical solvers implemented in double precision,
for instance SeDuMi~\cite{Sturm98usingsedumi} and SDPA~\cite{Yamashita10SDPA}, or arbitrary-precision solvers such as~SDPA-GMP~\cite{Nakata10GMP}.
However, such numerical solvers provide approximate
nonnegativity certificates, so that the polynomial $V$, and the matrices $\G_1, \G_2$ do not exactly satisfy the equality constraints of \eqref{eq:lyapunovSOS}. 


\subsection{From approximate to exact SOS certificates}
\label{sec:prelim_exact}
To overcome these feasibility issues, one can rely on the certification framework from  \cite{magron2018exact}, implemented in the RealCertify software library \cite{magron2018realcertify}.
For certification purpose, we are mainly interested in polynomials with rational coefficients, which belong to the interior of the SOS cone $\Sigma[\x]$, denoted by $\mathring{\Sigma}[\x]$.
Given a polynomial $g \in \mathring{\Sigma}[\x]$ of degree $2d$ with full support, there always exists by \cite[Proposition 5.5]{choi1995sums} a positive definite Gram matrix $\G$ associated to $g$, i.e., such that $g = \v_d(\x)^T \G \v_d(\x)$. 
In this case, one can apply the so-called \texttt{intsos} algorithm, stated in \cite[Algorithm 1]{magron2018exact}, to compute exact weighted rational SOS decompositions for polynomials with rational coefficients, lying in $\mathring{\Sigma}[\x]$.
Namely, we obtain as output a decomposition $g = c_1 g_1^2 + \dots c_r g_r^2$, with $c_1,\dots,c_r$ being nonnegative rational numbers and $g_1,\dots,g_r$ being polynomials with rational coefficients.

For the sake of self-containedness, let us briefly recall how \texttt{intsos} works. 
The algorithm takes as input a polynomial $g$ of degree $2 d$ with rational coefficients and full support, and accuracy parameters $\delta$ and $\delta_{c}$. 
Then, the procedure consists of a perturbation step and an absorption step.
The perturbation step consists of finding a small enough positive rational number $\varepsilon$,  such that the perturbed polynomial $g_{\varepsilon}(\x) := g(\x) - \varepsilon \sum_{\alpha \in \N^n_d} \x^{2 \alpha}$ also lies in $\mathring{\Sigma}[\x]_{2 d}$.

Then, \texttt{intsos} calls an external SDP solver with accuracy $\delta$ to find an approximate Gram matrix $\tilde{\G} \succ 0$ associated to $g_{\varepsilon}$, such that $g_{\varepsilon} (\x) \simeq \v_d(\x)^T \tilde{\G} \v_d(\x)$. 
Here the size (and rank) of $\tilde{\G}$ is $r = \binom{n+d}{n}$.
Then, we obtain an approximate Cholesky's decomposition $\L \L^T$ of $\tilde{\G}$ with accuracy $\delta_c$, and for all $i=1,\dots,r$, we compute the inner product of the $i$-th row of $\L$ by 
$\v_d(\x)$ to obtain a polynomial $s_i$.
For large enough accuracy $\delta$ of the SDP solver and accuracy $\delta_c$ of the  Cholesky's decomposition, the coefficients of the remainder $u =g_{\varepsilon} - \sum_{i=1}^r s_i^2$ become small enough to be ``absorbed'': that is $u + \varepsilon \sum_{\alpha \in \N^n_d} \x^{2 \alpha}$ admits itself an SOS decomposition, and so does $g = u + \varepsilon \sum_{\alpha \in \N^n_d} \x^{2 \alpha} + \sum_{i=1}^r s_i^2$.
If \texttt{intsos} fails to perform this absorption step, then it repeats the same procedure after increasing the accuracy of the SDP solver and Cholesky's decomposition.

We now present our algorithm, called \texttt{exactLyapunov}, similar to \texttt{intsos}, to compute exact weighted rational SOS Lyapunov functions for a given continuous-time polynomial system (\ref{def:dyn-sys}).

\begin{algorithm}
\caption{\texttt{exactLyapunov}}
\label{alg:exactLyapunov}
\begin{algorithmic}[1]

\Require polynomials $f_1,\dots, f_n \in  \mathbb{Q}[\x]$ of degree less than $d$
\Ensure lists $\clist_1$, $\clist_2$ of numbers in $\mathbb{Q}$ and lists $\slist_1$, $\slist_2$ of polynomials in $\mathbb{Q}[\x]$
\State $f \gets (f_1,\dots, f_n)$
\State $k=2$ \Comment{the minimal degree of the SOS Lyapunov function}
\State ok := false
\While {not ok} \label{line:ki}
\If {\eqref{eq:lyapunovSOS} is strictly feasible, i.e., returns a solution $(V, \G_1, \G_2)$ with $\G_1, \G_2 \succ 0$} ok := true
\Else $\ k \gets k + 2$
\EndIf
\EndWhile \label{line:kf}
\State $\clist_1$, $\slist_1$ = \texttt{intsos}$(V)$ \label{line:V}
\State $\clist_2$, $\slist_2$ = \texttt{intsos}$(- \grad V \cdot f)$ \label{line:gradV}
\State \Return $\clist_1$, $\clist_2$, $\slist_1$, $\slist_2$
\end{algorithmic}
\end{algorithm}

\begin{prop}
\label{prop:exactLyapunov}
Let $f = (f_1,\dots, f_n)$ be a vector of rational polynomials of degree at most $d$, with $d$ being odd. 
Assume that \eqref{eq:lyapunovSOS} admits a strictly feasible solution for some $k \in \N$.
Then \texttt{exactLyapunov}$(f_1,\dots,f_n)$ terminates and returns an exact rational weighted SOS Lyapunov function for the continuous time system (\ref{def:dyn-sys}). 
\end{prop}

\begin{proof}
The loop going from line~\lineref{line:ki} to line~\lineref{line:kf} increases the search degree for the SOS Lyapunov function. 
By assumption, \eqref{eq:lyapunovSOS} admits a strictly feasible solution for some  large enough $k \in \N$, so this loop  eventually terminates. 
In particular, we obtain an SOS Lyapunov function $V$ of degree $2 k$ satisfying the constraints of  \eqref{eq:lyapunovSOS}.
Then, one applies the perturbation-absorption procedure \texttt{intsos} to obtain an exact rational weighted SOS decomposition of $V$ at  line~\lineref{line:V}, and of $(- \grad V \cdot f)$ at line~\lineref{line:gradV}.
By \cite[Proposition 3.5]{magron2018exact},  the procedure \texttt{intsos} always terminates in both cases, since $V$ and $- \grad V \cdot f$ have associated positive definite Gram matrices. 
The outputs $\clist_1$ and $\slist_1$ respectively contain a list of nonnegative rational numbers $[c_1,\dots,c_r]$ and rational polynomials $[s_1,\dots,s_r]$ such that $V = c_1 s_1^2 + \dots + c_r s_r^2$. 
Similarly, we obtain an exact rational weighted SOS decomposition for $(- \grad V \cdot f)$.
\end{proof}


\section{Minlog system}
\label{sec:minlog}
In this section, we investigate how to formally prove stability of polynomial systems, by verifying the output of \texttt{exactLyapunov} within $\minlog$. We first introduce some preliminary concepts of the proof system \texttt{Minlog}, the notion of totality and the concept of inductive definitions, before we give the details how the formalization works in constructive setting. In this work the analysis and formalization are done within of the framework of constructive analysis with special focus on the concrete and efficient computational realization as well as formal verification. The constructive analysis by itself is a mathematical analysis, which deals with computational content of mathematical proofs in an informal manner. In order to obtain the certified programs there is no way around of the usage of higher-order programming languages \cite{Buss1998-handbook,geuvers2000constructive,fleuriot2000mechanization,harrison2012theorem,Krebbers2013-realanalysis}.  
\subsection{Necessary Preliminaries}
One of the possible interpretations of constructive analysis is due to Brouwer, Heyting and Kolmogorov (sometimes called BHK interpretation). 
Constructive analysis states the logical rules in a way which can be algorithmically interpreted as a mathematical proof \cite{Buss1998-handbook}. Let $P$ and $Q$ be some arbitrary formulas. The most important BHK rules are:
\begin{itemize}
	\item Proving the implication $P \Rightarrow Q$ means finding an algorithm that transforms a proof of $P$ into a proof of $Q$.
	\item To prove $\exists \gamma \spc P[\gamma]$, an algorithm computing $\gamma$ and the proof that $P[\gamma]$ holds are required.
	\item To prove $\forall \gamma \in A \spc P[\gamma]$, an algorithm is required that converts $\gamma \in A$ into a proof of $P[\gamma]$.
\end{itemize}
In practice, proof assistant software are  concrete realizations of the BHK interpretation. 

The interactive proof system \texttt{Minlog} has been developed with the special intention of  extracting programs from mathematical proofs. The computational model consists of the theory of the partial continuous functionals in the sense of Scott-Ershov 
and implements the theory of computable functionals \cite{Scott1982-domains}.
The atomic constructions are done by introduction of free algebras as types which are given by its constructors. 
For example, the free algebra for natural numbers is given by the constructors zero and successor, such that all other constructed elements  are also natural numbers.


Generally the types are build from the object types, type variables and type constants using algebra type formation and arrow type formation.
The later arrow type formation is denoted by ``$\rho\Rightarrow\sigma$'', where $\rho$ and $\sigma$ are two types.
The object types are usually types constructed by some free algebra variables. The type variables are specific variables of an unknown type. Consider for example a free algebra of lists. 
It is clear that the basic list manipulation like appending an element to some list or getting an element from the list should work independently of the type of the elements. 
Therefore the list manipulations are written in terms of type variables. 
The arrow type formation  $\rho\Rightarrow\sigma$ constructs a new type from the previously defined base types $\rho$ and $\sigma$, e.g., the successor of a natural number of type $\texttt{nat}$ is viewed as a function of type $\texttt{nat}\Rightarrow\texttt{nat}$.
The computation is implemented efficiently by the usage of the so-called normalization-by-evaluation. The reason of usage of typed theories is that the normalization-by-evaluation technique, which is applied to the proof terms, needs typed formulas.  As an example consider the construction of the type \texttt{nat} of natural numbers $\mathbb{N}$ by stating that \texttt{Zero} is of type \texttt{nat} and \texttt{Succ} is a mapping of type \texttt{nat$\Rightarrow$nat}. The normalization of the term 1+2 would result in \texttt{Succ(Succ(Succ Zero))} which is clearly of type $\texttt{nat}$ . The following types for the numbers with constructors and destructors are contained in the standard \texttt{Minlog} library:
\begin{itemize}
	\item \texttt{pos:One,SZero,SOne} (positive numbers $\mathbf{P}$ (binary))
	\item \texttt{nat:Zero,Succ} (natural numbers $\mathbf{N}$ )
	\item \texttt{list:Nil,Cons} 
	\item \texttt{int:IntP,IntZero,IntN} (integer numbers $\mathbf{Z}$)
	\item \texttt{rat:RatConstr} and destructors \texttt{RatN,RatD} (rational numbers $\mathbf{R}$)
	\item \texttt{rea:RealConstr} and destructors\\ \texttt{RealSeq,RealMod} (real numbers $\mathbf{R}$) 


\end{itemize} 
From this point onwards we use boldface symbols for numbers whenever we want to reason formally in \minlog.


 The  user-defined computation and corresponding rewriting rules are assigned to the so-called program constants. Usually the program constants are recursively written typed programs. Consider for example the multiplication of two natural numbers. The corresponding program constant \texttt{NatTimes} is of type \texttt{nat$\Rightarrow$nat$\Rightarrow$nat} and the computational rules are given recursively as:
\begin{align*}
&\texttt{NatTimes}\spc  n \spc \texttt{Zero} \rightarrow  \texttt{Zero} \\
&\texttt{NatTimes} \spc n \spc \texttt{Succ} \spc m \rightarrow \texttt{NatTimes}\spc n \spc m + n.
\end{align*}
As an additional example consider the strict inequality relation $<$ again for the natural numbers. The program constant \texttt{NatLt} is of type \texttt{nat$\Rightarrow$nat$\Rightarrow$bool} and the computational rules are:
\begin{align*}
&\texttt{NatLt}\spc   n \spc \texttt{Zero} \rightarrow  \texttt{False} \\
&\texttt{NatLt} \spc 0 \spc \texttt{Succ} \spc n \rightarrow \texttt{True}\\
&\texttt{NatLt} \spc \texttt{Succ} \spc n \spc \texttt{Succ} \spc m \rightarrow \texttt{NatLt} \spc  n \spc m.
\end{align*}
Usually the computer proofs incorporate sometimes trivial goals and therefore must be proven automatically. Thus the so-called rewriting rules are important for the automation of theorem proving. For example assume that one has proved that $\forall a \in \mathbf{Q} \spc a+0=a$. Then one can add the machine-generated computation rule $a+0 \rightarrow a$ to the library of  rewriting rules. Every time the normalization process is done, if some rational $a+0$ occurs then  the zero is canceled automatically.
\subsection{The notion of totality}
In order to establish the connection between proofs and programs, each program constant in \texttt{Minlog} has to fulfill the totality requirement, i.e., the termination of the program constant under any correct input. For a comprehensive description of totality, we refer the interested reader to \cite[Chapter~8.3]{Stoltenberg1994-domain} and \cite{Berger1993-total}.   


 In order to use the whole machinery of the computational model in \texttt{Minlog}, program constants are required to be total functions. The totality is usually proved by induction directly after the definition of the program constant. For example assume that for some natural numbers $n,m$, the program constant \texttt{NatPlus n m} is recursively defined as:
\begin{itemize}
	\item \texttt{NatPlus n Zero $\rightarrow$ $n$} 
	\item \texttt{NatPlus n Succ m $\rightarrow$ Succ(NatPlus n m)}.
\end{itemize}
The totality goal is: $\forall \texttt{n,m} \spc \texttt{TotalNat n} \rightarrow \texttt{TotalNat m} \rightarrow \texttt{TotalNat NatPlus n m}$. By induction on $\texttt{m}$ and by definition $\phantom{v}\texttt{TotalNat(NatPlus n Zero)} \phantom{v}$ is equivalent to $\texttt{TotalNat(n)}$ which follows from the assumptions. The induction step 
is $\texttt{TotalNat(n +Succ m)}$ which is equivalent by definition to $\texttt{TotalNat(Succ(n +m)}$. Finally the proof tree is split into two goals $\texttt{TotalNat(Succ(n)}$ and $\texttt{TotalNat(n+m)}$. Clearly the successor of a natural number is total and the second goal holds due to the induction hypothesis.
In the present work, the totality proofs of all program constants have been shown and implemented.

\subsection{Inductive definitions and theorem proving}
Before one starts to prove theorems with the help of the computer, one has to provide  definitions which are relevant for the problem statement. In \texttt{Minlog} it is done by introducing the so-called inductively defined predicate constant (IDPC). Informally speaking, these are formations of implications leading to the proposition to be defined. As an example assume that one wants to provide a definition for rational numbers which lie in a unit interval. Then the corresponding IDPC \texttt{RatInUnitI} is given by:
\begin{equation*}
\forall a \in \mathbf{Q} \spc 0 \le a \rightarrow a \le 1 \rightarrow \texttt{RatInUnitI} \spc a
\end{equation*}
which means semantically that to show $ \texttt{RatInUnitI} \spc a$ one has to show that $0 \le a$ and $a \le 1$. Additionally the definitions can also be provided inductively. As an example consider the definition of an even natural number given by the IDPC \texttt{Even}. Inductively one can proceed as follows:
\begin{align*}
&\texttt{Even Zero}\\
&\forall n \left( \texttt{Even} \spc n \rightarrow \texttt{Even} \spc n+2\right).
\end{align*}
Here the definition of the even natural number is given in terms of other natural number (in this case \texttt{Zero}) which is known to be even.

\subsection{Preliminaries for the formalization in constructive setting}
This section is based on the foundation of  constructive analysis with witnesses \cite{Schwichtenberg2006-constructive}. The witnesses are algorithms that certify certain statements, which are propagated through the mathematical proofs. For example, informally, a witness for a  strictly increasing continuous function $f:I \rightarrow \mathbb{R}$  would be a map $\tau: \mathbb{Q} \rightarrow \mathbb{Q}$ such that for any $\varepsilon > 0$ and $a,b \in I$, one has  $ \spc 0<\tau(\varepsilon)<f(x)-f(y)$ whenever $0<\varepsilon<y-x$. By relying on the axiom of choice \cite[p.~89]{Ye2011-SF} or the dependent choice axiom \cite{Schwichtenberg2006-inverting}, a continuous inverse function $f^{-1}$ can be constructed in terms of the witness $\tau$.

All presented concepts are implemented within of \texttt{Minlog}, and   furthermore we extended the standard library of \texttt{Minlog}  by the necessary notions and theorems for vectors, continuous functions on multiple variables, vector-valued functions as well as auxiliary theorems for real numbers.
The  most important object in the current setting is real number. In \texttt{Minlog} the reals are build by introducing type variables \texttt{as:$\mathbf{N}\Rightarrow\mathbf{Q}$} , \texttt{M:$\mathbf{P}\Rightarrow\mathbf{N}$}, defining the constructor by the arrow type formation \texttt{RealConstr:as$\Rightarrow$M$\Rightarrow\mathbf{R}$} and introducing the IDPC \texttt{Real}, i.e., a real $x$ is given by a regular Cauchy sequence of rationals $(a_n)_n$ (encoded by the type variable \texttt{as}), denoted as \texttt{Cauchy},  with a given monotone modulus $M$, denoted as \texttt{Mon}, that is $\forall p,m,n \spc |a_{n} - a_{m}|\le 2^{-p}$ whenever $M(p) \le n,m$ and $\forall p,q \spc p \le q \rightarrow M(p) \le M(q)$. In order to work with the IDPCs properly, the so-called elimination theorems have to be proven, which roughly speaking are deconstructors of the IDPCs and may be considered as an analogue of an informal statement: "By the definition of  proposition $P$, proposition $Q$ holds". The next theorem is an important example for such elimination theorems:     
\begin{thm}[\texttt{CauchyElim}] \label{thm:CauchyElim}
	\[
	\forall (a_n)_n,M \spc \mbox{\texttt{Cauchy}} \spc (a_n)_n \spc M \rightarrow \forall p,m,n \spc |a_{n} - a_{m}|\le 2^{-p} \,.
	\]
\end{thm}


A real $x=((a_n)_n,M)$ is called nonnegative if $\forall p \spc 0 \le a_{M(p)} + 2^{-p}$. For two real numbers $x$,$y$ denote $x \le y$ whenever the real $y-x$ is nonnegative. 
For a real $x$ denote $x \in_p \mathbf{R}^+ $ (alternately $0<_p x$) if $2^{-p} \le a_{M(p+1)}$. Two reals $x=((a_n)_n,M)$ and $y=((b_n)_n,N)$ are equal, if 
$\forall p\in \mathbf{P}\spc|a_{M(p+1)}-b_{N(p+1)}|\le 2^{-p}$. Note that the equality of the reals is undecidable. 
Two reals $x$,$y$ are called apart if $0<_p|x-y|$ which is usually denoted by $x\#y$.
The next two characterization theorems are often useful for operating with inequalities within  \texttt{Minlog}. 
\begin{thm}[\texttt{RealNNegChar}]\label{RealNNegChar}
	For a real $x=((a_n)_n,M)$ the following are equivalent:
	\begin{enumerate}
		\item $0 \le x$
		\item $ \forall p \spc \exists n_0 \spc \forall n\ge n_0 \spc -2^{-p}\le a_n$.
	\end{enumerate}
\end{thm}


\begin{thm}[\texttt{RealPosChar}]\label{RealPosChar}
	For a real $x=((a_n)_n,M)$ the following are equivalent:
	\begin{enumerate}
		\item $\exists p \spc x \in_p \mathbf{R}^+$ 
		\item $\exists q,p \spc \forall n\ge p \spc 2^{-q}\le a_n$.
	\end{enumerate}
\end{thm}


The next definition is the one for continuous functions. Note that a continuous function in the present context attains only rational values for  efficiency reasons. Later on, the application procedure will allow us to generate real numbers from such continuous functions.
\begin{defn}[\texttt{Cont}]\label{def:Cont}
	A uniformly continuous function $f:I \rightarrow \mathbf{R}$ on a compact interval $I$ with rational endpoints is given by
	\begin{enumerate}
		\item an approximation map $h_f :  (I \cap \mathbf{Q} ) \times \mathbf{N} \rightarrow \mathbf{Q} $ and a map $\alpha_f: \mathbf{P} \rightarrow \mathbf{N}$ such that  $(h_f(a,n))_n$ is a Cauchy sequence with modulus $\alpha_f$ for any $a \in (I \cap \mathbf{Q} )$
		\item a modulus  of continuity $\omega_f:\mathbf{P} \rightarrow \mathbf{P}$, which satisfies 
		\begin{equation*}
		|a-b| \le 2^{-\omega_f(p) +1} \rightarrow |h_f(a,n)-h_f(b,n)|\le 2^{-p}
		\end{equation*}		
	\end{enumerate}
	 for any $a,b \in (I \cap \mathbf{Q} )$ with monotone $\alpha_f$ and $\omega_f$.
\end{defn}


As an example, consider the representation of the function $f(x)=x^2$ on the interval $I = [0,2]$ in \texttt{Minlog}:\\
\texttt{ContConstr 0 2 ([a,n]a*a) ([p] 0) ([p] p+3)}
where \texttt{ContConstr} is a constructor of a function of type \texttt{cont} 
 which takes as the input the rational endpoints of the interval $I = [0,2]$, the approximating map \texttt{([a,n]a*a)} of type \texttt{$\left(\mathbf{Q} \Rightarrow \mathbf{N} \right) \Rightarrow \mathbf{Q} $}, the mapping $\alpha_f(p):=0$ of  type \texttt{$\mathbf{P} \Rightarrow \mathbf{N} $} as  modulus of convergence and the mapping $\omega_f(p):=p+3$ of type \texttt{$\mathbf{P} \Rightarrow \mathbf{P} $} as modulus of continuity. Along with the constructor, a proof of \texttt{Cont} has to be provided for this representation. That is, according to the definition of \texttt{Cont}
\begin{itemize}
	\item \texttt{([a,n]a*a)} is \texttt{Cauchy} with the modulus $\alpha_f(p):=0$ because the values of the sequence $(h_f(a_n,n))_n:=(a^2)_n$ are independent of $n$.
	\item The modulus of continuity $\omega_f(p):=p+3$ is correct and the corresponding goal  \[\forall a \in [0,2] \spc b \in [0,2] \spc p \in \mathbf{P} \spc \left( |a-b| \le 2^{-p-2} \rightarrow |a \cdot a-b \cdot b| \le 2^{-p}\right)\]
	holds because $|a \cdot a - b \cdot b| \le |(a-b)(a+b)| \le 2^{-p-2}\cdot|(a+b)| \le 2^{-p-2} \cdot 4 \le 2^{-p}$.\\
	Note that $(a+b) \le 4$ holds because   $a,b \in [0,2]$.
	\item $\alpha_f : p \mapsto 0$ and $\omega_f : p \mapsto p+3$ are both monotone.
\end{itemize}

 Note that the definition of the approximation map $h_f$ in \texttt{Cont} differs from the definition in \texttt{Minlog} in the sense that the approximation maps in \texttt{Minlog} are defined as programs of type $\left(\mathbf{Q}\Rightarrow \mathbf{N}\right) \Rightarrow \mathbf{Q}$, omitting the information of the interval $I$. However this information is provided on the proof level by attaching the assumptions that the rational inputs have to belong to the specific interval $I$. 
 

\begin{defn}[\texttt{Application}]
	The application of a uniformly continuous function $f:I\rightarrow \mathbf{R}$ to a real $x=((a_n)_n,M))$ in $I$, is defined to be a Cauchy sequence
	$(h_f(a_n,n))_n$ with modulus $\max(\alpha_f(p+2),M(\omega_f(p+1)-1).$
\end{defn}
\begin{thm}[\texttt{ApplicationCorr}]\label{thm:ApplicationCorr}
	For a uniformly continuous function $f:I\rightarrow \mathbb{R}$ and a real $x=((a_n)_n,M))$ in $I$, the Cauchy sequence emanating from \texttt{Application} and accompanied with the modulus $\max(\alpha_f(p+2),M(\omega_f(p+1)-1)$ is a real number.
\end{thm}

\begin{rem}
	The definition of the application incorporates the so-called extension theorem due to Bishop \cite[p.~91]{Bishop1985-constr-analysis}. In \texttt{Minlog} the application is realized by a program constant and allows one to address the lifting from rational mappings (see Definition \ref{def:Cont}) to the reals, which is at some point similar to the technique of the completion monad for the constructive reals in \texttt{Coq} \cite{OConnor2007-monad}. Additionally, note that for the continuous function $f:I \rightarrow \mathbf{R}$ and real numbers $x,y \in I$ under the usage of the program \texttt{Application}, one is able to prove:
	\[
		|x-y| \le 2^{-\omega_f(p)} \rightarrow |f(x)-f(y)|\le 2^{-p} \,,
	\]
	which coincides to another usual existing notions of continuity in effective analysis \cite{Bishop1985-constr-analysis,Ko1991,Kushner1999-markov,Weihrauch2012-computable-analys}.
\end{rem}
The vectors, functions on multiple variables and vector-valued functions can be formalized in a straightforward way by utilizing the algebra of lists. This representation allows to unfold the vector elements and their properties inductively.

\begin{defn}[Free algebra \texttt{list}]
		A list \texttt{xls} containing the elements of type $\alpha$ is either :
		\begin{itemize}
			\item \texttt{(Nil alpha)}, also called empty,
			\item constructed \texttt{xl::xls} with head \texttt{xl} and tail \texttt{xls}.
		\end{itemize}
\end{defn}
\begin{rem}
	A list contaning one element \texttt{xl} is denoted by $\texttt{xl:}$.	
\end{rem}
\begin{defn}[\texttt{CauchyVector}]
		Consider the lists $v:\texttt{list }$$(\mathbf N \Rightarrow \mathbf Q) $, $M_v:\texttt{list }$$(\mathbf P \Rightarrow \mathbf N) $. A pair $\left(v,M_v\right)$ is a \texttt{CauchyVector} if $|v|=|M_v|$ and either:
		\begin{itemize}
			\item $v$ and $M_v$ are empty,
			\item the pair (\texttt{head $v$}, \texttt{head $M_v$}) is \texttt{Cauchy} and the pair (\texttt{tail $v$}, \texttt{tail $M_v$}) is \texttt{CauchyVector}.
		\end{itemize}
\end{defn}

\begin{defn}[\texttt{MonVector}]
	A list $M_v:\texttt{list }$$(\mathbf P \Rightarrow \mathbf N) $ is  \texttt{MonVector} if either:
	\begin{itemize}
		\item $M_v$ are empty,
		\item the map \texttt{head $M_v$} is \texttt{Mon} and the list \texttt{tail $M_v$} is \texttt{MonVector}.
	\end{itemize}
\end{defn}

\begin{defn}[\texttt{RealVector}]
	For the lists $v:\texttt{list }$$(\mathbf N \Rightarrow \mathbf Q) $, $M_v:\texttt{list }$$(\mathbf P \Rightarrow \mathbf N) $ of length $n$,  a vector $\bold x=(v,M_v) \in \mathbf{R}^n$ if the pair $(v,M_v)$ is \texttt{CauchyVector} and $M_v$ is \texttt{MonVector}.
\end{defn}


Similarly to Definition \ref{def:Cont}, a function on multiple variables can be defined as follows:
\begin{defn}[\texttt{ContMV}]
	A uniformly continuous function of multiple variables $g:B \rightarrow \mathbf{R}$ on a compact ball $B=(c,R) \subset \mathbf{R}^n$ with rational center $ c$ and rational radius $R$ is given by
	\begin{enumerate}
		\item an approximation map $h_{g}:(B \cap \mathbf{Q}^n ) \times \mathbf{N} \rightarrow \mathbf{Q} $ and a map $\alpha_{g}: \mathbf{P} \rightarrow \mathbf{N}$ such that $(h_{g}( e,n))_n$ is a Cauchy sequence with modulus $\alpha_g$ for any $e \in (B \cap \mathbf{Q}^n )$
		\item a modulus $\omega_{g}:\mathbf{P} \rightarrow \mathbf{P}$ of continuity, which satisfies 
		\begin{equation*}
		\|e_1-e_2\| \le 2^{-\omega_{g}(p) +1} \rightarrow |h_{g}(e_1,n)-h_{g}(e_2,n)|\le 2^{p}
		\end{equation*}		
	\end{enumerate}
	for any $e_1,e_2 \in (B \cap \mathbf{Q}^n )$ with monotone $\alpha_{g}$ and $\omega_{g}$.
\end{defn}
\begin{rem}
	The choice of the norm $\| \cdot \|$ influences the values of moduli $\alpha_g$ and $\omega_g$. Usually, one would prefer to use either $1$-norm or $\infty$-norm for  computational purpose, since for rational vectors the values of the norm are also rationals, which is not the case, e.g., for the Euclidean norm.
\end{rem}
In summary, in this section we presented the most important concepts of the \texttt{Minlog} system and presented the formalization of constructive analysis in \texttt{Minlog} as well as the extension of the standard library which was necessary for our goal. In the next section we provide the proof-theoretical verification of the SOS Lyapunov functions in our framework.
\section{Formal Proofs of Stability for Polynomial Systems}
\label{sec:formal}
Once the SOS Lyapunov function has been obtained by  Algorithm \ref{alg:exactLyapunov}, the proof of the final formal verification goal involves dealing with nonlinear inequalities. In this section we make a brief remark on the general undecidability of the problem as well as provide technical details about the implementation of polynomials in \texttt{Minlog}. Finally we present the developed framework for the formal verification of Lyapunov functions by  introducing various inductive definitions for positive definite and nonnegative functions and by establishing the proof theoretical connections between these definitions.
  
It is an easy observation that proving the (asymptotic) stability for a general dynamical system is undecidable and the problem can be reduced to the statement of  Richardson's Theorem \cite{Richardson1969-undecidable}.

\begin{thm}[Richardson]
	Let $E$ be a  set of expressions representing real, single valued, partially defined functions of one real variable. Let $E^*$ be the set of functions represented by expressions in $E$. If $E^*$ contains $\log 2, \pi,e^x,\sin x$ then the problem ``{\textbf {given an expression $A$ in $E$, decide if there is a real number $x$ such that $A(x)$ is less than zero}}'' is undecidable. 
\end{thm}

\begin{rem}
	 Richardson's Theorem has important implications for computer algebra systems since it shows that for a certain class of functions the problem of deciding if there is a canonical form for a symbolic expression \cite{Caviness1970-canonical} or if there exists an elementary undefined integral of an elementary function \cite{Risch1970-integration} are both undecidable.
\end{rem}

\begin{rem}
	Usually when working with polynomials, one has to provide a proper representation of polynomials which encodes the coefficients and degree. Note that in the case of real number coefficients the degree of the polynomial is not computable \cite{Pauly2018-comparing}. In this work, the polynomials are given directly as the continuous functions with polynomial mapping $h_f$ with rational coefficients and by the modulus of convergence $\alpha(p)=0$. Thus such continuous functions are marked as polynomials and normalization-by-evaluation allows one often to work directly with the polynomial mapping $h_f$. Due to Theorem \ref{thm:ApplicationCorr} one is still able to address real number coefficients or inputs for the polynomials in this framework.
\end{rem}

Now in order to incorporate the notions of positive definite and nonnegative functions for the verification of the output of Algorithm \ref{alg:exactLyapunov} (\texttt{exactLyapunov}), the appropriate definitions and the equivalent IDPCs are provided as follows:
\begin{defn}[Nonnegative function (\texttt{ContNonNeg})]\label{def:NonNeg}
	A uniformly continuous function $g:B=(\bold c,R) \rightarrow \mathbb{R}$ given by $(h_{g},\alpha_{g},\omega_{g},(\bold c,R))$ is called nonnegative whenever 
	\begin{enumerate}
		\item $\|\bold c \| \le R$,
		\item $\forall \bold x \spc 0 \le g(\bold x)$,
	\end{enumerate}
\end{defn}
and the corresponding IDPC in \texttt{Minlog}:
\if{
\begin{minted}[linenos,frame=single,framesep=5pt]{scheme}
(add-ids
(list (list "ContNonNeg" (make-arity (py "contmv"))))
'("all fmv(ContMV g -> all xx(xx dim = g dimmv ->
RealVector xx -> InDom xx g -> 
 0 <<= (g xx) ) -> ContNonNeg g)" 
"ContNonNegIntro"))
\end{minted}
}\fi

\begin{Verbatim}[commandchars=\\\{\}]
\PYG{p}{(}\PYG{n+nf}{add\PYGZhy{}ids}
\PYG{p}{(}\PYG{n+nb}{list }\PYG{p}{(}\PYG{n+nb}{list }\PYG{l+s}{\PYGZdq{}ContNonNeg\PYGZdq{}} \PYG{p}{(}\PYG{n+nf}{make\PYGZhy{}arity} \PYG{p}{(}\PYG{n+nf}{py} \PYG{l+s}{\PYGZdq{}contmv\PYGZdq{}}\PYG{p}{))))}
\PYG{o}{\PYGZsq{}}\PYG{p}{(}\PYG{l+s}{\PYGZdq{}all fmv(ContMV g \PYGZhy{}\PYGZgt{} all xx(xx dim = g dimmv \PYGZhy{}\PYGZgt{}}
\PYG{l+s}{RealVector xx \PYGZhy{}\PYGZgt{} InDom xx g \PYGZhy{}\PYGZgt{}}
\PYG{l+s}{ 0 \PYGZlt{}\PYGZlt{}= (g xx) ) \PYGZhy{}\PYGZgt{} ContNonNeg g)\PYGZdq{}}
\PYG{l+s}{\PYGZdq{}ContNonNegIntro\PYGZdq{}}\PYG{p}{))}
\end{Verbatim}
\begin{rem}
	Due to the undecidability, the nonnegativity on real numbers is realized via the inductively defined predicate as opposed to the inequality on rational numbers, which is implemented as a simple program constant, by the following IDPC:
	\[
	\texttt{RealNNeg: $\forall x$ Real $x \rightarrow$  $\forall p$ 0 $\le$ $x$ seq($x$ mod $p$)$+2^{-p} \rightarrow$ RealNNeg $x$ } 
	\]
	in \minlog, where \texttt{seq} and \texttt{mod} are the display tokens for the deconstructors of the type \texttt{real}. 
\end{rem}

\begin{defn}[Positive Definite Function (\texttt{PosDef})]\label{def:PosDef}
	A uniformly continuous function $g:B=(c,R) \rightarrow \mathbb{R}$ given by $(h_{g},\alpha_{g},\omega_{g},(c,R))$ is called positive definite whenever 
	\begin{enumerate}
		\item $\|\bold c \| \le R$,
		\item $\bold x=0 \rightarrow g(\bold x) = 0$,
		\item $\forall \bold x \, (\exists \spc p \spc 0 \#_p \spc \bold x \rightarrow \exists \spc q \spc 0 <_q \spc g( \bold x))$,
	\end{enumerate}
\end{defn}
and the corresponding IDPC in \texttt{Minlog}:
\if{
\begin{minted}[linenos,frame=single,framesep=5pt]{scheme}
(add-ids
(list (list "PosDef" (make-arity (py "contmv"))))
'("all g(ContMV g -> (|_ g defC <<= g defB 
&0===g (ZerosReal (g defC dim)) )&
all xx(RealVector xx -> InDom xx g -> exl p RealLt 0 (|_ xx) p -> 
exl q(RealLt 0 (g xx) q)) -> PosDef g)" 
"PosDefIntro"))	
\end{minted}
}\fi

\begin{Verbatim}[commandchars=\\\{\}]
\PYG{p}{(}\PYG{n+nf}{add\PYGZhy{}ids}
\PYG{p}{(}\PYG{n+nb}{list }\PYG{p}{(}\PYG{n+nb}{list }\PYG{l+s}{\PYGZdq{}PosDef\PYGZdq{}} \PYG{p}{(}\PYG{n+nf}{make\PYGZhy{}arity} \PYG{p}{(}\PYG{n+nf}{py} \PYG{l+s}{\PYGZdq{}contmv\PYGZdq{}}\PYG{p}{))))}
\PYG{o}{\PYGZsq{}}\PYG{p}{(}\PYG{l+s}{\PYGZdq{}all g(ContMV g \PYGZhy{}\PYGZgt{} (|\PYGZus{} g defC \PYGZlt{}\PYGZlt{}= g defB}
\PYG{l+s}{\PYGZam{}0===g (ZerosReal (g defC dim)) )\PYGZam{}}
\PYG{l+s}{all xx(RealVector xx \PYGZhy{}\PYGZgt{} InDom xx g \PYGZhy{}\PYGZgt{} exl p RealLt 0 (|\PYGZus{} xx) p \PYGZhy{}\PYGZgt{}}
\PYG{l+s}{exl q(RealLt 0 (g xx) q)) \PYGZhy{}\PYGZgt{} PosDef g)\PYGZdq{}}
\PYG{l+s}{\PYGZdq{}PosDefIntro\PYGZdq{}}\PYG{p}{))}
\end{Verbatim}

While doing theorem proving, we can ease the process of deciding whether a given function is positive definite by assuming the existence of a witness $\eta$ which actually derives the positive definiteness of function if all inputs of $g$ are considered to be rational vectors.


\begin{defn}[\texttt{PosDefRatWit}]\label{def:PosDefRatWit} A uniformly continuous function $g:B=(c,R) \rightarrow \mathbb{R}$ given by $(h_g,\alpha_{g},\omega_{g},(c,R))$ is called positive definite with rational witness $\eta: \mathbf{P} \rightarrow \mathbf{P}$ whenever 
	\begin{enumerate}
		\item $\|\bold c \| \le R$,
		\item $\bold e=0 \rightarrow g(\bold e) = 0$,
		\item $\forall \bold e,p \, ( 0 \#_p \spc \bold e \rightarrow  0 <_{\eta(p)} \spc g(\bold e))$,  
	\end{enumerate}
\end{defn}
and the corresponding IDPC in \texttt{Minlog}:
\if{
\begin{minted}[linenos,frame=single,framesep=5pt]{scheme}
(add-ids
(list (list "PosDefRatWit" (make-arity (py "contmv") 
(py "(pos=>pos)"))))
'("all g,eta(ContMV g -> 
(|_ g defC <<= g defB
& 0===g (ZerosRealVector (g defC dim)))
& all e( e dim = g dim -> all p (RealLt 0 (||_ e) p ->  
(RealLt 0 (g e) (eta p)))) -> PosDefRatWit g eta)" 
"PosDefRatWitIntro"))
\end{minted}
}\fi
\begin{Verbatim}[commandchars=\\\{\}]
\PYG{p}{(}\PYG{n+nf}{add\PYGZhy{}ids}
\PYG{p}{(}\PYG{n+nb}{list }\PYG{p}{(}\PYG{n+nb}{list }\PYG{l+s}{\PYGZdq{}PosDefRatWit\PYGZdq{}} \PYG{p}{(}\PYG{n+nf}{make\PYGZhy{}arity} \PYG{p}{(}\PYG{n+nf}{py} \PYG{l+s}{\PYGZdq{}contmv\PYGZdq{}}\PYG{p}{)}
\PYG{p}{(}\PYG{n+nf}{py} \PYG{l+s}{\PYGZdq{}(pos=\PYGZgt{}pos)\PYGZdq{}}\PYG{p}{))))}
\PYG{o}{\PYGZsq{}}\PYG{p}{(}\PYG{l+s}{\PYGZdq{}all g,eta(ContMV g \PYGZhy{}\PYGZgt{}}
\PYG{l+s}{(|\PYGZus{} g defC \PYGZlt{}\PYGZlt{}= g defB}
\PYG{l+s}{\PYGZam{} 0===g (ZerosRealVector (g defC dim)))}
\PYG{l+s}{\PYGZam{} all e( e dim = g dim \PYGZhy{}\PYGZgt{} all p (RealLt 0 (||\PYGZus{} e) p \PYGZhy{}\PYGZgt{}}
\PYG{l+s}{(RealLt 0 (g e) (eta p)))) \PYGZhy{}\PYGZgt{} PosDefRatWit g eta)\PYGZdq{}}
\PYG{l+s}{\PYGZdq{}PosDefRatWitIntro\PYGZdq{}}\PYG{p}{))}
\end{Verbatim}

\begin{thm}\label{thm:PosDef}
	$\forall g,\eta \, (\texttt{PosDefRatWit} \spc g \spc \eta \rightarrow \texttt{PosDef} \spc g)$.
\end{thm}
\begin{proof}
	Without loss of generality consider $g$ as a uniformly continuous function of two variables. The first two properties are clear. To show the third property one may assume that for all $x$ there exists $p$ such that $0 \#_p \spc \bold x$. It remains to show that $\exists \spc q \spc 0 <_q \spc g(\bold x))$. We claim that $q=\eta(p)+2$.
The goal $0 <_{\eta(p)+2} g(\bold x)$ is equivalent to $g(\bold x) \in_{\eta(p)+2} \mathbf{R}^+$. By applying  Theorem \ref{RealPosChar}  (\texttt{RealPosChar}) backwards with the parameter $\alpha_{g}(\eta(p)+1) $,  the new goal is
\[
\forall n (\alpha_{g}(\eta(p)+1) \le n \rightarrow \frac{1}{2^{\eta(p)+1}} \le h_{g(e_1(n),e_2(n),n)}).
\]
Now, for any $n$ such that $\alpha_{g}(\eta(p)+1) \le n$, the intermediate step is
\begin{align*}
\frac{1}{2^{\eta(p)}} - \frac{1}{2^{\eta(p)+1} } \le & h_{g}(e_1(n),e_2(n),n) + h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4)) \\
& - h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4)) \,.
\end{align*}
Observe that
\begin{align*}
\frac{1}{2^{\eta(p)}}  \le h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4))
\end{align*}
holds by assumption after unfolding  Definition \ref{def:PosDefRatWit} \ $\texttt{PosDefRatWit} \spc g \spc \eta$.
In addition,
\begin{align*}
 - \frac{1}{2^{\eta(p)+1} } \le h_{g}(e_1(n),e_2(n),n) - h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4))
\end{align*}
holds since
\begin{align*}
  &h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4)) - h_{g}(e_1(n),e_2(n),n) \le\\
  &|h_{g}(e_1(n),e_2(n),\alpha_{g}(\eta(p)+4)) - h_{g}(e_1(n),e_2(n),n)| \le  \frac{1}{2^{\eta(p)+1} } \,,
\end{align*}
by applying Theorem \ref{thm:CauchyElim} (\texttt{CauchyElim}) with the parameter $\alpha_{g}$ and by the fact that $\frac{1}{2^{p}} \le |e_1(n)| + |e_2(n)|$ whenever $\alpha_{g}(\eta(p)+1) \le n$.
\end{proof}

The next section focuses on examples where we certify formally the stability of continuous/discrete-time dynamical systems.
\section{Examples for verification of Lyapunov functions}
\label{sec:benchs}
Here we combine the previous formalization framework from Section \ref{sec:formal} with the hybrid numeric/symbolic algorithm \texttt{exactLyapunov}, stated in Algorithm  \ref{alg:exactLyapunov} from Section~\ref{sec:prelim_exact}. 
Namely, for a given vector $f$ of polynomials, we first compute a weighted rational SOS  function $V$ with \texttt{exactLyapunov} and then formally prove that $V$ is a positive definite form. 
If we succeed in proving that $-\grad V \cdot f$ is a nonnegative (resp. positive definite) function then the system is formally proved to be stable (resp. asymptotically stable).\\

{\itshape{Example 1:}}
We start with the following two-dimensional example of continuous-time polynomial system:
\begin{align}
\label{ex:1}
\dot{x}_1 & = -x_1^3 + x_2 \,, \nonumber \\
\dot{x}_2 & = - x_1 - x_2 \,.
\end{align}
A quadratic SOS Lyapunov function $V_2(x_1,x_2) = 1.115 (x_1^2+x_2^2)$ is obtained via \texttt{exactLyapunov}.

Then, one can formally verify within $\minlog$ via Theorem \ref{thm:PosDef} that the following (in)-equalities hold:
\begin{align*}
V_2(x_1,x_2) - (x_1^2+x_2^2) & =  0.115 (x_1^2+x_2^2) \geq 0 \,, \\
- \grad V_2 \cdot f & = 223/100 (x_1^4 + x_2^2) \,,
\end{align*}

so that the proof system can assert that  $V_2$ and $-\grad V_2 \cdot f$ are \texttt{PosDef}. In order to do that it suffices  to show that the function $W(x_1,x_2)=x_1^2+x_2^2$ is \texttt{PosDef}. Without loss of generality assume that the function $W$ is defined on the unit ball $B(0,1)$ together with the following functions:
\begin{itemize}
	\item $h_{W}(e_1,e_2)=e_1^2+e_2^2$,
	\item $\alpha(p)=0$,
	\item $\omega(p)=p+3$.
\end{itemize}
To obtain the map $\eta(p)$, required to prove  \texttt{PosDefRatWit}, we proceed by means of  program extraction. That is the first goal of the proof is:
\begin{align}
	\forall e(\exists p \spc 0 \#_p \bold e \rightarrow \exists q \spc 0<_q W(\bold e)) \,.
\end{align}
Observe that the assumption $0 \#_p \bold e$ is equivalent to the fact that $\exists p \in \mathbf{P} \spc 2^{-p} \le |e_1|+|e_2|$ which is equivalent to $0<|e_1|+|e_2|$ which is in turn equivalent to $0<|e_1| \vee 0 < |e_2|$. The last two equivalences follow from the fact that $e_1,e_2 \in \mathbf{Q}$. Additionally, the normalization-by-evaluation converts the goal $\exists q \spc 0<_q W(\bold e)$  to the goal $\exists q \spc 2^{-q} \le e_1 \cdot e_1 + e_2 \cdot e_2$,  also due to the fact that $e_1,e_2 \in \mathbf{Q}$ and $\alpha(p)=0$.
We may now proceed by case distinction:\\
\noindent\hspace*{5mm}{\itshape{Case 1:}} $0<|e_1|$\\
\noindent\hspace*{11mm} We have $p'_1 $ such that $2^{-p'_1} \le |e_1|$ and set $q=p'_1+p'_1$.\\
\noindent\hspace*{5mm}{\itshape{Case 2:}} $0 \nless e_1 \wedge 0<|e_2|$\\
\noindent\hspace*{11mm}Similarly to {\itshape{Case 1:}} we have $p'_2$  such that $2^{-p'_2} \le |e_1|$ and set $q=p'_2+p'_2$.\\
\noindent\hspace*{5mm}{\itshape{Case 3:}} $0 \nless e_1 \wedge 0 \nless |e_2|$\\
\noindent\hspace*{11mm}The goal follows from {\itshape{ex-falso-quodlibet}} since $0<|e_1| \vee 0 < |e_2|$.\\
Now let $\eta: \mathbf{P} \rightarrow \mathbf{P}$ be a program which generates $q$ for any $p$ such that $0 \#_p \bold e$. In \texttt{Minlog} the extracted $\eta$ is represented by: 
\if{
\begin{minted}[linenos,frame=single,framesep=5pt]{scheme}
[e0,p1]
[if e0 ([rls2] [if rls2 1 ([a3,rls4] [if rls4 1 ([a5,rls6] 
[if rls6 [if (a3*a3+a5*a5) ([k7,p8] [if k7
([p9](Rec pos=>pos)p8 1 ([p10]PosS)([p10,p11]PosS(PosS p11)))
1([p9]1)])]([a7,rls8]1)])])])]
\end{minted}
}\fi
\begin{Verbatim}[commandchars=\\\{\}]
\PYG{p}{[}\PYG{n+nv}{e0,p1}\PYG{p}{]}
\PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{e0} \PYG{p}{([}\PYG{n+nv}{rls2}\PYG{p}{]} \PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{rls2} \PYG{l+m+mi}{1} \PYG{p}{([}\PYG{n+nv}{a3,rls4}\PYG{p}{]} \PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{rls4} \PYG{l+m+mi}{1} \PYG{p}{([}\PYG{n+nv}{a5,rls6}\PYG{p}{]}
\PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{rls6} \PYG{p}{[}\PYG{k}{if }\PYG{p}{(}\PYG{n+nf}{a3*a3+a5*a5}\PYG{p}{)} \PYG{p}{([}\PYG{n+nv}{k7,p8}\PYG{p}{]} \PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{k7}
\PYG{p}{([}\PYG{n+nv}{p9}\PYG{p}{](}\PYG{n+nf}{Rec} \PYG{n+nv}{pos=\PYGZgt{}pos}\PYG{p}{)}\PYG{n+nv}{p8} \PYG{l+m+mi}{1} \PYG{p}{([}\PYG{n+nv}{p10}\PYG{p}{]}\PYG{n+nv}{PosS}\PYG{p}{)([}\PYG{n+nv}{p10,p11}\PYG{p}{]}\PYG{n+nv}{PosS}\PYG{p}{(}\PYG{n+nf}{PosS} \PYG{n+nv}{p11}\PYG{p}{)))}
\PYG{l+m+mi}{1}\PYG{p}{([}\PYG{n+nv}{p9}\PYG{p}{]}\PYG{l+m+mi}{1}\PYG{p}{)])]([}\PYG{n+nv}{a7,rls8}\PYG{p}{]}\PYG{l+m+mi}{1}\PYG{p}{)])])])]}
\end{Verbatim}
which takes as input a rational vector $\bold e_0$ and a positive number $p1$. If $0\#_{p_1} \bold e_0$, it generates a positive number $\eta(p_1)$ such that $0<_{\eta(p_1)} f(\bold e_0)$.
Trivially, \texttt{PosDefRatWit} $W$ $\eta$ follows from the derivation above. By Theorem  \ref{thm:PosDef}, 
the function $W$ is also \texttt{PosDef}. 
Similarly, ignoring the factor $223/100$, for $\grad V_2 \cdot f$, we obtain the following extracted $\eta(p)$:
\if{
\begin{minted}[linenos,frame=single,framesep=5pt]{scheme}
[e0,p1][if e0 ([rls2][if rls2 1 
([a3,rls4][if rls4 1 ([a5,rls6][if rls6 
[if (a3*a3*a3*a3+a5*a5) ([k7,p8]
[if k7 ([p9](Rec pos=>pos)p8 1([p10]PosS)([p10,p11]PosS(PosS p11)))
 1 ([p9]1)])] ([a7,rls8]1)])])])]
\end{minted}
}\fi
\begin{Verbatim}[commandchars=\\\{\}]
\PYG{p}{[}\PYG{n+nv}{e0,p1}\PYG{p}{][}\PYG{k}{if }\PYG{n+nv}{e0} \PYG{p}{([}\PYG{n+nv}{rls2}\PYG{p}{][}\PYG{k}{if }\PYG{n+nv}{rls2} \PYG{l+m+mi}{1}
\PYG{p}{([}\PYG{n+nv}{a3,rls4}\PYG{p}{][}\PYG{k}{if }\PYG{n+nv}{rls4} \PYG{l+m+mi}{1} \PYG{p}{([}\PYG{n+nv}{a5,rls6}\PYG{p}{][}\PYG{k}{if }\PYG{n+nv}{rls6}
\PYG{p}{[}\PYG{k}{if }\PYG{p}{(}\PYG{n+nf}{a3*a3*a3*a3+a5*a5}\PYG{p}{)} \PYG{p}{([}\PYG{n+nv}{k7,p8}\PYG{p}{]}
\PYG{p}{[}\PYG{k}{if }\PYG{n+nv}{k7} \PYG{p}{([}\PYG{n+nv}{p9}\PYG{p}{](}\PYG{n+nf}{Rec} \PYG{n+nv}{pos=\PYGZgt{}pos}\PYG{p}{)}\PYG{n+nv}{p8} \PYG{l+m+mi}{1}\PYG{p}{([}\PYG{n+nv}{p10}\PYG{p}{]}\PYG{n+nv}{PosS}\PYG{p}{)([}\PYG{n+nv}{p10,p11}\PYG{p}{]}\PYG{n+nv}{PosS}\PYG{p}{(}\PYG{n+nf}{PosS} \PYG{n+nv}{p11}\PYG{p}{)))}
 \PYG{l+m+mi}{1} \PYG{p}{([}\PYG{n+nv}{p9}\PYG{p}{]}\PYG{l+m+mi}{1}\PYG{p}{)])]} \PYG{p}{([}\PYG{n+nv}{a7,rls8}\PYG{p}{]}\PYG{l+m+mi}{1}\PYG{p}{)])])])]}
\end{Verbatim}

Overall, we have formally proved that the system is asymptotically stable.\\
{\itshape{Example 2:}}
Let us consider another continuous-time polynomial system from \cite[\S~4.2]{sostools}:
\begin{align}
\label{ex:2}
\dot{x}_1 & = - x_1^3 - x_1 x_3^2 \,, \nonumber \\
\dot{x}_2 & = - x_2 - x_1^2 x_2 \,, \\
\dot{x}_3 & = - x_3 - \frac{3 x_3}{x_3^2 + 1} x_1 + 3 x_1^2 x_3  \,. \nonumber
\end{align}
A homogeneous quadratic SOS Lyapunov $V_3(x_1,x_2) =  5.5489 x_1^2+ 4.1068 x_2^2 + 1.7945 x_3^2$ is computed with \texttt{exactLyapunov}.
Then, one formally verifies within $\minlog$ that the following (in)-equalities hold:
\begin{align*}
V_3(x_1,x_2,x_3) - (x_1^2+x_2^2+x_3^2) & =  4.5489 x_1^2+ 3.1068 x_2^2 + 0.7945 x_3^2  \geq 0 \,, \\
-\grad V_3 \cdot (x_3^2+1) f & = 55489/5000 \big(569705563  
(2 x_2  /55489)^2 \\
& + 
995750105 ( 2/55489  x_3  - 495266 / 24893752625 x_1^2 x_3   )^2 \\
& + (x_1^2 + \frac{ 119248 x_3^2}{1387225 })^2 +
569705563 (2 x_1 x_2 / 55489)^2 \\
& + 98987825902 (x_1 x_3  /277445)^2 +
569705563 (2 x_2 x_3  / 55489)^2 \\
& +
38212870902320486809 (x_3^2 /10996532575)^2 \\
& + 
375452958989124261245 (x_1^2  x_3 /24893752625)^2 \\ 
& +
569705563 (2 x_1 x_2 x_3 / 55489)^2 +
91778806 (x_1 x_3^2  /55489)^2 
\big) \geq 0
\,,
\end{align*}

so that the proof system can assert that  $V_3$ is \texttt{PosDef} 
and that $-\grad V_3 \cdot f$ is \texttt{NonNeg}. 
 For this later goal, Theorem \ref{RealNNegChar} (RealNNegChar) is applied backwards and the resulting goal is read as:
 \begin{align*}
 \forall p \spc \exists n_0 \spc \forall n \ge n_0 \spc -2^{p} \le  & \spc 55489/5000 \big(569705563  
 (2 e_2(n) /55489)^2 \\
 & + 
 995750105 ( 2/55489  e_3(n)  - 495266 / 24893752625 e_1(n)^2 e_3(n)   )^2 \\
 & + (e_1(n)^2 + \frac{ 119248 e_3(n)^2}{1387225 })^2 +
 569705563 (2 e_1(n) e_2(n) / 55489)^2 \\
 & + 98987825902 (e_1(n) e_3(n)  /277445)^2 +
 569705563 (2 e_2(n) e_3(n)  / 55489)^2 \\
 & +
 38212870902320486809 (e_3(n)^2 /10996532575)^2 \\
 & + 
 375452958989124261245 (e_1(n)^2  x_3 /24893752625)^2 \\ 
 & +
 569705563 (2 e_1(n) e_2(n) e_3(n) / 55489)^2 +
 91778806 (e_1(n) e_3(n)^2  /55489)^2 
 \big) \,.
 \end{align*}
 For any $p$, after introducing $n_0=0$,  the goal is trivial by proving the nonnegativity of each summand separately.
Overall, we have formally proved  that the system is stable.\\
{\itshape{Example 3:}} Eventually, we take the following  two-dimensional discrete-time rational system:
\begin{align}
\label{ex:3}
x_{t+1} & = \frac{y_t}{1+x_t^2} \,, \nonumber \\
y_{t+1} & = \frac{x_t}{1+y_t^2} \,.
\end{align}
In a similar way as for continuous-time polynomial systems, one can compute a homogeneous quadratic SOS Lyapunov $V_2(x,y) = x^2+y^2$ with \texttt{exactLyapunov}, then formally verify within $\minlog$ that the following (in)-equalities hold:

\begin{align*}
V_2(x,y)  & =   x^2+y^2 \succ 0 \,, \\
(1+x^2)^2(1+y^2)^2 (V_2(x,y) - V_2 (f(x,y)))  & = x^6 y^4 + x^4 y^6 + 2 x^6 y^2   \\
& + 4 x^4 y^4 + 2 x^2 y^6 + 5 x^4 y^2  \\
& + 5 x^2 y^4 + 4 x^2 y^2   \geq 0 \,.
\end{align*}

\if{
\begin{align*}
V_2(x,y)  & =   x^2+y^2 \succ 0 \\
(1+x^2)^2(1+y^2)^2 (V_2(x,y) - V_2 (f(x,y)))  & = (2 x y - 2/3 x^3 y - 2/3 x y^3)^2 \\
& + 69 (1/2 x^2 y - 9/161 x^2 y^3 )^2 \\
& + 69 (1/3 x y^2 - 9/161 x^3 y^2)^2 \\
& + 14 (1/3 x^3 y - 155/588 x y^3)^2 \\
& + 75/7 (x^2 y^2)^2 + 1599/2744 (x y^3)^2 \\
& + 884/1127 (x^3 y^2)^2 +  884/1127 (x^2 y^3)^2 
 \geq 0 \,,
\end{align*}
}\fi

Doing so, the proof system can assert that $V_2$ is \texttt{PosDef} and $V_2(x,y)-V_2(f(x,y))$ is \texttt{NonNeg}.
For the goal \texttt{NonNeg} $V_2(x,y)-V_2(f(x,y))$ one relies on Theorem \ref{RealNNegChar} (RealNNegChar) backwards, again for the sake of easing the derivation and  obtain the goal 
\begin{align*}
\forall p \spc \exists n_0 \spc \forall n \ge n_0 \spc -2^{p} \le  & \spc e_1(n)^6 e_2(n)^4 + e_1(n)^4 e_2(n)^6 + 2 x^6 e_2(n)^2  \\
& + 4 e_1(n)^4 e_2(n)^4 + 2 e_1(n)^2 e_2(n)^6 + 5 e_1(n)^4 e_2(n)^2  \\
& + 5 e_1(n)^2 e_2(n)^4 + 4 e_1(n)^2 e_2(n)^2 \,.
\end{align*}
As for the second example, the goal is easy to show by proving the nonnegativity of each summand separately and yields a formal stability proof for this last system.


\if{
Let us consider the continuous-time rational system from~\cite[Example~2.1]{ahmadi2011converse}.
\begin{align}
\label{ex:2}
\dot{x_1} & = -0.15 x_1^7 + 200 x^6_1 x_2 - 10.5x ^5_1 x^2_2 - 807 x^4_1 x^3_2 + 14 x^3_1 x^4_2 \nonumber \\
& \quad + 600 x^2_1 x^5_2 - 3.5 x^1 x^6_2 + 9 x^7_2 \nonumber \\
\dot{x_2} & = - 9 x^7_1 - 3.5 x^6_1 x_2 - 600 x^5_1 x^2_2 + 14 x^4_1 x^3_2 + 807 x^3_1 x^4_2 - 10.5 x^2_1 x^5_2\\
& \quad - 200 x_1 x^6_2 - 0.15 x^7_2 \nonumber \,.
\end{align}

An SOS Lyapunov function $V_8$ of degree 8 for~\eqref{ex:2} can be computed via semidefinite programming and reads as follows:
\begin{align*}
V_8(x_1,x_2) & =  0.02 x^8_1 + 0.015 x^7_1 x_2  + 1.743 x^6_1 x^2_2 - 0.106 x^5_1 x^3_2 -3.517 x^4_1 x^4_2 + 0.106 x^3_1 x^5_2 \\
& \quad + 1.743 x^2_1 x^6_2 - 0.015 x_1 x^7_2 + 0.02 x^8_2 
\,.
\end{align*}
}\fi

\section{Conclusion}
We provide  a fully computer-assisted framework to formally prove the stability of a given continuous/discrete-time polynomial system. 
For this, our framework consists of proving formally that a given function is a Lyapunov function, by means of exact SOS certificates.

As a proof of concept, we illustrate our formalization framework on simple examples with up to three variables. 

In order to apply such techniques to higher dimensional systems, one needs to overcome the typical scalability issues related to SOS-based optimization on polynomial systems. 
One possible remedy is to exploit the sparsity of the input data, either when there are few correlations between the variables \cite{tacchi2019approximating} or when there are few involved  monomials terms \cite{tssos,chordal-tssos,cstssos}.
An alternative to SOS is to rely on other nonnegativity certificates, such as sums of nonnegative circuits or sums of arithmetic-geometric-exponentials, which can be computed with more efficient classes of conic programming, such as second-order cone programming \cite{wang2019second}. 
Exact certificates can also be obtained in this case \cite{magron2019exact}.
An interesting track of further research is to derive formally certified outer/inner approximations of sets of interest arising in the context of dynamical systems, such as backward/forward reachable sets or maximal invariants.

\paragraph{\textbf{Acknowledgements}.} 
The second author was supported by the FMJH Program PGMO (EPICS project) and  EDF, Thales, Orange et Criteo, as well as from the Tremplin ERC Stg Grant ANR-18-ERC2-0004-01 (T-COPS project).
This work has benefited from the Tremplin ERC Stg Grant ANR-18-ERC2-0004-01 (T-COPS project), the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Actions, grant agreement 813211 (POEMA) as well as from the AI Interdisciplinary Institute ANITI funding, through the French ``Investing for the Future PIA3'' program under the Grant agreement n$^{\circ}$ANR-19-PI3A-0004.

\bibliographystyle{plain}
\input{FormalSOSLF.bbl}

\end{document}
{"config": "arxiv", "file": "2006.09884/FormalSOSLFarxiv.tex"}