Computation of Bounds for Polynomial Dynamic Systems

1. Introduction

Lyapunov’s direct method enables stability analysis without knowing the explicit solution of a (usually nonlinear) differential equation [1,2]. This method was originally introduced for stability analysis of equilibrium points, but has been extended to numerous related problems in recent decades. In addition to control theory applications, e.g. with control Lyapunov functions, Lyapunov’s method can also be used to estimate basins of attraction, trapping regions and attractors [3,4,5].

If a nonlinear dynamic system has a chaotic attractor, then the geometry of this attractor usually cannot be described exactly in a straightforward manner. For many questions, such as estimating the attractor’s dimension [6,7], it is necessary or helpful to narrow down this attractor in the state space. Such an containment can be realized via positive invariant sets, which are described by sublevel sets of a Lyapunov-like function [3,4,5]. In this case, one starts from a positive definite Lypunov-like candidate function and must show that its time derivative along the dynamic system in question is negative outside the sublevel set. There are two difficulties here: firstly, choosing the correct candidate function, and secondly, showing that its total time derivative is negative for large values. In particular, the second step usually requires extensive experience and intuition.

The positive definiteness of the candidate function can usually be ensured or easily verified. The property that the time derivative of the candidate function is negative in certain regions can be described by combining inequalities with quantifiers. For polynomial systems with a polynomial Lyapunov candidate function, these inequalities are polynomial as well and the resulting formulas with quantifiers can be solved with a technique known as quantifier elimination (QE) [8,9,10,11].

The calculation of bounds for positively invariant sets using QE has already been successfully applied to the Lorenz system, the Lorenz-Haken system, and a system going back to Rössler [12,13]. This approach was successful because the Lyapunov candidate functions used were fairly simple. Minor generalizations of the Lyapunov candidate functions would quickly result in a situation where the required calculations could no longer be performed within a reasonable amount of time. For further applications of the proposed approach, it is therefore important to reduce the computational effort. In [14], a method was developed to simplify the underlying quantifier elimination problem when searching for suitable Lyapunov functions. This method was reduced to the two-dimensional case and focused on stability analysis of equilibrium points. This paper extends the approach to systems of arbitrary finite dimension. In addition to analyzing the stability of equilibrium points, it is also possible to determine bounds for positively invariant sets.

The paper is structured as follows. In Section 2 we provide mathematical preliminaries. The Lorenz system is discussed in Section 3 as a motivational example, where the choice of the Lyapunov-like candidate function has an impact on the computability. The main result is presented in Section 4. In Section 5 we illustrate the approach on Chua’s system with cubic nonlinearity. The results are discussed in Section 6.

2. Mathematical Preliminaries

2.1. Nonlinear Systems and Positive Invariant Sets

Consider a nonlinear autonomous system

$$\dot{x}=f\left(x\right)$$

(1)

of ordinary differential equations with the vector field $f:{\mathbb{R}}^n\to {\mathbb{R}}^n$ with the state vector x. For the computational methods derived in the paper we assume that the vector field f is polynomial, i.e., its components are multivariate polynomials in the variables $x_1,\dots ,x_n$ with rational coefficients. From this follows that f is smooth and locally Lipschitz. Consequently, for every initial value $x\left(0\right)\in {\mathbb{R}}^n$ there exists an at least locally defined unique solution $x\left(t\right)$. The existence of a global solution will become apparent later, when bounds of the solution have been shown.

A subset $\mathrm{\Omega }\subseteq {\mathbb{R}}^n$ is called positive invariant with respect to (1) if every solution of (1) starting in $\mathrm{\Omega }$ stays in $\mathrm{\Omega }$ for all future times, i.e.,

$$\forall t\ge 0:x\left(0\right)\in \mathrm{\Omega }\Rightarrow x\left(t\right)\in \mathrm{\Omega }.$$

(2)

Such sets are also called trapping regions [15]. Positive invariant sets can be approximated using Lyapunov-like candidate functions $V:{\mathbb{R}}^n\to \mathbb{R}$ and a constant $\gamma >0$ by

$$\mathrm{\Omega }=\{x\in {\mathbb{R}}^n:V\left(x\right)\le \gamma \},$$

(3)

where the subset $\mathrm{\Omega }\subset \mathbb{R}$ is bounded by the level set $V\left(x\right)=\gamma$. If all solutions of (1) starting with initial values $x\left(0\right)\in {\mathbb{R}}^n\setminus \mathrm{\Omega }$ with $V\left(x\right)>\gamma$ tend to $V\left(x\right)\le \gamma$, the set $\mathrm{\Omega }$ contains an attractor. In this case, the solution of (1) is bounded by

$$\underset{t\to \infty }{lim\; sup}V\left(x\left(t\right)\right)\le \gamma .$$

(4)

Without explicit knowledge of the solution $x(·)$ of (1), the bound $\gamma$ describing the subset (3) can be computed by

$$\forall x\in {\mathbb{R}}^n:V\left(x\right)>\gamma \Rightarrow \dot{V}\left(x\right)<0$$

(5)

using Lyapunov techniques, where the derivative $\dot{V}$ is the directional derivative of V along the dynamics of system (1). This derivative can be written as a Lie derivative of the scalar field V along the vector field f, i.e.,

$$\dot{V}\left(x\right)=L_{f}V\left(x\right)={\displaystyle \frac{\partial V\left(x\right)}{\partial x}}f\left(x\right),$$

(6)

see [2]. Alternatively, the bound $\gamma$ can be calculated using the differential inequality

$$\exists q>0\forall x\in {\mathbb{R}}^n:\dot{V}\left(x\right)\le -q·(V\left(x\right)-\gamma )$$

(7)

as suggested in [16]. The result is an interval for $\gamma$, from which the minimum or infimum is selected as the best bound.

