DeepAI

# Computer-assisted proofs for Lyapunov stability via Sums of Squares certificates and Constructive Analysis

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.

• 3 publications
• 17 publications
• 5 publications
12/05/2021

### Windmills of the minds: an algorithm for Fermat's Two Squares Theorem

The two squares theorem of Fermat is a gem in number theory, with a spec...
03/11/2019

### How far away must forced letters be so that squares are still avoidable?

We describe a new non-constructive technique to show that squares are av...
11/11/2017

### Real-number Computability from the Perspective of Computer Assisted Proofs in Analysis

We present an interval approach to real-number computations. In some asp...
05/04/2022

### Constructive Analysis in the Agda Proof Assistant

Proof assistant software has recently been used to verify proofs of majo...
04/06/2021

### Extraction of a computer-certified ODE solver

Reliably determining system trajectories is essential in many analysis a...
03/12/2021

### Packing Squares into a Disk with Optimal Worst-Case Density

We provide a tight result for a fundamental problem arising from packing...
02/16/2022

### A Survey of Semen Quality Evaluation in Microscopic Videos Using Computer Assisted Sperm Analysis

The Computer Assisted Sperm Analysis (CASA) plays a crucial role in male...

## 1. Introduction

Computer-assisted theorem proving – performed by special software called proof assistants – have been successfully used for important problems, including the formalization of Kepler’s conjecture [28]

and a proof of the odd order theorem in group theory

[26]. 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 second method of Lyapunov, also called the 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 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 [25], we are interested in reasoning with the so-called exact real arithmetic. For this purpose the standard library of 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.

### Related works

In the context of dynamical systems, [12] formally proves the correctness of a program implementing the numerical resolution of a wave equation, with Frama-C [19]. In [31]

, the authors formalize variational equations related to the solutions of ordinary differential equations in Isabelle/HOL

[50]. Beyond formalization of numerical analysis procedures, several frameworks are currently available to perform computer-assisted proofs of mathematical analysis. The work in [1] proposes formalized mathematical structures in Coq [18] for function spaces. Earlier developments have been pursued towards the formalization of real analysis in several proof assistants, including Coq [46], HOL-Light, [30], Isabelle/HOL [21] or ACL2 [22], as well as constructive real analysis [24, 35]. Proof assistants also play an important role in formal verification of control systems. Zou [78] used Isabelle for implementation of their control system formal verification based on Hybrid Hoare logic. The Why3 platform [11] coupled with MATLAB/Simulink was used in [3, 4] to perform simple stability checks of linear discrete systems with quadratic Lyapunov functions. Gao [23] mentioned several proof assistant software tools to be considered for implementation of formal verification of control systems, including Coq [18], HOL-Light [29], Isabelle/HOL and Lean [20]. The most closely related work is due to [2] where the authors present the framework for the certified programs for robots based on constructive real numbers [52]. 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 [61] the author formally proves the soundness of a control function for the inverted pendulum within the classical logic. The author in [61] also mentions the impossibility of generating and running programs with that framework.

In the present work the interactive proof system Minlog [64, 7] 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 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

[7, 8]. The conversion of proofs to instances handled by the soundness theorem can be done automatically [7, 47, 5]. This kind of approach is not used in, e.g., Coq [18], 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 [7].

Consider the autonomous continuous-time dynamical system

 (1.1) ˙x=f(x)

where is a locally Lipschitz function with equilibrium point at origin and . A related result [33] states that the equilibrium of the system (1.1) is asymptotically stable if there is some positive definite homogeneous function such that the derivative of 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 is called a Lyapunov function. For polynomial systems, i.e., when is a polynomial, these two sufficient positivity/nonnegativity conditions can be strengthened by computing a positive definite homogeneous function which is a sum of squares (SOS) of polynomials and such that is SOS [53]. In this case, we refer to as an 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

[37, 53]. In [37], 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) [70]

, 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

[27], control theory [13]

or parameter estimation

[67]. Available software libraries such as SeDuMi [68], SDPA [76] or Mosek [48] 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 moment-SOS or Lasserre’s hierarchy. The primal is a 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

[9]

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

. 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 [57]. For dynamical systems, this hierarchy have been extended multiple times, for example to obtain converging approximations for optimal control problems [38], reachable sets [43] or invariant measures [40].

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 [16, 42] or in the multivariate case [56, 32, 44]. In particular, the rounding-projection algorithm from [56] outputs a weighted rational SOS decompositon of polynomials lying in the interior of the SOS cone. More recently, [44] 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 [41].

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 [16] or formal verification of real inequalities [39] within proof assistants such as Coq [18] or HOL-Light [29].