Example 1. 

Quadratic forms

$$V\left(x\right)=x^{T}Px$$

(8)

with a positive definite symmetric matrix $P=P^T\in {\mathbb{R}}^{n\times n}$ are often used as a candidate function for V. When the matrix P contains free parameters, the definiteness can be ensured. For example, by Sylvester’s criterion or a Cholesky decomposition [17,18,19]. A construction method for more general candidate functions is based on sum of squares (SOS) decomposition [20], but leads to a purely numerical procedure. However, it is much more difficult to solve (5) or (7) for the bound γ. In particular, formulas (5) or (7) are only solvable if $\dot{V}\left(x\right)<0$ for large $\parallel x\parallel >0$.

2.2. Quantifier Elimination

The formulas (5) and (7) specify the bound $\gamma$ for the positive invariant set (3). These formulas contain quantifiers, namely the universal quantifier ∀ and the existential quantifier ∃. These formulas can be solved by quantifier elimination.

Let us briefly introduce some notations and concepts [9]. The starting point for our considerations are polynomial equations or inequalities over the reals of the form

$$\phi (X_1,\dots ,X_k)⋉0\mathrm{with}⋉\in \{=,\ne ,<,>,\le ,\ge \}$$

(9)

called atomic formulas in the variables $X=(X_1,\dots ,X_k)$. A Boolean combination of finitely many atomic formulas using logical operators such as $\neg ,\wedge ,\vee ,\Rightarrow ,\dots$ is called a quantifier-free formula. A quantifier-free formula describes a semi-algebraic set. Let $F(X,Y)$ be a quantifier-free formula in the variables X and the additional variables $Y=(Y_1,\dots ,Y_l)$. A prenex formula has the form

$$G(X,Y)=\left(Q_1{Y}_1\right)\dots \left(Q_l{Y}_l\right)F(X,Y)$$

(10)

with quantifiers $Q_1,\dots ,Q_l\in \{\forall ,\exists \}$. The variables Y are called quantified variables. The remaining variables not bound by a quantifier are called free variables. The Tarski-Seidenberg theorem [8,9,21] states, that for every prenex formula $G(X,Y)$ there exists an equivalent quantifier-free formula $H\left(X\right)$ solely depending on the free variables X. This transformation from a prenex formula to an equivalent quantifier-free formula is called quantifier elimination (QE).

Over the past half-century, a number of different methods for performing quantifier elimination have been developed. The best-known methods are the cylindrical algebraic decomposition (CAD) [22,23], the virtual substitution (VS) [24,25] and the real root classification [26]. However, all these methods are very computationally intensive. In order to apply QE to practical problems, it is very important to find an appropriate mathematical formulation.

Several software packages are available for solving non-trivial QE problems. The open source package QEPCAD B [27] uses CAD and is available in the repositories of common Linux distributions. The package REDLOG [28] is based on the computer algebra system REDUCE and contains implementations of CAD as well as VS. The computed quantifier-free expressions can be simplified with SLFQ [29]. Commercial computer algebra systems such as Maple or Mathematica also have procedures and libraries for carrying out QE, e.g. [30,31].

Example 2. 

Consider the real biquadratic polynomial

$$g\left(x\right)=x^4+p{x}^2+q$$

(11)

with the parameters p and q. The question of the existence of at least one real root can be described by the following prenex formula

$$\exists x:g\left(x\right)=0,$$

(12)

where x is a quantified variable and $p,q$ are free variables. Applying QE to (12) leads to the equivalent quantifier-free formula

$$4q-p^2\le 0\wedge (p<0\vee q\le 0).$$

(13)

Similarly, one can ask under what conditions on the parameters the polynomial only takes positive function values:

$$\forall x:g\left(x\right)>0.$$

(14)

Carrying out QE to (14) leads to the equivalent quantifier-free formula

$$q>0\wedge (p>0\vee 4q-p^2>0).$$

(15)

3. Motivational Example

The Lorenz system is one of the most popular nonlinear systems exhibiting chaotic dynamics. The system is governed by the differential equations

$$\begin{array}{ccc}{\dot{x}}_1\hfill & =& s(x_2-x_1)\hfill \\ {\dot{x}}_2\hfill & =& r{x}_1-x_2-x_1{x}_3\hfill \\ {\dot{x}}_3\hfill & =& x_1{x}_2-b{x}_3\hfill \end{array}$$

(16)

with positive parameters $s,r,b>0$. To calculate bounds for the Lorenz attractor, numerous Lyapunov candidate functions were investigated [16,32,33,34]. Two well-known candidate functions were given in [5] and [4, Appendix C]:

$$\begin{array}{ccc}\hfill V_a\left(x\right)& =& x_1^2+x_2^2+{(x_3-r-s)}^2,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill V_b\left(x\right)& =& V_3\left(x\right)=r{x}_1^2+s{x}_2^2+s{(x_3-2r)}^2.\hfill \end{array}$$

(18)

The level sets of these functions describe a sphere and an ellipsoid, respectively, both shifted with respect to the $x_3$-axis. Another common feature is that the coefficients of the monomials $x_2^2$ and $x_3^2$ are the same within each equation. If we ignore this commonality for the time being, we can consider the functions (17) and (18) as special cases of the Lyapunov candidate function

$$V_c\left(x\right)=p_1{x}_1^2+p_2{x}_2^2+p_3{(x_3-x_{30})}^2,$$

(19)

with the coefficients $p_1,p_2,p_3>0$ and $x_{30}\in \mathbb{R}$, cf. [12]. The time derivative of $V_c$ along (16) is given by

$$\begin{array}{ccc}{\dot{V}}_c\left(x\right)\hfill & =& -2\left(p_2-p_3\right)x_1{x}_2{x}_3-2p_{1}s{x}_1^2-2p_2{x}_2^2-2b{p}_3{x}_3^2\hfill \\ \hfill & & +2\left(p_{2}r+p_{1}s\right)x_1{x}_2-2\left(p_3{x}_1{x}_2-b{p}_3{x}_3\right)x_{30}.\hfill \end{array}$$

(20)

The free parameters should be chosen so that ${\dot{V}}_c\left(x\right)<0$ is possible for large $\parallel x\parallel \gg 0$. The cubic term $-2\left(p_2-p_3\right)x_1{x}_2{x}_3$ is the one with the highest degree. However, the monomial $x_1{x}_2{x}_3$ has an odd degree. This monomial can be both positive and negative for large $\parallel x\parallel$. However, if we choose $p_2=p_3$, the derivative (20) can become negative for large $\parallel x\parallel$ due to the quadratic terms. This explains the above-mentioned commonality in the candidate functions (17) and (18).

If we now ignore the displacement along the $x_3$-axis, then (19) represents a special case of a quadratic form (8), where the matrix P is a diagonal matrix. The question arises whether a better bound on the invariant set is possible with a dense matrix

$$P=\left(\begin{array}{ccc}{p}_{11}& p_{12}& p_{13}\\ p_{12}& p_{22}& p_{23}\\ p_{13}& p_{23}& p_{33}\end{array}\right).$$

(21)

The matrix P is positive definite if and only if

$$a_{11}>0\wedge a_{11}{a}_{22}-a_{12}^2>0\wedge det\left(P\right)>0$$

(22)

according to Sylvester’s criterion [17,18].

The Lie derivative of (8) with the generic matrix (21) is a polynomial in the variables $x_1,x_2,x_3$ with total degree 3. The terms with the highest degree are

$$2p_{13}{x}_1^2{x}_2+2p_{23}{x}_1{x}_2^2-2p_{12}{x}_1^2{x}_3-2p_{23}{x}_1{x}_3^2+2(p_{33}-p_{22})x_1{x}_2{x}_3.$$

(23)

Since the highest degree is an odd number, the associated monomials

$$x_1^2{x}_2,x_1{x}_2^2,x_1^2{x}_3,x_1{x}_3^2,x_1{x}_2{x}_3$$

(24)

are indefinite w.r.t. their sign for large $\parallel x\parallel$. If these monomials occur in $\dot{V}$, the condition $\dot{V}\left(x\right)<0$ for $\parallel x\parallel \gg 0$ will not be fulfilled. Therefore, the entries of the matrix (21) have to be chosen such that the terms (23) vanish. This is achieved by choosing the coefficients as follows

$$p_{12}=0\wedge p_{13}=0\wedge p_{23}=0\wedge p_{22}=p_{33}$$

(25)

resulting in the diagonal matrix

$$P=\left(\begin{array}{ccc}{p}_{11}& 0& 0\\ 0& p_{22}& 0\\ 0& 0& p_{22}\end{array}\right)$$

with $p_{11},p_{22}>0$.

Remark 1. 

The choice (25) corresponds exactly to the coefficients in (19). This means that the dense matrix (21) does not provide any advantage for evaluating the global attraction of (3), i.e., for the behavior of the system for large $\parallel x\parallel$. This does not rule out the possibility that a dense matrix P could also have advantages for the more precise determination of bounds for invariant sets for smaller $\parallel x\parallel$.

Remark 2. 

For the candidate Lyapunov functions (17), (18), (19) as well as (8) with (21), bounds γ can theoretically be determined exactly (with reference to the respective candidate function) in finite time by applying QE to the formulas (27) and (31). For the candidate function (17), (18), symbolic expressions for the bound γ as a function of the symbolic parameters $s,r,b$ can be obtained on a standard PC in a few seconds, whereby case distinctions occur for different parameter ranges [12]. When using the candidate function (19), the situation is different. Neither with symbolic nor with numerically specified system parameters $s,r,b$ did QE deliver a solution, even after several days of computing time. Only when numerical values were specified for the parameters $p_1,p_2,p_3,p_{30},\gamma$, thereby creating a decision problem, did QE deliver a result. This also meant that it was not possible to directly verify whether the more general approach (8) with (21) offers any advantages.

Problem 1. 

Remark 2 implies that the candidate function must be selected in such a way that, on the one hand, there are enough degrees of freedom to meet the definiteness requirements, but on the other hand, there are not too many free parameters to be able to check these definiteness conditions with QE in reasonable time.

4. Results

Section 4.1 deals with the use of QE to validate the local definiteness of functions, as used in the stability analysis of steady states. For this, Section 4.2 describes a necessary condition that requires significantly less computing power. In this process, the terms of lowest degree of a polynomial are evaluated consecutively. If one proceeds in a similar manner in reverse, starting with terms of the highest degree and working downwards, one obtains conditions for Lyapunov candidate functions that are suitable for estimating positive invariant sets, see Section 4.3.

4.1. Description of Local Positive Definiteness by Prenex Formulas

To formulate the stability conditions, we need the concept of positive definiteness of a Lyapunov candidate function V and negative definiteness or semi-definiteness of its time derivative $\dot{V}$. The algorithms refer to both V and $\dot{V}$, so we will use the designation $W\in \{V,\dot{V}\}$ in the following. Since we are dealing with polynomials in this paper, all functions are globally defined and continuously differentiable.

Definition 1. 

A function $W:{\mathbb{R}}^n\to \mathbb{R}$ is called locally positive definite or positive definite on a neighborhood $\mathcal{U}$ around the origin if

$$W\left(0\right)=0\wedge \forall x\in \mathcal{U}\setminus \left\{0\right\}:W\left(x\right)>0.$$

(26)