### 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 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 [41] tool available outside of the proof assistant MinLog and verified inside. Our first contribution, presented in Section 2.3, is to provide an algorithm, called exactLyapunov to compute exact weighted rational SOS Lyapunov functions for a given continuous-time polynomial system (1.1). 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 3 is dedicated to the presentation of the Minlog formalization framework and its extension of the standard library. Our further contribution, provided in Section 4, is to derive a full framework to formally prove the stability of such systems, by verifying the output of exactLyapunov within the Minlog system. We illustrate our formal framework with several benchmarks in Section 5.

## 2. Preliminary Background and Results

### 2.1. Positive polynomials, sums of squares

Let (resp. ) be the field of rational (resp. real) numbers. With , let us denote by (resp. ) the space of real-valued -variate polynomials (resp. of degree at most ) in the variable . Given a real symmetric matrix , the notation (resp. ) means that is positive semidefinite (resp. positive definite), that is has only nonnegative (resp. positive) eigenvalues. Let us note the set of sums of squares (SOS) of polynomials and the set of SOS polynomials of degree at most . An -variate polynomial belongs to if there exist a nonzero and such that . An example of SOS polynomial is . For each multi-index , we use the notation , and we write . The support of a polynomial written as above is the set . For our previous example, the support is . We say that has full support if its support is .

It is trivial that the condition of to be SOS ensures that is nonnegative over . 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 [58] and the references therein for more details. The condition that is equivalent to the existence of a symmetric matrix such that , where

is the vector of all monomials of degree at most

. The length of this vector is equal to . The matrix is called a Gram matrix [17] associated to and computing boils down to solving a semidefinite program (SDP for short) [70]. We refer to this technique as the Gram matrix method. The SOS decomposition is a by-product of a decomposition of the form , with being a lower triangular matrix with nonnegative real entries on the diagonal and being a diagonal matrix with nonnegative real entries, such that .

### 2.2. SOS Lyapunov functions

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 , let us now consider a continuous polynomial system (1.1). Let be the maximal degree of the polynomials . The Gram matrix method can be used to certify that such a system with equilibrium at the origin is stable by computing a Lyapunov polynomial function . More precisely, one would like to obtain such that (1) is a positive definite form and (2) is nonnegative.