If (26) holds on $\mathcal{U}={\mathbb{R}}^n$, the function is called globally positive definite.

In the same manner we define positive semi-definite with $W\left(x\right)\ge 0$, negative definite with $W\left(x\right)<0$ and negative semi-definite with $W\left(x\right)\le 0$ for all $x\in \mathcal{U}$. Clearly, a function W is negative definite (semi-definite) if $-W$ is positive definite (semi-definite). A function W is called indefinite if the function is neither positive semi-definite nor negative semi-definite, i.e., for every neighborhood $\mathcal{U}$ there exists two points $x,y\in \mathcal{U}\setminus \left\{0\right\}$ with $W\left(x\right)>0$ and $W\left(y\right)<0$.

We assume that the function W is a multivariate polynomial without a constant term, which ensures $W\left(0\right)=0$. If we characterize the neighborhood $\mathcal{U}=\{x\in {\mathbb{R}}^n:\parallel x{\parallel }_{\infty }<\epsilon \}$ with $\epsilon >0$ using the maximum norm, the test for a given candidate function W regarding local positive definiteness on $\mathcal{U}$ can be formulated with the following prenex formula:

$$\exists \epsilon >0\forall x_1\dots \forall x_n:|x_1|\ge \epsilon \vee \dots \vee |x_n|\ge \epsilon \vee (x_1=0\wedge \dots \wedge x_n=0)\vee W\left(x\right)>0.$$

(27)

If W is fully specified without free parameters, the prenex formula (27) constitutes a decision problem, where QE leads either to true or false. If W contains free parameters (e.g. (8) with a symbolic matrix (21)), QE would result in a quantifier-free formula with conditions on the free parameters. In principle, after applying QE, the prenex formula (27) would provide a conclusive test for positive definiteness or a characterization of the free parameter values required for it. However, it is not clear whether, for a given function, QE can be performed within a reasonable time due to the computational effort involved.

Consider a region ${\mathcal{D}}_i\subset {\mathbb{R}}^n\setminus \left\{0\right\}$ of the state space defined by

$$|x_i|\ge |x_j|\mathrm{for}j\ne i.$$

(28)

In this region, the substitution

$$(x_1,\dots ,x_n)\leftarrow ({\lambda }_1{x}_i,\dots ,{\lambda }_{i-1}{x}_i,x_i,{\lambda }_i{x}_i,\dots ,{\lambda }_{n-1}{x}_i)$$

(29)

with ${\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]$ defines the mapping

$$(x_1,\dots ,x_n)\mapsto \left(x_i,{\lambda }_1={\displaystyle \frac{x_1}{x_i}},\dots ,{\lambda }_{i-1}={\displaystyle \frac{x_{i-1}}{x_i}},{\lambda }_i={\displaystyle \frac{x_{i+1}}{x_i}},\dots ,{\lambda }_{n-1}={\displaystyle \frac{x_n}{x_i}}\right)$$

and describes a change of coordinates between $(x_1,\dots ,x_n)\in {\mathcal{D}}_i$ and $(x_i,{\lambda }_1,\dots ,{\lambda }_{n-1})\in \mathbb{R}\setminus \left\{0\right\}\times {[-1,1]}^{n-1}$. The Lyapunov candidate function W can be written on ${\mathcal{D}}_i$ in the new coordinates as follows:

$$W_i(x_i,{\lambda }_1,\dots ,{\lambda }_{n-1})=W({\lambda }_1{x}_i,\dots ,{\lambda }_{i-1}{x}_i,x_i,{\lambda }_i{x}_i,\dots ,{\lambda }_{n-1}{x}_i).$$

(30)

The regions ${\mathcal{D}}_1,\dots ,{\mathcal{D}}_n$ describe a decomposition of the state space ${\mathcal{D}}_1\cup \dots \cup {\mathcal{D}}_n={\mathbb{R}}^n\setminus \left\{0\right\}$ as depicted in Figure 1 for $n=3$. Condition (27) can be tested for each region separately using the new coordinates by

$$\exists \epsilon >0\forall z\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]:\left|z\right|\ge \epsilon \vee z=0\vee W_i(z,{\lambda }_1,\dots ,{\lambda }_{n-1})>0$$

(31)

for $i=1,\dots ,n$. Note that the prenex formula (31) still has the same number of quantified variables as (27). The restriction to one sector could in some cases lead to a simplification of QE, but on the other hand, the degree of some terms will increase due to the multiplication of $\lambda$ and x in (29), which can lead to a higher computational effort.

Remark 3. 

If the auxiliary variables ${\lambda }_1,\dots ,{\lambda }_{n-1}$ are fixed (or assumed to have already been eliminated), then the function $W_i$ is a univariate polynomial in the variable z. The definiteness of a single term $a_k{z}^k$ of a univariate polynomial can be easily determined. For an even number k, the term is positive or negative definite if and only if $a_k>0$ or $a_k<0$, respectively. For an odd number k, the term with $a_k\ne 0$ is indefinite. Transforming a multivariate polynomial into a univariate polynomial therefore simplifies the definiteness test.

4.2. Necessary Condition

If the order of consecutive quantified variables of the same quantifier is swapped, the logical statement in question does not change. Changing the order of variables with different quantifiers, on the other hand, changes the statement. Swapping the existential quantifier with the universal quantifier results in the following implication:

$$\exists \epsilon \forall \lambda :F\Rightarrow \forall \lambda \exists \epsilon :F$$

(32)

Now we apply this implication to (31) and obtain the following necessary condition

$$\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]\underset{{\displaystyle =:H}}{\underbrace{\exists \epsilon >0\forall z:\left|z\right|\ge \epsilon \vee z=0\vee W_i(z,{\lambda }_1,\dots ,{\lambda }_{n-1})>0}}$$

(33)

with the prenex formula H. The function $W_i$ can be written as a polynomial

$$W_i(z,{\lambda }_1,\dots ,{\lambda }_{n-1})=a_0\left(\lambda \right)+a_1\left(\lambda \right)z+a_2\left(\lambda \right)z^2+\cdots +a_M\left(\lambda \right)z^M$$

(34)

in the variable z with coefficients $a_0,\dots ,a_M\in \mathbb{R}[{\lambda }_1,\dots ,{\lambda }_{n-1}]$ from the ring of real polynomials in $\lambda =({\lambda }_1,\dots ,{\lambda }_{n-1})$. Without loss of generality we assume that the constant term vanishes, i.e., $a_0\left(\lambda \right)=0$. For fixed values $\lambda \in {\mathbb{R}}^{n-1}$, the prenex formula H simply states the local positive definiteness of the function $z\mapsto W_i(z,\lambda )$, i.e., $W_i(·,\lambda )$ has a local minimum at $z=0$. For univariate polynomials, this means that the term of lowest degree is even with a positive coefficient. From (33), the formula

$$\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]:\left[(a_1=0\wedge a_2>0)\vee (a_1=0\wedge \dots \wedge a_3=0\wedge a_4>0)\vee \cdots \right]$$

(35)

is obtained. Formula (35) is therefore a necessary condition for $W_i$ to be locally positive definite. Necessary condition for $W_i$ to be locally negative definite as well as locally positive or negative semi-definite can be formulated similarly.

Remark 4. 

For practical implementation, it would be useful to check the various conditions for the coefficients step by step:

$$\begin{array}{c}{a}_1=0\wedge a_2>0\hfill \\ a_1=0\wedge \dots \wedge a_3=0\wedge a_4>0\hfill \\ a_1=0\wedge \dots \wedge a_5=0\wedge a_6>0\hfill \\ ⋮\hfill \end{array}$$

If the first condition is fulfilled, then the function is locally positive definite and the algorithm is terminated. Otherwise, the subsequent conditions must be checked.

Remark 5. 

Although relation (35) is only necessary for (31), it can be used to restrict the degrees of freedom in the choice of parameters for the Lyapunov candidate function and thus possibly eliminate superfluous terms. From a computational point of view, the formulation (35) has several advantageous compared to (26). On the one hand, the number of quantified variables has been reduced by two. In addition, there is no change between different quantifier types, which is particularly advantageous for performing QE with VS [25]. On the other hand, formulation (35) is based only on polynomials in λ, so that omitting z results in terms of lower degree and thus allows for faster computation.

4.3. Description of Definiteness at Large Values

To check the stability of an equilibrium (typically at the origin), local definiteness conditions must be checked in the neighborhood of this equilibrium. These conditions apply to small values of $\parallel x\parallel$. For determining bounds for positively invariant sets, however, the behavior for large $\parallel x\parallel$ is of interest. The following form of definiteness is essentially complementary to local definiteness [35].

Definition 2. 

A function $W:{\mathbb{R}}^n\to \mathbb{R}$ is called  positive definite at infinity  or  positive definite at large values  if there exists a number $r>0$ with

$$\forall x:\parallel x\parallel >r\Rightarrow W\left(x\right)>0.$$

(36)

A function $W:{\mathbb{R}}^n\to \mathbb{R}$ is called (positive) coercive [36] if $W\left(x\right)\to \infty$ for $\parallel x\parallel \to \infty$.

In a similar way, the terms positive semi-definite, negative definite and negative semi-definite at infinity can be introduced. Clearly, a polynomial coercive function is positive definite at infinity. The use of a coercive function is recommended in itself, since the sets (3) defined by the preimages of such a function are compact.

Example 3. 

The function $W:{\mathbb{R}}^2\to \mathbb{R}$ with $W\left(x\right)=1+{(x_1-x_2)}^2$ is positive definite at infinity, but not coercive.

A minor modification of (27) leads to a test of a given candidate function W regarding positive definiteness at infinity with the following prenex formula:

$$\exists \epsilon >0\forall x_1\cdots \forall x_n:\left(|x_1|\le \epsilon \wedge \cdots \wedge |x_n|\le \epsilon \right)\vee W\left(x\right)>0.$$

(37)

As in Section 4.1, the Lyapunov-like candidate function W can be rewritten on a region ${\mathcal{D}}_i$ according to (30). Then, condition (37) leads to

$$\exists \epsilon >0\forall z\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]:\left|z\right|\le \epsilon \vee W_i(z,{\lambda }_1,\dots ,{\lambda }_{n-1})>0$$

(38)

for $i=1,\dots ,n$. Applying the implication (32) yields the necessary condition

$$\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]\underset{{\displaystyle =:J}}{\underbrace{\exists \epsilon >0\forall z:\left|z\right|\le \epsilon \vee W_i(z,{\lambda }_1,\dots ,{\lambda }_{n-1})>0}}$$

(39)

with the prenex formula J. The prenex formula J simply states the positive definiteness of the function $z\mapsto W_i(z,\lambda )$ at infinity. As before, we consider $W_i$ as an univariate polynomial (34) of degree M in the variable z. For an even degree M, the prenex formula

$$\begin{array}{c}\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]:\hfill \\ \left[a_M>0\vee (a_M=0\wedge a_{M-1}=0\wedge a_{M-2}>0)\vee \cdots \right]\hfill \end{array}$$

(40)

is a necessary condition for $W_i$ to be positive definite at infinity. For an odd degree M, the prenex formula

$$\begin{array}{c}\forall {\lambda }_1,\dots ,{\lambda }_{n-1}\in [-1,1]:\hfill \\ \left[(a_M=0\wedge a_{M-1}>0)\vee (a_M=0\wedge \dots \wedge a_{M-2}=0\wedge a_{M-3}>0)\vee \dots \right]\hfill \end{array}$$