Regarding (1), since we expect to be a homogeneous polynomial, we look for a decomposition , where is the vector of all monomials of degree equal to . Regarding (2), we assume that is odd, so that the degree of is even, i.e., that is a natural integer (otherwise 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 , and such that:

 (SDPk)

The optimization variables of the above feasibility program (SDP) are the coefficients of and the entries of the matrices and . Any feasible solution of this SDP ensures that is a positive definite form and is nonnegative, yielding a stability proof for the initial system (1.1).

To solve (SDP), one relies on numerical solvers implemented in double precision, for instance SeDuMi [68] and SDPA [76], or arbitrary-precision solvers such as SDPA-GMP [49]. However, such numerical solvers provide approximate nonnegativity certificates, so that the polynomial , and the matrices do not exactly satisfy the equality constraints of (SDP).

### 2.3. From approximate to exact SOS certificates

To overcome these feasibility issues, one can rely on the certification framework from [44], implemented in the RealCertify software library [41]. For certification purpose, we are mainly interested in polynomials with rational coefficients, which belong to the interior of the SOS cone , denoted by . Given a polynomial of degree with full support, there always exists by [17, Proposition 5.5] a positive definite Gram matrix associated to , i.e., such that . In this case, one can apply the so-called intsos algorithm, stated in [44, Algorithm 1], to compute exact weighted rational SOS decompositions for polynomials with rational coefficients, lying in . Namely, we obtain as output a decomposition , with being nonnegative rational numbers and being polynomials with rational coefficients.

For the sake of self-containedness, let us briefly recall how intsos works. The algorithm takes as input a polynomial of degree with rational coefficients and full support, and accuracy parameters and . Then, the procedure consists of a perturbation step and an absorption step. The perturbation step consists of finding a small enough positive rational number , such that the perturbed polynomial also lies in . Then, intsos calls an external SDP solver with accuracy to find an approximate Gram matrix associated to , such that . Here the size (and rank) of is . Then, we obtain an approximate Cholesky’s decomposition of with accuracy , and for all , we compute the inner product of the -th row of by to obtain a polynomial . For large enough accuracy of the SDP solver and accuracy of the Cholesky’s decomposition, the coefficients of the remainder become small enough to be “absorbed”: that is admits itself an SOS decomposition, and so does . If 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 exactLyapunov, similar to intsos, to compute exact weighted rational SOS Lyapunov functions for a given continuous-time polynomial system (1.1).

###### Proposition 2.1.

Let be a vector of rational polynomials of degree at most , with being odd. Assume that (SDP) admits a strictly feasible solution for some . Then exactLyapunov terminates and returns an exact rational weighted SOS Lyapunov function for the continuous time system (1.1).

###### Proof.

The loop going from line line:ki to line line:kf increases the search degree for the SOS Lyapunov function. By assumption, (SDP) admits a strictly feasible solution for some large enough , so this loop eventually terminates. In particular, we obtain an SOS Lyapunov function of degree satisfying the constraints of (SDP). Then, one applies the perturbation-absorption procedure intsos to obtain an exact rational weighted SOS decomposition of at line line:V, and of at line line:gradV. By [44, Proposition 3.5], the procedure intsos always terminates in both cases, since and have associated positive definite Gram matrices. The outputs and respectively contain a list of nonnegative rational numbers and rational polynomials such that . Similarly, we obtain an exact rational weighted SOS decomposition for . ∎

## 3. Minlog system

In this section, we investigate how to formally prove stability of polynomial systems, by verifying the output of exactLyapunov within Minlog. We first introduce some preliminary concepts of the proof system 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 [14, 24, 21, 30, 35].

### 3.1. 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 [14]. Let and be some arbitrary formulas. The most important BHK rules are:

• Proving the implication means finding an algorithm that transforms a proof of into a proof of .

• To prove , an algorithm computing and the proof that holds are required.

• To prove , an algorithm is required that converts into a proof of .

In practice, proof assistant software are concrete realizations of the BHK interpretation.

The interactive proof system 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 [65]. 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 “”, where and 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 constructs a new type from the previously defined base types and , e.g., the successor of a natural number of type nat is viewed as a function of type . 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 nat of natural numbers by stating that Zero is of type nat and Succ is a mapping of type natnat. The normalization of the term 1+2 would result in Succ(Succ(Succ Zero)) which is clearly of type nat . The following types for the numbers with constructors and destructors are contained in the standard Minlog library:

• pos:One,SZero,SOne (positive numbers (binary))

• nat:Zero,Succ (natural numbers )

• list:Nil,Cons

• int:IntP,IntZero,IntN (integer numbers )

• rat:RatConstr and destructors RatN,RatD (rational numbers )

• rea:RealConstr and destructors
RealSeq,RealMod (real numbers )

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 NatTimes is of type natnatnat and the computational rules are given recursively as:

 NatTimesnZero→Zero NatTimesnSuccm→% NatTimesnm+n.

As an additional example consider the strict inequality relation again for the natural numbers. The program constant NatLt is of type natnatbool and the computational rules are:

 NatLtnZero→False NatLt0Succn→True NatLtSuccnSuccm→NatLtnm.

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 . Then one can add the machine-generated computation rule to the library of rewriting rules. Every time the normalization process is done, if some rational occurs then the zero is canceled automatically.

### 3.2. The notion of totality

In order to establish the connection between proofs and programs, each program constant in 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 [66, Chapter 8.3] and [6]. In order to use the whole machinery of the computational model in 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 , the program constant NatPlus n m is recursively defined as:

• NatPlus n Zero

• NatPlus n Succ m Succ(NatPlus n m).

The totality goal is: . By induction on m and by definition is equivalent to TotalNat(n) which follows from the assumptions. The induction step is TotalNat(n +Succ m) which is equivalent by definition to TotalNat(Succ(n +m). Finally the proof tree is split into two goals TotalNat(Succ(n) and 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.

### 3.3. 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 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 RatInUnitI is given by:

 ∀a∈Q0≤a→a≤1→% RatInUnitIa

which means semantically that to show one has to show that and . Additionally the definitions can also be provided inductively. As an example consider the definition of an even natural number given by the IDPC Even. Inductively one can proceed as follows:

 Even Zero ∀n(Evenn→Evenn+2).

Here the definition of the even natural number is given in terms of other natural number (in this case Zero) which is known to be even.

### 3.4. Preliminaries for the formalization in constructive setting

This section is based on the foundation of constructive analysis with witnesses [62]. 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 would be a map such that for any and , one has whenever . By relying on the axiom of choice [77, p. 89] or the dependent choice axiom [63], a continuous inverse function can be constructed in terms of the witness .

All presented concepts are implemented within of Minlog, and furthermore we extended the standard library of 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 Minlog the reals are build by introducing type variables as: , M:, defining the constructor by the arrow type formation RealConstr:asM and introducing the IDPC Real, i.e., a real is given by a regular Cauchy sequence of rationals (encoded by the type variable as), denoted as Cauchy, with a given monotone modulus , denoted as Mon, that is whenever and . 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 , proposition holds”. The next theorem is an important example for such elimination theorems:

###### Theorem 3.1 (CauchyElim).
 ∀(an)n,M{Cauchy}(an)nM→∀p,m,n|an−am|≤2−p.

A real is called nonnegative if . For two real numbers , denote whenever the real is nonnegative. For a real denote (alternately ) if . Two reals and are equal, if . Note that the equality of the reals is undecidable. Two reals , are called apart if which is usually denoted by . The next two characterization theorems are often useful for operating with inequalities within Minlog.

###### Theorem 3.2 (RealNNegChar).

For a real the following are equivalent:

1. .

###### Theorem 3.3 (RealPosChar).

For a real the following are equivalent:

1. .

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.

###### Definition 3.4 (Cont).

A uniformly continuous function on a compact interval with rational endpoints is given by

1. an approximation map and a map such that is a Cauchy sequence with modulus for any

2. a modulus of continuity , which satisfies

 |a−b|≤2−ωf(p)+1→|hf(a,n)−hf(b,n)|≤2−p

for any with monotone and .

As an example, consider the representation of the function on the interval in Minlog:
ContConstr 0 2 ([a,n]a*a) ([p] 0) ([p] p+3) where ContConstr is a constructor of a function of type cont which takes as the input the rational endpoints of the interval , the approximating map ([a,n]a*a) of type , the mapping of type as modulus of convergence and the mapping of type as modulus of continuity. Along with the constructor, a proof of Cont has to be provided for this representation. That is, according to the definition of Cont

• ([a,n]a*a) is Cauchy with the modulus because the values of the sequence are independent of .

• The modulus of continuity is correct and the corresponding goal

 ∀a∈[0,2]b∈[0,2]p∈P(|a−b|≤2−p−2→|a⋅a−b⋅b|≤2−p)

holds because .
Note that holds because .

• and are both monotone.

Note that the definition of the approximation map in Cont differs from the definition in Minlog in the sense that the approximation maps in Minlog are defined as programs of type , omitting the information of the interval . However this information is provided on the proof level by attaching the assumptions that the rational inputs have to belong to the specific interval .

###### Definition 3.5 (Application).

The application of a uniformly continuous function to a real in , is defined to be a Cauchy sequence with modulus

###### Theorem 3.6 (ApplicationCorr).

For a uniformly continuous function and a real in , the Cauchy sequence emanating from Application and accompanied with the modulus is a real number.

###### Remark 3.7.

The definition of the application incorporates the so-called extension theorem due to Bishop [10, p. 91]. In Minlog the application is realized by a program constant and allows one to address the lifting from rational mappings (see Definition 3.4) to the reals, which is at some point similar to the technique of the completion monad for the constructive reals in Coq [51]. Additionally, note that for the continuous function and real numbers under the usage of the program Application, one is able to prove:

 |x−y|≤2−ωf(p)→|f(x)−f(y)|≤2−p,

which coincides to another usual existing notions of continuity in effective analysis [10, 34, 36, 75].

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.

###### Definition 3.8 (Free algebra list).

A list xls containing the elements of type is either :

• (Nil alpha), also called empty,

• constructed xl::xls with head xl and tail xls.

###### Remark 3.9.

A list contaning one element xl is denoted by xl:.

###### Definition 3.10 (CauchyVector).

Consider the lists , . A pair is a CauchyVector if and either:

• and are empty,

• the pair (head , head ) is Cauchy and the pair (tail , tail ) is CauchyVector.

###### Definition 3.11 (MonVector).

A list is MonVector if either:

• are empty,

• the map head is Mon and the list tail is MonVector.

###### Definition 3.12 (RealVector).

For the lists , of length , a vector if the pair is CauchyVector and is MonVector.

Similarly to Definition 3.4, a function on multiple variables can be defined as follows:

###### Definition 3.13 (ContMV).

A uniformly continuous function of multiple variables on a compact ball with rational center and rational radius is given by

1. an approximation map and a map such that is a Cauchy sequence with modulus for any

2. a modulus of continuity, which satisfies

 ∥e1−e2∥≤2−ωg(p)+1→|hg(e1,n)−hg(e2,n)|≤2p

for any with monotone and .

###### Remark 3.14.

The choice of the norm influences the values of moduli and . Usually, one would prefer to use either -norm or -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.

In summary, in this section we presented the most important concepts of the Minlog system and presented the formalization of constructive analysis in 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.

## 4. Formal Proofs of Stability for Polynomial Systems

Once the SOS Lyapunov function has been obtained by Algorithm 1, 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 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 [59].

###### Theorem 4.1 (Richardson).

Let be a set of expressions representing real, single valued, partially defined functions of one real variable. Let be the set of functions represented by expressions in . If contains then the problem “given an expression in , decide if there is a real number such that is less than zero” is undecidable.

###### Remark 4.2.

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 [15] or if there exists an elementary undefined integral of an elementary function [60] are both undecidable.

###### Remark 4.3.

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 [55]. In this work, the polynomials are given directly as the continuous functions with polynomial mapping with rational coefficients and by the modulus of convergence . Thus such continuous functions are marked as polynomials and normalization-by-evaluation allows one often to work directly with the polynomial mapping . Due to Theorem 3.6 one is still able to address real number coefficients or inputs for the polynomials in this framework.

Now in order to incorporate the notions of positive definite and nonnegative functions for the verification of the output of Algorithm 1 (exactLyapunov), the appropriate definitions and the equivalent IDPCs are provided as follows:

###### Definition 4.4 (Nonnegative function (ContNonNeg)).

A uniformly continuous function given by is called nonnegative whenever

1. ,

2. ,

and the corresponding IDPC in Minlog: [commandchars=
{}] (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”))

###### Remark 4.5.

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:

 RealNNeg: ∀x Real x→ ∀p 0 ≤ x seq(x mod p)+2−p→ RealNNeg x

in Minlog, where seq and mod are the display tokens for the deconstructors of the type real.

###### Definition 4.6 (Positive Definite Function (PosDef)).

A uniformly continuous function given by is called positive definite whenever

1. ,

2. ,

3. ,

and the corresponding IDPC in Minlog:

[commandchars=
{}] (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”))

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 which actually derives the positive definiteness of function if all inputs of are considered to be rational vectors.

###### Definition 4.7 (PosDefRatWit).

A uniformly continuous function given by is called positive definite with rational witness whenever

1. ,

2. ,

3. ,

and the corresponding IDPC in Minlog: [commandchars=
{}] (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”))

.

###### Proof.

Without loss of generality consider 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 there exists such that . It remains to show that . We claim that . The goal is equivalent to . By applying Theorem 3.3 (RealPosChar) backwards with the parameter , the new goal is

 ∀n(αg(η(p)+1)≤n→12η(p)+1≤hg(e1(n),e2(n),n)).

Now, for any such that , the intermediate step is

 12η(p)−12η(p)+1≤ hg(e1(n),e2(n),n)+hg(e1(n),e2(n),αg(η(p)+4)) −hg(e1(n),e2(n),αg(η(p)+4)).

Observe that

 12η(p)≤hg(e1(n),e2(n),αg(η(p)+4))

holds by assumption after unfolding Definition 4.7  . In addition,

 −12η(p)+1≤hg(e1(n),e2(n),n)−hg(e1(n),e2(n),αg(η(p)+4))

holds since

 hg(e1(n),e2(n),αg(η(p)+4))−hg(e1(n),e2(n),n)≤ |hg(e1(n),e2(n),αg(η(p)+4))−hg(e1(n),e2(n),n)|≤12η(p)+1,

by applying Theorem 3.1 (CauchyElim) with the parameter and by the fact that whenever . ∎

The next section focuses on examples where we certify formally the stability of continuous/discrete-time dynamical systems.