(41)

is a necessary condition for $W_i$ to be positive definite at infinity.

Remark 6. 

As stated in Remark 4, it may be useful for practical implementation to check the conditions on the coefficients step by step. If M is even, the conditions are obtained from  (40) as follows:

$$\begin{array}{c}{a}_M>0\hfill \\ a_M=0\wedge a_{M-1}=0\wedge a_{M-2}>0\hfill \\ a_M=0\wedge \dots \wedge a_{M-3}=0\wedge a_{M-4}>0\hfill \\ ⋮\hfill \end{array}$$

If the first condition is met, the algorithm can be terminated. Otherwise, the subsequent conditions would need to be examined.

5. Chua’s Circuit with Cubic Nonlinearity

In Section 5.1, we derive a model for Chua’s circuit in which the piecewise-linear nonlinearity is replaced by a cubic polynomial. The conditions for the definiteness of functions at infinity derived in Section 4.3 are applied in Section 5.2 to derive conditions on the structure of a Lyapunov-type candidate function for the Chua system. Admissible ranges for the parameter values of the Lyapunov-like candidate function are determined in Section 5.3. The computation of bounds for the containment of positively invariant sets is performed in Section 5.4.

5.1. System Model

Chua’s circuit is an electronic circuit that can generate a chaotic attractor known as double scroll when appropriate parameter values are used [37,38]. The circuit diagram is shown in Figure 2. Kirchhoff’s circuit laws yields the following differential equations

$$\begin{array}{ccc}{\dot{x}}_1\hfill & =& {\displaystyle \frac{G}{C_1}}(x_2-x_1)-{\displaystyle \frac{1}{C_1}}\varphi \left(x_1\right)\hfill \\ {\dot{x}}_2\hfill & =& {\displaystyle \frac{G}{C_2}}(x_1-x_2)+{\displaystyle \frac{1}{C_2}}{x}_3\hfill \\ {\dot{x}}_3\hfill & =& -{\displaystyle \frac{1}{L}}{x}_2\hfill \end{array}$$

(42)

with $x={(x_1,x_2,x_3)}^T={(v_1,v_2,i)}^T$ and the conductance $G=1/R$ and $C_1,C_2,L>0$. In [37,39], the authors used the normalized parameter values $C_1=1/9$, $C_2=1$, $L=1/7$ and $G=0.7$.

The nonlinearity $\varphi$ of the circuit is a voltage controlled current source known as Chua’s diode and labeled D in the circuit diagram. Typically, the nonlinearity is a 3-segment piecewise-linear function

$$\varphi \left(x_1\right)=m_0{x}_1+{\displaystyle \frac{1}{2}}(m_1-m_0)\left(\right|x_1+1|-|x_1-1\left|\right).$$

(43)

In [39], the parameter values $m=-0.5$ and $m_1=-0.8$ were used. The replacement of the non-smooth function (43) by cubic polynomial

$$\varphi \left(x_1\right)=a{x}_1+b{x}_1^3\mathrm{with}a<0\wedge b>0$$

(44)

was suggested in [40]. In [41, p. 143], the parameters a and b were chosen such that the functions (43) and (44) assume the smallest distance in the sense of quadratically integrable functions on the interval $[-2,2]$:

$$a=-{\displaystyle \frac{4183}{5120}}\approx -0.81699\mathrm{and}b={\displaystyle \frac{189}{4096}}\approx 0.046143.$$

(45)

Figure 3 shows both characteristics, which are in very good agreement over the interval $[-2,2]$, but diverge outside this range.

To illustrate the system’s behavior we carried out a numerical simulation with SageMath [42] using the parameters given above. Using the initial conditions $x_1\left(0\right)=x_2\left(0\right)=x_3\left(0\right)=0.1$, the transients seem to converge to the chaotic attractor as shown in Figure 4 (left). With the initial conditions $x_1\left(0\right)=x_2\left(0\right)=x_3\left(0\right)=1$, the solution, however, converges towards a limit cycle. This transient response is shown in Figure 4 (right). The axes scalings in Figure 4 show that the limit cycle lies further out in the state space than the chaotic attractor.

5.2. Structure of the Quadratic Lyapunov-like Candidate Function

Consider the Chua system (42) with the cubic nonlinearity (44) and the additional constraints

$$G>0\wedge C_1>0\wedge C_2>0\wedge L>0\wedge a<0\wedge b>0$$

(46)

on the parameters. To derive bounds on a positive invariant set we use a Lyapunov-like candidate function V as a quadratic form (8) with the matrix (21). The function V is positive definite if and only if (22) holds. The positive definiteness of V can thus be easily ensured. The question of whether or how the entries $p_{ij}$ of the matrix P can or should be chosen so that the time derivative $\dot{V}$ becomes negative definite or at least semi-definite for large values is much more difficult.

The time derivative of the Lyapunov-like candidate function V has the form

$$\dot{V}\left(x\right)={\displaystyle \frac{1}{C_1{C}_{2}L}}W\left(x\right),$$

(47)

where the numerator is a multivariate polynomial W in the variables $x_1,x_2,x_3$ and the parameters. The denominator is a positive constant. The definiteness of $\dot{V}$ can therefore be determined by the definiteness of W. The expression for W is very extensive and is therefore not given here.

Now, we will decompose the state space into three regions ${\mathcal{D}}_1,{\mathcal{D}}_2,{\mathcal{D}}_3$ as already depicted in Figure 1. Since the vector field f is a polynomial of degree 3 and V is quadratic in the variables $x_1,x_2,x_3$, the Lie derivative (6) is a polynomial of degree 4. Therefore, the functions $W_1,W_2,W_3$ defined by (30) on the three regions are polynomials of the degree 4 in the variable z having the form

$$W_i(z,{\lambda }_1,{\lambda }_2)=a_0\left(\lambda \right)+a_1\left(\lambda \right)z+a_2\left(\lambda \right)z^2+a_3\left(\lambda \right)z^3+a_4\left(\lambda \right)z^4$$

(48)

with $\lambda =({\lambda }_1,{\lambda }_2)$ and $a_0,\dots ,a_4\in \mathbb{R}\left[\lambda \right]$.

For large $\parallel x\parallel \gg 0$, only the terms of degree 4 are relevant. Therefore, the terms of degree 4 should not be positive for all ${\lambda }_1,{\lambda }_2\in [-1,1]$. This corresponds to the property of positive semi-definiteness of $W_i$ at infinity. According to Remark 6 one would therefore check the leading coefficient regarding to its sign first. Together with (46), we therefore obtain the prenex formulas

$$\forall {\lambda }_1,{\lambda }_2:\mathrm{Cond}.\left(46\right)\wedge \left[(-1\le {\lambda }_1\le 1\wedge -1\le {\lambda }_2\le 1)\Rightarrow a_4({\lambda }_1,{\lambda }_2)\le 0\right]$$

(49)

for $i=1,2,3$ similar to (40). The QE was carried out with REDLOG. The resulting quantifier-free expressions were combined using conjunction (AND) and subsequently simplified with SLFQ leading to

$$p_{11}\ge 0\wedge p_{12}=0\wedge p_{13}=0.$$

(50)

Note that if we replace the condition $a_4({\lambda }_1,{\lambda }_2)\le 0$ in (49) by the strict inequality $a_4({\lambda }_1,{\lambda }_2)<0$, QE would lead to the result $\mathtt{false}$.

The first condition of (50) is consistent with the required positive definiteness of the matrix P. The second and third conditions constrain the structure of the matrix P and thus lead to

$$P=\left(\begin{array}{ccc}{p}_{11}& 0& 0\\ 0& p_{22}& p_{23}\\ 0& p_{23}& p_{33}\end{array}\right).$$

(51)

Compared to the dense matrix P from Eq. (21) this reduces the number of free parameters from 6 to 4.

5.3. Parametrization of the Quadratic Lyapunov-like Candidate Function

Consider the matrix P with the spezial structure according to (51). The matrix P is positive definite if and only if

$$p_{11}>0\wedge p_{22}>0\wedge p_{33}>0\wedge p_{22}{p}_{33}-p_{23}^2>0.$$

(52)

With the simplification resulting from the structure (51), the numerator polynomial of (47) can be written as

$$\begin{array}{ccc}\hfill W\left(x\right)& =& -2C_{2}Lb{p}_{11}{x}_1^4-2C_{2}GL{p}_{11}{x}_1^2-2C_{2}La{p}_{11}{x}_1^2+2C_{2}GL{p}_{11}{x}_1{x}_2\hfill \\ & & +2C_{1}GL{p}_{22}{x}_1{x}_2-2C_{1}GL{p}_{22}{x}_2^2+2C_{1}GL{p}_{23}{x}_1{x}_3-2C_{1}GL{p}_{23}{x}_2{x}_3\hfill \\ & & -2C_1{C}_2{p}_{23}{x}_2^2+2C_{1}L{p}_{22}{x}_2{x}_3-2C_1{C}_2{p}_{33}{x}_2{x}_3+2C_{1}L{p}_{23}{x}_3^2.\hfill \end{array}$$

This polynomial function has no terms of degree 3. For the function W to be negative definite at infinity, the monomials $x_1^4$, $x_2^2$, $x_3^2$ must have negative coefficients:

$$\begin{array}{cc}\hfill x_1^4:& -2C_{2}Lb{p}_{11}\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill x_2^2:& -2(GL{p}_{22}+C_2{p}_{23})C_1\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill x_3^2:& \phantom{-}2C_{1}L{p}_{23}\hfill \end{array}$$

(55)

These coefficients are negative if

$$p_{11}>0\wedge p_{23}<0\wedge GL{p}_{22}+C_2{p}_{23}>0.$$

(56)

The second order mixed terms could also have an impact on the definiteness:

$$\begin{array}{cc}\hfill x_1{x}_2:& \phantom{+}2\left(C_2{p}_{11}+C_1{p}_{22}\right)GL\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill x_1{x}_3:& \phantom{+}2C_{1}GL{p}_{23}\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill x_2{x}_3:& -2\left(GL{p}_{23}-L{p}_{22}+C_2{p}_{33}\right)C_1\hfill \end{array}$$

(59)

The calculation of a bound for the delimitation of positively invariant sets would be simplified if these terms could be set to zero. Because of $p_{11}>0\wedge p_{22}>0$ according to (52) the coefficient in (57) cannot become zero. Because $p_{23}<0$ after (56), the coefficient in (58) cannot become zero either. The coefficient in (59) becomes zero if and only if

$$GL{p}_{23}-L{p}_{22}+C_2{p}_{33}=0,$$

(60)

which leads to a further simplification of W:

$$\begin{array}{ccc}\hfill W\left(x\right)& =& -2C_{2}Lb{p}_{11}{x}_1^4-2C_{2}GL{p}_{11}{x}_1^2-2C_{2}La{p}_{11}{x}_1^2\hfill \\ & & +2C_{2}GL{p}_{11}{x}_1{x}_2+2C_{1}GL{p}_{22}{x}_1{x}_2-2C_{1}GL{p}_{22}{x}_2^2\hfill \\ & & +2C_{1}GL{p}_{23}{x}_1{x}_3-2C_1{C}_2{p}_{23}{x}_2^2+2C_{1}L{p}_{23}{x}_3^2.\hfill \end{array}$$