## 5. Examples for verification of Lyapunov functions

Here we combine the previous formalization framework from Section 4 with the hybrid numeric/symbolic algorithm exactLyapunov, stated in Algorithm 1 from Section 2.3. Namely, for a given vector of polynomials, we first compute a weighted rational SOS function with exactLyapunov and then formally prove that is a positive definite form. If we succeed in proving that is a nonnegative (resp. positive definite) function then the system is formally proved to be stable (resp. asymptotically stable).
Example 1: We start with the following two-dimensional example of continuous-time polynomial system:

 ˙x1 =−x31+x2, (5.1) ˙x2 =−x1−x2.

A quadratic SOS Lyapunov function is obtained via exactLyapunov. Then, one can formally verify within Minlog via Theorem 4.8 that the following (in)-equalities hold:

so that the proof system can assert that and are PosDef. In order to do that it suffices to show that the function is PosDef. Without loss of generality assume that the function is defined on the unit ball together with the following functions:

• ,

• ,

• .

To obtain the map , required to prove PosDefRatWit, we proceed by means of program extraction. That is the first goal of the proof is:

 (5.2) ∀e(∃p0#pe→∃q0

Observe that the assumption is equivalent to the fact that which is equivalent to which is in turn equivalent to . The last two equivalences follow from the fact that . Additionally, the normalization-by-evaluation converts the goal to the goal , also due to the fact that and . We may now proceed by case distinction:
Case 1:
We have such that and set .
Case 2:
Similarly to Case 1: we have such that and set .
Case 3:
The goal follows from ex-falso-quodlibet since .
Now let be a program which generates for any such that . In Minlog the extracted is represented by: [commandchars=
{}] [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)])])])] which takes as input a rational vector and a positive number . If , it generates a positive number such that . Trivially, PosDefRatWit follows from the derivation above. By Theorem 4.8, the function is also PosDef. Similarly, ignoring the factor , for , we obtain the following extracted : [commandchars=
{}] [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)])])])]