For $C_1,C_2,G,L,a,b$ we used the parameter values given in Section 5.1. In addition, we normalized $p_{11}=1$. The conditions (52), (56), and (60) lead to

$$p_{23}<0\wedge p_{22}+10p_{23}>0\wedge p_{22}-p_{23}^2>0\wedge 70p_{33}-10p_{22}+7p_{23}=0.$$

If we add existence quantifiers for $p_{22}$ and $p_{23}$ to the expression, then after QE we obtain the condition $p_{33}>0$. For further calculation, we choose $p_{33}=1$. Eliminating $p_{23}$ yields the condition

$$p_{22}-7<0\wedge 107p_{22}-700>0.$$

on $p_{22}$. This means that $p_{22}$ lies in the (relatively small) interval

$$6.542\approx {\displaystyle \frac{700}{107}}<p_{22}<7.$$

For the parameter, we choose the arithmetic mean of the two limits:

$$p_{22}={\displaystyle \frac{1449}{214}}\approx 6.771.$$

With $p_{11},p_{22},p_{33}$ specified, we obtain the condition

$$107p_{23}+35=0$$

on the remaining parameter resulting in

$$p_{23}=-{\displaystyle \frac{35}{107}}\approx -0.327.$$

(61)

In summary, the following parameters are used for the matrix (51):

$$p_{11}=1,p_{22}={\displaystyle \frac{1449}{214}}\approx 6.77,p_{23}=-{\displaystyle \frac{35}{107}}\approx -0.327,p_{33}=1.$$

(62)

5.4. Computation of the Bounds

After specifying the Lyapunov-like candidate function V according to (8) with (51) and (62), an attempt can be made to calculate the bound $\gamma$ for a positive invariant set (3). Similar to (7) we use the prenex formula

$$\exists q\forall x_1,x_2,x_3:q>0\wedge W\left(x\right)\le -q·[V\left(x\right)-\gamma ],$$

(63)

where only the denominator W of $\dot{V}$ was used to obtain a polynomial formulation (see (47)). Applying QE to (63) should provide conditions for $\gamma$ in the form of a quantifier-free formula. Unfortunately, even after several days of computing time, the attempt with QE was unable to deliver a result, so the calculation was aborted.

However, if we replace the symbolic variable $\gamma$ with a specific positive number, the prenex formula (63) becomes a decision problem with the possible results true or false, which could be solved on a standard PC in a few hundred milliseconds. In combination with a bisection method, the bound $\gamma$ can therefore be calculated with arbitrary precision. Carrying out QE on (63) with $\gamma =2021.851$ yields true  whereas $\gamma =2021.8509$ yields false. The levels sets $\{x:{\mathbb{R}}^3:V\left(x\right)=\gamma \}$ of this bound are ellipses shown in Figure 5 with red color. In addition, the transients already depicted in Figure 4 are also shown using blue and green color for the chaotic attractor and the stable limit cycle, respectively. When specifying the elements of the matrix P, arbitrary decisions were also made to determine $p_{22}$ and $p_{33}$. Therefore, other parameter values for the admissible ranges were generated using a random number generator, and we attempted to find a bound $\gamma$ using the prenex formula (63). The levels resulting from these bounds are shown with gray color in Figure 5.

6. Discussion

Although the quantifier elimination is a powerful tool to compute definite bounds for attractors, its high computational effort restricts the applicability to rather simple systems and Lyapunov candidates. A slight simplification or reformulation of the problem can reduce the computation time significantly. Therefore, the construction of appropriate quantified formulas is an important topic. In this contribution we have shown a state transformation that is particularly useful for the proof of bounds on positive invariant sets using Lyapunov functions.

By the change of order for the quantifiers necessary conditions for the ansatz parameters can be found much more easily, that is, with less computational effort. Every such reduction for the parameter space in return simplifies the original problem. That the quantifier elimination problem must be solved for each region $D_i$, $i=1,\dots ,n$ separately is of minor importance.

From (31) to (33) or from (38) to (39) the existential quantifier for $\epsilon$ was swapped several times with the universal quantifiers for ${\lambda }_1,\dots ,{\lambda }_{n-1}$. There are a lot of intermediate steps, precisely, one for each partition of the set $\{1,2,\dots ,n-1\}$ into two, of gradually changing difficulty that could be considered: Let $A,B,C\subset \{1,\dots ,n-1\}$ be disjoint sets such that

$$A\cup B\cup C=\{1,\dots ,n-1\}.$$

Then, according to (32),

$$\underset{j\in A}{\forall }{\lambda }_j\exists \epsilon \underset{j\in B\cup C}{\forall }{\lambda }_j:F\Rightarrow \underset{j\in A\cup C}{\forall }{\lambda }_j\exists \epsilon \underset{j\in B}{\forall }{\lambda }_j:F.$$

Applying this implication to (31) (or (38)) reduces to the problem of positive definiteness (at infinity) in dimension $1+\left|B\right|$.

All these partitions and implications are shown in Figure 6 for $n=4$. For the example in Section 5 we have only used the quantified formulas with quantifiers shown on the left corresponding to (38) and on the right corresponding to (39). In general one could get additional restrictions on the parameter space for the Lyapunov candidate function by traversing the graph from the right to the left.

It should be noted that there are in general still a lot of Lyapunov functions that are positive definite (at infinity) with a negative definite (at infinity) derivative along the system dynamics for a particular parametric ansatz. This could also be observed for the example in Section 5. Each of these Lyapunov-like functions yields a positive invariant set, where one candidate function may lead to a better bound than the other. In principle, one could compute the intersection of all these bounds, which is, however, an even more sophisticated semi-algebraic problem.