Overall, we have formally proved that the system is asymptotically stable.
Example 2: Let us consider another continuous-time polynomial system from [54, § 4.2]:

 ˙x1 =−x31−x1x23, (5.3) ˙x2 =−x2−x21x2, ˙x3 =−x3−3x3x23+1x1+3x21x3.

A homogeneous quadratic SOS Lyapunov is computed with exactLyapunov. Then, one formally verifies within Minlog that the following (in)-equalities hold:

 V3(x1,x2,x3)−(x21+x22+x23) =4.5489x21+3.1068x22+0.7945x23≥0, −gradV3⋅(x23+1)f =55489/5000(569705563(2x2/55489)2 +995750105(2/55489x3−495266/24893752625x21x3)2 +(x21+119248x231387225)2+569705563(2x1x2/55489)2 +98987825902(x1x3/277445)2+569705563(2x2x3/55489)2 +38212870902320486809(x23/10996532575)2 +375452958989124261245(x21x3/24893752625)2 +569705563(2x1x2x3/55489)2+91778806(x1x23/55489)2)≥0,

so that the proof system can assert that is PosDef and that is NonNeg. For this later goal, Theorem 3.2 (RealNNegChar) is applied backwards and the resulting goal is read as:

 ∀p∃n0∀n≥n0−2p≤ 55489/5000(569705563(2e2(n)/55489)2 +995750105(2/55489e3(n)−495266/24893752625e1(n)2e3(n))2 +(e1(n)2+119248e3(n)21387225)2+569705563(2e1(n)e2(n)/55489)2 +98987825902(e1(n)e3(n)/277445)2+569705563(2e2(n)e3(n)/55489)2 +38212870902320486809(e3(n)2/10996532575)2 +375452958989124261245(e1(n)2x3/24893752625)2 +569705563(2e1(n)e2(n)e3(n)/55489)2+91778806(e1(n)e3(n)2/55489)2).

For any , after introducing , the goal is trivial by proving the nonnegativity of each summand separately. Overall, we have formally proved that the system is stable.
Example 3: Eventually, we take the following two-dimensional discrete-time rational system:

 xt+1 =yt1+x2t, (5.4) yt+1 =xt1+y2t.

In a similar way as for continuous-time polynomial systems, one can compute a homogeneous quadratic SOS Lyapunov with exactLyapunov, then formally verify within Minlog that the following (in)-equalities hold:

 V2(x,y) =x2+y2≻0, (1+x2)2(1+y2)2(V2(x,y)−V2(f(x,y))) =