Exact Asymptotic Stability Analysis and Region-of-Attraction Estimation for Nonlinear Systems

Abstract

We address the problem of asymptotic stability and region-of-attraction analysis of nonlinear dynamical systems. A hybrid symbolic-numeric method is presented to compute exact Lyapunov functions and exact estimates of regions of attraction of nonlinear systems efficiently. A numerical Lyapunov function and an estimate of region of attraction can be obtained by solving an (bilinear) SOS programming via BMI solver, then the modified Newton refinement and rational vector recovery techniques are applied to obtain exact Lyapunov functions and verified estimates of regions of attraction with rational coefficients. Experiments on some benchmarks are given to illustrate the efficiency of our algorithm.

1. Introduction

Lyapunov's stability theory, which concern is with the behavior of trajectories of differential systems, plays an important role in analysis and synthesis of nonlinear continuous systems. Lyapunov functions can be used to verify the asymptotic stability, including locally asymptotic stability and globally asymptotic stability. In the literature, there has been lots of work on constructing Lyapunov functions [1–12]. For example, [4, 5] used the linear matrix inequality (LMI) method to compute quadratic Lyapunov functions. In [8, 13], based on sum of squares (SOS) decomposition, a method is proposed to construct high-degree numerical polynomial Lyapunov functions. Reference [9] proposed a new method for computing polynomials Lyapunov functions by solving semialgebraic constraints via the tool DISCOVERER [14]. Reference [15] constructed Lyapunov functions beyond polynomials by using radial basis functions. In [11], the Gröbner based method is used to choose the parameters in Lyapunov functions in an optimal way.

Since the analysis of asymptotic stability is not sufficient for safety critical systems, the analysis of the region of attraction (ROA) of an asymptotically stable equilibrium point is a topic of significant importance. The ROA is the set of initial states from which the system converges to the equilibrium point. Computing exact regions of attraction (ROAs) for nonlinear dynamical systems is very hard if not impossible; therefore, researchers have focused on finding estimates of the actual ROAs. There are many well-established techniques for computation of ROAs [5, 16–22]. Among all methods, those based on Lyapunov functions are dominant in the literature. These methods compute both a Lyapunov function as a stability certificate and the sublevel sets of this Lyapunov function, which give rise to estimates of the ROA. In [17], the authors computed quadratic Lyapunov functions to optimize the volume of an ellipsoid contained in ROA by using nonlinear programming. Reference [5] employed LMI based method to compute the optimal quadratic Lyapunov function for maximizing the volume of an ellipsoid contained in the ROA of odd polynomial systems. In [22], based on SOS decomposition, a method is presented to search for polynomial Lyapunov functions that enlarge a provable ROA of nonlinear polynomial system. Reference [18] used discretization (in time) to flow invariant sets backwards along the flow of the vector field, obtaining larger and larger estimates for the ROA. Reference [21] employed quantifier elimination (QE) method via QEPCAD to find Lyapunov functions for estimating the ROA.

Taking advantage of the efficiency of numerical computation and the error-free property of symbolic computation, in this paper, we present a hybrid symbolic-numeric algorithm to compute exact Lyapunov functions for asymptotic stability and verified estimates of the ROAs of continuous systems. The algorithm is based on SOS relaxation [23] of a parametric polynomial optimization problem with LMI or bilinear matrix inequality (BMI) constraints, which can be solved directly using LMI or BMI solver in MATLAB, such as SOSTOOLS [24], YALMIP [25], SeDuMi [26], PENBMI [27], and exact SOS representation recovery techniques presented in [28, 29]. Unlike the numerical approaches, our method can yield exact Lyapunov functions and verified estimates of ROAs, which can overcome the unsoundness in analysis of nonlinear systems caused by numerical errors [30]. In comparison with some symbolic approaches based on qualifier elimination technique, our approach is more efficient and practical, because parametric polynomial optimization problem based on SOS relaxation with fixed size can be solved in polynomial time theoretically.

The rest of the paper is organized as follows. In Section 2, we introduce some definitions and notions about Lyapunov stability and ROA. Section 3 is devoted to transform the problem of computing Lyapunov functions and estimates of ROAs into a parametric program with LMI or BMI constraints. In Section 4, a symbolic-numeric approach via SOS relaxation and exact rational vector recovery is proposed to compute Lyapunov functions and estimates of ROAs, and an algorithm is described. In Section 5, experiments on some benchmarks are shown to illustrate our algorithm on asymptotic stability and ROA analysis. At last, we conclude our results in Section 6.

2. Lyapunov Stability and Region of Attraction

In this section, we will present the notion of Lyapunov stability and regions of attraction (ROAs) of dynamical systems.

Consider the autonomous system 𝐱̇=𝐟(𝐱), (1) where 𝐟∶ℝ^𝑛→ℝ^𝑛 is continuous and 𝐟 satisfies the Lipschitz condition. Denote by 𝜙(𝑡;𝐱₀) the solution of (1) corresponding to the initial state 𝐱₀=𝐱(0), evaluated at time 𝑡>0.

A vector 𝐱∈ℝ^𝑛 is an equilibrium point of the system (1) if 𝐟(𝐱)=𝟎. Since any equilibrium point can be shifted to the origin 𝟎 via a change of variables, without loss of generality, we may always assume that the equilibrium point of interest occurs at the origin.

Lyapunov theory is concerned with the behavior of the solution 𝜙(𝑡;𝐱₀) where the initial state 𝐱₀ is not at the equilibrium 𝟎 but is “close” to it.

Definition 1 (Lyapunov stability). The equilibrium point 𝟎 of (1) is (i)stable, if for any 𝜖>0, there exists 𝛿=𝛿(𝜖) such that 𝐱₀<𝛿⟹𝜙𝑡;𝐱₀<𝜖,∀𝑡>0, (2) here ‖⋅‖ denotes any norm defined on ℝ^𝑛; (ii)unstable, if it is not stable, (iii)asymptotically stable, if it is stable, and 𝛿 can be chosen such that 𝐱₀<𝛿⟹lim_𝑡→∞𝜙𝑡;𝐱₀=𝟎; (3)(iv)globally asymptotically stable if it is stable, and for all 𝐱₀∈ℝ^𝑛, lim_𝑡→∞𝜙𝑡;𝐱₀=𝟎. (4)

Intuitively, the equilibrium point 𝟎 is stable if all solutions starting near the origin (meaning that the initial conditions are in a neighborhood of the origin) remain near the origin for all time. The equilibrium point 𝟎 is asymptotically stable if all solutions starting at nearby points not only stay nearby but also converge to the equilibrium point as time approaches infinity. And the equilibrium point 𝟎 is globally asymptotically stable if it is asymptotically stable for all initial conditions 𝐱₀∈ℝ^𝑛. The stability is an important property in practice, because it means arbitrarily small perturbations of the initial state about the equilibrium point 𝟎 result in arbitrarily small perturbations in the corresponding solution trajectories of (1).

A sufficient condition for stability of the zero equilibrium is the existence of a Lyapunov function, as shown in the following theorem.

Theorem 2 ([31, Theorem 4.1]). Let 𝐷⊂ℝ^𝑛 be a domain containing the equilibrium point 𝟎 of (1). If there exists a continuously differentiable function 𝑉∶𝐷→ℝ such that 𝑉(𝟎)(𝐱){𝟎},𝑉̇(𝐱)=0,𝑉>0in𝐷⧵∶=𝜕𝑉(𝐱)𝜕𝐱⋅𝐟≤0in𝐷, (5)(6) then the origin is stable. Moreover, if 𝑉̇(𝐱)<0in𝐷⧵{𝟎}, (7) then the origin is asymptotically stable.

A function 𝑉(𝐱) satisfying the conditions (5) and (6) in Theorem 2 is commonly known as a Lyapunov function. And we can verify globally asymptotic stability of system (1) by using Lyapunov functions stated as follows.

Theorem 3 ([31, Theorem 4.2]). Let the origin be an equilibrium point for (1). If there exists a continuously differentiable function 𝑉∶ℝ^𝑛→ℝ such that 𝑉(𝟎)(𝐱)‖𝐱‖(𝐱)𝑉̇(𝐱)=0,𝑉>0∀𝐱≠𝟎,⟶∞⟹𝑉⟶∞,<0∀𝐱≠𝟎, (8)(9)(10) then the origin is globally asymptotically stable.

Remark that a function satisfying condition (9) is said to be radially unbounded.

It is known that globally asymptotic stability is very desirable, but in many applications it is difficult to achieve. When the equilibrium point is asymptotically stable, we are interested in determining how far from the origin the trajectory can be and still converge to the origin as 𝑡 approaches ∞. This gives rise to the following definition.

Definition 4 (region of attraction). The region of attraction (ROA) Ω of the equilibrium point 𝟎 is defined as Ω={𝐱∈ℝ^𝑛∣lim_𝑡→∞𝜙(𝑡;𝐱)=𝟎}.

The ROA of the equilibrium point 𝟎 is a collection of all points 𝐱 such that any trajectory starting at initial state 𝐱 at time 0 will be attracted to the equilibrium point. In the literature, the terms “domain of attraction” and “attraction basin” are also used instead of “region of attraction”.

Definition 5 (positively invariant set). A set 𝑀⊂ℝ^𝑛 is called a positively invariant set of the system (1), if 𝐱₀∈𝑀 implies 𝜙(𝑡;𝐱₀)∈𝑀 for all 𝑡≥0. Namely, if a solution belongs to a positively invariant set 𝑀 at some time instant, then it belongs to 𝑀 for all future time.

In general, finding the exact ROA analytically might be difficult or even impossible. Usually, Lyapunov functions are used to find underestimates of the ROA, that is, sets contained in the region of attraction, as stated in the following theorem.

Theorem 6 ([20, Theorem 10]). Let 𝐷⊂ℝ^𝑛 be a domain containing the equilibrium point 𝟎 of (1). If there exists a continuously differentiable function 𝑉∶𝐷→ℝ satisfying (5) and (7), then any region Ω_𝑐={𝐱∈ℝ^𝑛∣𝑉(𝐱)≤𝑐} with 𝑐≥0 such that Ω_𝑐⊆𝐷 is a positively invariant set contained in the ROA of the equilibrium 𝟎.

Theorem 6 describes an approach to compute estimates of the ROA through Lyapunov functions. More specifically, in the case of asymptotic stability, if Ω_𝑐={𝐱∈ℝ^𝑛∣𝑉(𝐱)≤𝑐} is bounded and contained in 𝐷, then every trajectory starting in Ω_𝑐 remains in Ω_𝑐 and approaches the origin as 𝑡→∞. Thus, Ω_𝑐 is an estimate of the ROA. Remark that when the origin is globally asymptotically stable, the region of attraction is the whole space ℝ^𝑛.

3. Problem Reformulation

In this section, we will transform the problem of the asymptotic stability and ROA analysis of system (1) to a parametric program with LMI or BMI constraints. In the sequel, we suppose that system (1) is a polynomial dynamical system with 𝑓_𝑖(𝐱)∈ℝ[𝐱] for 1≤𝑖≤𝑛, where ℝ[𝐱] denotes the ring of real polynomials in the variables 𝐱.

3.1. Asymptotic Stability

Firstly, we consider the asymptotic stability of system (1). As shown in Theorem 2, the existence of a Lyapunov function 𝑉(𝐱) which satisfies the conditions (5) and (7) is a certificate for asymptotical stability of the equilibrium point 𝟎, and the problem of computing such a 𝑉(𝐱) can be transformed into the following problem: (𝐱)[𝐱](𝟎)𝑉(𝐱){𝟎},𝑉̇(𝐱)=find𝑉∈ℝs.t.𝑉=0,>0∀𝐱∈𝐷⧵𝜕𝑉𝐟(𝐱){𝟎}.𝜕𝐱<0∀𝐱∈𝐷⧵ (11) In general, 𝐷 can be an arbitrary neighborhood of the equilibrium point 𝟎. However, to simplify calculation, in practice, we can assume 𝐷={𝐱∈ℝ^𝑛∶𝑔(𝐱)≤0}, where 𝑔(𝐱)∈ℝ[𝐱] is chosen to be 𝑔(𝐱)=∑^𝑛_𝑖=1𝑥²_𝑖−𝑟² for 1≤𝑖≤𝑛, with a given constant 𝑟∈ℝ₊, for instance, 𝑟=10⁻².

Let us first predetermine a template of Lyapunov functions with the given degree 𝑑. We assume that 𝑉(𝐱)=_𝛼𝑐_𝛼𝐱^𝛼, (12) where 𝐱^𝛼=𝑥^𝛼1₁,…,𝑥^𝛼𝑛_𝑛, 𝛼=(𝛼₁,…,𝛼_𝑛)∈ℤ^𝑛_≥0 with ∑^𝑛_𝑖=1𝛼_𝑖≤𝑑, and 𝑐_𝛼∈ℝ are parameters. We can rewrite 𝑉(𝐱)=𝐜^𝑇⋅𝑇(𝐱), where 𝑇(𝐱) is the (column) vector of all terms in 𝑥₁,…,𝑥_𝑛 with total degree ≤𝑑, and 𝐜∈ℝ^𝜈, with 𝜈=_𝑛^𝑛+𝑑, is the coefficient vector of 𝑉(𝐱). In the sequel, we write 𝑉(𝐱) as 𝑉(𝐱,𝐜) for clarity.

For a polynomial 𝜑(𝐱)∈ℝ[𝐱], we say that 𝜑(𝐱) is positive definite (resp., positive semidefinite), if 𝜑(𝐱)>0 for all 𝐱∈ℝ^𝑛⧵{𝟎} (resp., 𝜑(𝐱)≥0 for all 𝐱∈ℝ^𝑛). Observe that, if 𝜑(𝐱) is a sum of squares (SOS), then 𝜑(𝐱) is globally nonnegative. And, to ensure positive definiteness of 𝜑(𝐱), we can use a polynomial 𝑙(𝐱) of the form: 𝑙(𝐱)=^𝑛_𝑖=1𝜖_𝑖𝑥^𝑘_𝑖, (13) where 𝜖_𝑖∈ℝ₊ and 𝑘 is assumed to be even. Clearly, the condition that 𝜑(𝐱)−𝑙(𝐱) is an SOS polynomial guarantees the positive definiteness of 𝜑(𝐱). Therefore, to ensure positive definiteness of the second and third constraints in (11), we can use two polynomials 𝑙₁(𝐱), 𝑙₂(𝐱) of the form of (13). It is notable that SOS programming can be applied to determine the nonnegativity of a multivariate polynomial over a semialgebraic set. Consider the problem of verifying whether the implication ^𝑚_𝑖=1𝑝_𝑖(𝐱)≥0⟹𝑞(𝐱)≥0 (14) holds, where 𝑝_𝑖(𝐱)∈ℝ[𝐱] for 1≤𝑖≤𝑚 and 𝑞(𝐱)∈ℝ[𝐱]. According to Stengle's Positivstellensatz, Schmüdgen's Positivstellensatz, or Putinar's Positivstellensatz [32], if there exist SOS polynomials 𝜎_𝑖∈ℝ[𝐱] for 𝑖=0,…,𝑚, such that 𝑞(𝐱)=𝜎₀(𝐱)+^𝑚_𝑖=1𝜎_𝑖(𝐱)𝑝_𝑖(𝐱), (15) then the assertion (14) holds. Therefore, the existence of SOS representations provides a sufficient condition for determining the nonnegativity of 𝑞(𝐱) over {𝐱∈ℝ^𝑛∶∧^𝑚_𝑖=1𝑝_𝑖(𝐱)≥0}.

Based on the above observation, the problem (11) can be further transformed into the following SOS programming: find𝐜∈ℝ^𝜈()𝑉()s.t.𝑉𝟎,𝐜=0,𝐱,𝐜−𝑙₁(𝐱)=𝜌₁(𝐱)−𝜎₁(𝐱)𝑔(𝐱),−𝜕𝑉𝐟(𝐱)𝜕𝐱−𝑙₂(𝐱)=𝜌₂(𝐱)−𝜎₂(𝐱)𝑔(𝐱), (16) where 𝜎_𝑗(𝐱) and 𝜌_𝑗(𝐱) are SOS in ℝ[𝐱] for 𝑗=1,2. Moreover, the degree bound of those unknown SOS polynomials 𝜎_𝑗, 𝜌_𝑗 is exponential in 𝑛, deg(𝑉), deg(𝐟), and deg(𝑔). In practice, to avoid high computational complexity, we simply set up a truncated SOS programming by fixing a priori (much smaller) degree bound 2𝑒, with 𝑒∈ℤ₊, of the unknown SOS polynomials. Consequently, the existence of a solution 𝐜 of (16) can guarantee the asymptotic stability of the given system.

Similarly, the problem of computing the Lyapunov function 𝑉(𝐱) for globally asymptotic stability of system (1), that satisfies the conditions in Theorem 3 can be rewritten as the following problem: (𝐱)[𝐱](𝟎)𝑉(𝐱)‖𝐱‖(𝐱)𝑉̇(𝐱)find𝑉∈ℝs.t.𝑉=0,>0∀𝐱≠𝟎,⟶∞⟹𝑉⟶∞,<0∀𝐱≠𝟎. (17) By introducing two polynomials 𝑙_𝑖(𝐱), 𝑖=1,2 of the form (13), the condition, that 𝑉(𝐱,𝐜)−𝑙₁(𝐱) is SOS, guarantees the radially unboundedness of 𝑉(𝐱,𝐜). Consequently, the problem (17) can be further transformed into the following SOS programming: find𝐜∈ℝ^𝜈()𝑉()s.t.𝑉𝟎,𝐜=0,𝐱,𝐜−𝑙₁(𝐱)=𝜌₁(𝐱),−𝜕𝑉𝐟(𝐱)𝜕𝐱−𝑙₂(𝐱)=𝜌₂(𝐱), (18) where 𝜌_𝑗(𝐱) are SOS polynomials in ℝ[𝐱] for 𝑗=1,2. The decision variables are the coefficients of all the polynomials appearing in (18), such as 𝑉(𝐱,𝐜),𝑙_𝑖(𝐱), and 𝜌_𝑗(𝐱).

3.2. Region of Attraction

Suppose that the equilibrium point 𝟎 is asymptotically stable. In this section, we will consider how to find a large enough underestimate of the ROA. In the case where the equilibrium point 𝟎 is globally asymptotically stable, the ROA becomes the whole space ℝ^𝑛.

Our idea of computing the estimate of ROA is similar to that of the algorithm in [20, Section 4.2.2]. Suppose that 𝐷 is a semialgebraic set 𝐱∈ℝ^𝑛∣𝑉(𝐱)≤1 (19) given by a Lyapunov function 𝑉(𝐱). In order to enlarge the computed positively invariant set contained in the ROA, we define a variable sized region 𝑃_𝛽=𝐱∈ℝ^𝑛∣𝑝(𝐱)≤𝛽, (20) where 𝑝(𝐱)∈ℝ[𝐱] is a fixed and positive definite polynomial, and maximize 𝛽 subject to the constraint 𝑃_𝛽⊆𝐷.

Fix a template of 𝑉 of the form (12). The problem of finding an estimate 𝑃_𝛽 of the ROA can be transformed into the following polynomial optimization problem: max_{𝐜∈ℝ𝜈}𝛽()𝑉(){𝟎},𝑉̇()=s.t.𝑉𝟎,𝐜=0,𝐱,𝐜>0∀𝐱∈𝐷⧵𝐱,𝐜𝜕𝑉𝐟(𝐱){𝟎},𝑝(𝐱)()𝜕𝐱<0∀𝐱∈𝐷⧵≤𝛽⊨𝑉𝐱,𝐜≤1. (21) By introducing two polynomials 𝑙_𝑖(𝐱), 𝑖=1,2 of the form (13) and based on SOS relaxation, the problem (21) can be further transformed into the following SOS programming: max_{𝐜∈ℝ𝜈}𝛽()𝑉()s.t.𝑉𝟎,𝐜=0,𝐱,𝐜−𝑙₁(𝐱)=𝜌₁(𝐱)−𝜎₁(𝐱𝑉()),−)(𝐱,𝐜−1𝜕𝑉𝐟(𝐱)𝜕𝐱−𝑙₂(𝐱)=𝜌₂(𝐱)−𝜎₂(𝐱𝑉()),())(𝐱,𝐜−11−𝑉𝐱,𝐜=𝜌₃(𝐱)−𝜎₃(𝐱)𝑝(𝐱)−𝛽, (22) where 𝜎_𝑖(𝐱) and 𝜌_𝑖(𝐱) are SOS in ℝ[𝐱] for 1≤𝑖≤3. In (22), the decision variables are 𝛽 and the coefficients of all the polynomials appearing in (22), such as 𝑉(𝐱,𝐜), 𝜎_𝑖(𝐱), and 𝜌_𝑖(𝐱). Since 𝛽 and the coefficients of 𝑉(𝐱,𝐜) and 𝜎_𝑖(𝐱) are unknown, some nonlinear terms that are products of these coefficients, will occur in the second, third, and fourth constraints of (22), which yields a nonconvex bilinear matrix inequalities (BMI) problem. We will discuss in Section 4 how to handle the SOS programming (22) directly using the BMI solver or iterative method.

4. Exact Certificate of Sum of Squares Decomposition

According to Theorems 2, 3, and 6, the key of asymptotic stability analysis and ROA estimation lies in finding a real function 𝑉(𝐱) and a constant 𝛽 satisfying the desired conditions. We will present a symbolic-numeric hybrid method, based on SOS relaxation and rational vector recovery, to compute the exact polynomial 𝑉(𝐱) and constant 𝛽.

4.1. Approximate Solution from LMI or BMI Solver

In this section, we will discuss how to handle the SOS programmings (16), (18), and (22) directly using the LMI and BMI solver or iterative method.

Using Gram matrix representation method [23] (also known as square matrix representation (SMR) [33]), a polynomial 𝜓(𝐱) of degree 2𝑚 is SOS if and only if there exists a positive semidefinite matrix 𝑄 such that 𝜓(𝐱)=𝑍(𝐱)^𝑇𝑄𝑍(𝐱), (23) where 𝑍(𝐱) is a monomial vector in 𝐱 of degree 𝑚. Thus, the SOS programming (16) is equivalent to the following semidefinite program (SDP) problem: inf_{𝑊[𝑗],𝑄[𝑗]}²_𝑗=1Trace𝑊^[𝑗]+𝑄^[𝑗](𝟎)𝑉()s.t.𝑉=0,𝐱,𝐜−𝑙₁(𝐱)=𝑚₁(𝐱)^𝑇⋅𝑊^[1]⋅𝑚₁(𝐱)−𝑝₁(𝐱)^𝑇⋅𝑄^[1]⋅𝑝₁(𝐱)𝑔(𝐱),−𝜕𝑉𝐟(𝐱)𝜕𝐱−𝑙₂(𝐱)=𝑚₂(𝐱)^𝑇⋅𝑊^[2]⋅𝑚₂(𝐱)−𝑝₂(𝐱)^𝑇⋅𝑄^[2]⋅𝑝₂(𝐱)𝑔(𝐱), (24) where all the matrices 𝑊^[𝑗],𝑄^[𝑗] are symmetric and positive semidefinite, and ∑²_𝑗=1Trace(𝑊^[𝑗]+𝑄^[𝑗]) denotes the sum of traces of all these matrices, which acts as a dummy objective function commonly used in SDP for optimization problem with no objective function.

Similarly, the SOS programming (18) can be rewritten as the following SDP problem: inf_{𝑊[1],𝑊[2]}Trace𝑊^[1]+Trace𝑊^[2](𝟎)𝑉()s.t.𝑉=0,𝐱,𝐜−𝑙₁(𝐱)=𝑚₁(𝐱)^𝑇⋅𝑊^[1]⋅𝑚₁(𝐱),−𝜕𝑉𝐟(𝐱)𝜕𝐱−𝑙₂(𝐱)=𝑚₂(𝐱)^𝑇⋅𝑊^[2]⋅𝑚₂(𝐱), (25) where matrices 𝑊^[1], 𝑊^[2] are symmetric and positive semidefinite.

Many Matlab packages of SDP solvers, such as SOSTOOLS [24], YALMIP [25], and SeDuMi [26], are available to solve the problem (24) and (25) efficiently.

Now, let us consider the problem (22). The following example shows how to transform nonlinear parametric polynomial constraints into a BMI problem.

Example 7. To find a polynomial 𝜑(𝑥) satisfying 𝜑(𝑥)≥0∧−𝑥²+1≥0⊨2𝑥(𝑑𝜑/𝑑𝑥)≥0, it suffices to find 𝜑(𝑥) such that 2𝑥𝑑𝜑𝑑𝑥=𝜙₀(𝑥)+𝜙₁(𝑥)1−𝑥²+𝜙₂(𝑥)𝜑(𝑥), (26) where 𝜙₀(𝑥), 𝜙₁(𝑥), 𝜙₂(𝑥) are SOSes. Suppose that deg(𝜑)=1, deg(𝜙₀)=2 and deg(𝜙₁)=deg(𝜙₂)=0 and that 𝜑(𝑥)=𝑢₀+𝑢₁𝑥, 𝜙₁=𝑢₂ and 𝜙₂=𝑣₁, with 𝑢₀,𝑢₁,𝑢₂,𝑣₁∈ℝ parameters. From (1) we have 𝜙₀(𝑥)=𝑢₂𝑥²+2𝑢₁−𝑢₁𝑣₁𝑥−𝑢₂−𝑢₀𝑣₁, (27) whose square matrix representation (SMR) [33] is 𝜙₀(𝑥)=𝑍^𝑇𝑄𝑍, where 𝑄=−𝑢₂−𝑢₀𝑣₁𝑢₁−12𝑢₁𝑣₁𝑢₁−12𝑢₁𝑣₁𝑢₂1𝑥,𝑍=. (28) Since all the 𝜙_𝑖(𝑥) are SOSes, we have 𝑢₂≥0, 𝑣₁≥0 and 𝑄≽0. Therefore, the original constraint is translated into a BMI constraint: ℬ𝑢₀,𝑢₁,𝑢₂,𝑣₁=𝑢₁0000000000010010+𝑢₂1000000000−100001+𝑣₁0000010000000000+𝑢₀𝑣₁0000000000−100000+𝑢₁𝑣₁100000000000−2100−20≽0. (29)

Similar to Example 7, the SOS programming (22) can be transformed into a BMI problem of the form ()inf−𝛽s.t.ℬ𝐮,𝐯=𝐴₀+^𝑚_𝑖=1𝑢_𝑖𝐴_𝑖+^𝑘_𝑗=1𝑣_𝑗𝐴_𝑚+𝑗+_{1≤𝑖≤𝑚1≤𝑗≤𝑘}𝑢_𝑖𝑣_𝑗𝐵_𝑖𝑗≽0, (30) where 𝐴_𝑖, 𝐵_𝑖𝑗 are constant symmetric matrices, 𝐮=(𝑢₁,…,𝑢_𝑚), 𝐯=(𝑣₁,…,𝑣_𝑘) are parameter coefficients of the SOS polynomials occurring in the original SOS programming (22).

Many methods can be used to solve the BMI problem (30) directly, such as interior-point constrained trust region methods [34] an augmented Lagrangian strategy [35]. PENBMI solver [27] is a Matlab package to solve the general BMI problem, whose idea is based on a choice of penalty/barrier function Φ_𝑝 that penalizes the inequality constraints. This function satisfies a number of properties [35] such that, for any 𝑝>0, we have ℬ(𝐮,𝐯)≽0⟺Φ_𝑝(ℬ(𝐮,𝐯))≽0. (31) This means that, for any 𝑝>0, the problem (30) is equivalent to the following augmented problem: inf−𝛽s.t.Φ_𝑝(ℬ(𝐮,𝐯))≽0. (32) The associated Lagrangian of (32) can be viewed as a (generalized) augmented Lagrangian of (30): 𝐹𝐮,𝐯,𝑈,𝑝=−𝛽+𝑇𝑟𝑈Φ_𝑝(ℬ(𝐮,𝐯)), (33) where 𝑈 is the Lagrangian multiplier associated with the inequality constraint. Remark that 𝑈 is a real symmetric matrix of the same order as the matrix operator ℬ. For more details please refer to [35].

Alternatively, observing (30), ℬ(𝐮,𝐯) involves no crossing products like 𝑢_𝑖𝑢_𝑗 and 𝑣_𝑖𝑣_𝑗. Taking this special form into account, an iterative method can be applied by fixing 𝐮 and 𝐯 alternatively, which leads to a sequential convex LMI problem [20]. Remark that although the convergence of the iterative method cannot be guaranteed, this method is easier to implement than PENBMI solver and can yield a feasible solution (𝐮,𝐯) efficiently in practice.

4.2. Exact SOS Recovery

Since the SDP (LMI) or BMI solvers in Matlab are running in fixed precision, applying the techniques in Section 4.1 only yields numerical solutions of (16), (18), and (22). In the sequel, we will propose an improved algorithm based on a modified Newton refinement and rational vector recovery technique, to compute exact solutions of polynomial optimization problems with LMI or BMI constraints.

Without loss of generality, we can reduce the problems (16), (18), and (22) to the following problem: find𝐜∈ℝ^𝜈s.t.𝑉₁(𝐱,𝐜)=𝐦₁(𝐱)^𝑇⋅𝑊^[1]⋅𝐦₁(𝐱),𝑉₂(𝐱,𝐜)=𝐦₂(𝐱)^𝑇⋅𝑊^[2]⋅𝐦₂(𝐱)+𝐦₃(𝐱)^𝑇⋅𝑊^[3]⋅𝐦₃(𝐱)⋅𝑉₃(),𝑊𝐱,𝐜_[𝑖]≽0,𝑖=1,2,3, (34) where the coefficients of the polynomials 𝑉_𝑖(𝐱,𝐜), 1≤𝑖≤3 are affine in 𝐜. Note that (34) involves both LMI and BMI constraints. After solving the SDP system (34) by applying the techniques in Section 4.1, the numerical vector 𝐜 and the numerical positive semidefinite matrices 𝑊^[𝑖], 𝑖=1,2,3 may not satisfy the conditions in (34) exactly, that is, 𝑉₁(𝐱,𝐜)≈𝐦₁(𝐱)^𝑇⋅𝑊^[1]⋅𝐦₁(𝐱),𝑉₂(𝐱,𝐜)≈𝐦₂(𝐱)^𝑇⋅𝑊^[2]⋅𝐦₂(𝐱)+𝐦₃(𝐱)^𝑇⋅𝑊^[3]⋅𝐦₃(𝐱)⋅𝑉₃(𝐱,𝐜), (35) as illustrated by the following example.

Example 8. Consider the following nonlinear system: 𝑥̇₁𝑥̇₂=−0.5𝑥₁+𝑥₂+0.24999𝑥²₁+0.083125𝑥³₁+0.0205295𝑥⁴₁+0.0046875𝑥⁵₁+0.0015191𝑥⁶₁−𝑥₂+𝑥₁𝑥₂−0.49991𝑥³₁+0.040947𝑥⁵₁. (36) We want to find a certified estimate of the ROA. In the associated SOS programming (22) with BMI constraints, we suppose 𝑝(𝑥₁,𝑥₂)=𝑥²₁+𝑥²₂. When 𝑑=2, we obtain 𝑉𝑥₁,𝑥₂=3.112937368𝑥²₁+3.112937368𝑥²₂,𝛽=0.32124. (37) However, 𝑉(𝑥₁,𝑥₂) and 𝛽 cannot satisfy the conditions in (21) exactly, because there exists a sample point (15/32,85/256) such that the third condition in (21) cannot be satisfied. Therefore, Ω_𝑉=𝑥₁,𝑥₂∈ℝ²∣𝑉𝑥₁,𝑥₂≤1 (38) is not an estimate of the ROA of this system.

In our former papers [36, 37], we applied Gauss-Newton iteration and rational vector recovery to obtain exact solutions that satisfy the constraints in (34), exactly. However, these techniques may fail in some cases, as shown in [38, Example 2]. The reason may lie in that we recover the vector 𝐜 and the associated positive semidefinite matrices separately. Here, we will compute exact solutions of (34) by using a modified Newton refinement and rational vector recovery technique [38], which applies on the vector 𝐜 and the associated positive semidefinite matrices simultaneously. The main ideal is as follows.

We first convert 𝑊^[3] to a nearby rational positive semidefinite matrix 𝑊^∼[3] by nonnegative truncated PLDL^TP^T-decomposition. In practice, 𝑊^[3] is numerical diagonal matrix; in other words, the off-diagonal entries are very tiny and the diagonal entries are nonnegative. Therefore, by setting the small entries of 𝑊^[3] to zeros, we easily get the nearby rational positive semidefinite matrix 𝑊^∼[3]. We then apply Gauss-Newton iteration to refine 𝐜, 𝑊^[1], and 𝑊^[2] simultaneously with respect to a given tolerance 𝜏.

Now, we will discuss how to recover from the refined 𝐜,𝑊^[1],𝑊^[2], the rational vector 𝐜^∼, and rational positive semidefinite matrices 𝑊^∼[1] and 𝑊^∼[2] that satisfy exactly𝑉₁𝐱,𝐜^∼−𝐦₁(𝐱)^𝑇⋅𝑊^∼[1]⋅𝐦₁(𝐱)𝑉=0,₂𝐱,𝐜^∼−𝐦₂(𝐱)^𝑇⋅𝑊^∼[2]⋅𝐦₂(𝐱)−𝑉₃𝐱,𝐜^∼𝐦₃(𝐱)^𝑇⋅𝑊^∼[3]⋅𝐦₃(𝐱)=0. (39) Since the equations in (39) are affine in entries of 𝐜^∼ and 𝑊^∼₁, 𝑊^∼₂, one can define an affine linear hyperplane ℒ=𝐜,𝑊^[1],𝑊^[2]∣𝑉₁(𝐱,𝐜)−𝐦₁(𝐱)^𝑇⋅𝑊^[1]⋅𝐦₁(𝐱)𝑉=0,₂(𝐱,𝐜)−𝐦₂(𝐱)^𝑇⋅𝑊^[2]⋅𝐦₂(𝐱)−𝑉₃(𝐱,𝐜)𝐦₃(𝐱)^𝑇⋅𝑊^∼[3]⋅𝐦₃(𝐱)=0. (40) Note that the hyperplane (40) can be constructed from a linear system 𝐴𝐲=𝐛, where 𝐲 consists of the entries of 𝐜, 𝑊^[1] and, 𝑊^[2]. If 𝐴 has full row rank, such a hyperplane is guaranteed to exist. Then, for a given bound 𝐷 of the common denominator, the rationalized SOS solutions of (34) can be computed by orthogonal projection if the matrices 𝑊^[1],𝑊^[2] have full ranks with respect to 𝜏, or by the rational vector recovery otherwise. Finally, we check whether the matrices 𝑊^∼[𝑖] are positive semidefinite for 𝑖=1,2. If so, the rational vector 𝐜^∼ and rational positive semidefinite matrices 𝑊^∼[𝑖]1≤𝑖≤3 can satisfy the conditions of problem (34) exactly. For more details the reader can refer to [38].

4.3. Algorithm

The results in Sections 4.1 and 4.2 yield an algorithm to find exact solutions to the problem (34).

Algorithm 9. Verified Parametric Optimization Solver.
Input(i)A polynomial optimization problem of the form (34).(ii)𝐷∈ℤ_>0: the bound of the common denominator.(iii)𝑒∈ℤ_≥0: the degree bound 2𝑒 of the SOSes used to construct the SOS programming. (iv)𝜏∈ℝ_>0: the given tolerance.
Output(i)The verified solution 𝐜^∼ of (34) with the 𝑊^∼[𝑖],1≤𝑖≤3 positive semidefinite.(1)(Compute the numerical solutions) Apply LMI or BMI solver to compute numerical solutions of the associated polynomial optimization problem (34). If this problem has no feasible solutions, return “we can't find solutions of (34) with the given degree bound 2𝑒.” Otherwise, obtain 𝐜 and 𝑊^[𝑖]⪸0,1≤𝑖≤3.(2)(Compute the verified solution 𝐜^∼) (2.1) Convert 𝑊^[3] to a nearby rational positive semidefinite matrix 𝑊^∼[3] by non-negative truncated PLDL^TP^T-decomposition. (2.2) For the tolerance 𝜏, apply the modified Newton iteration to refine 𝐜 and 𝑊^[1], 𝑊^[2]. (2.3) Determine the singularity of 𝑊^[1] and 𝑊^[2] with respect to 𝜏. Then, for a given common denominator, the rational vector 𝐜^∼, and rational matrices 𝑊^∼[1], 𝑊^∼[2] can be obtained by orthogonal projection if 𝑊^[1], 𝑊^[2] are of full rank, or by rational vector recovery method otherwise. (2.4) Check whether the matrices 𝑊^∼[1], 𝑊^∼[2] are positive semidefinite. If so, return 𝐜^∼, and 𝑊^∼[𝑖], 1≤𝑖≤3. Otherwise, return “we can't find the solutions of (34) with the given degree bound”.

5. Experiments

Let us present some examples of asymptotic stability and ROAs analysis of nonlinear systems.

Example 10 ([8, Example 1]). Consider a nonlinear continuous system ̇(𝐱)=𝑥̇𝐱=𝐟¹𝑥̇₂𝑥̇₃𝑥̇₄𝑥̇₅𝑥̇₆=−𝑥³₁+4𝑥³₂−6𝑥₃𝑥₄−𝑥₁−𝑥₂+𝑥³₅𝑥₁𝑥₄−𝑥₃+𝑥₄𝑥₆𝑥₁𝑥₃+𝑥₃𝑥₆−𝑥³₄−2𝑥³₂−𝑥₅+𝑥₆−3𝑥₃𝑥₄−𝑥³₅−𝑥₆. (41) Firstly, we will certify the locally asymptotic stability of this system. According to Theorem 2, we need to find a Lyapunov function 𝑉(𝐱), which satisfies all the conditions in Theorem 2. We can set up an associated SOS programming (16) with 𝑔(𝐱)=𝑥²₁+𝑥²₂+𝑥²₃−0.0001. When 𝑑=4, we obtain a feasible solution of the associated SDP system. Here we just list one approximate polynomial 𝑉(𝐱)=0.81569𝑥²₁+0.18066𝑥²₂+1.6479𝑥²₃+2.1429𝑥²₄+1.2996𝑥²₆+⋯+0.67105𝑥⁴₄+0.52855𝑥⁴₅. (42) Let the tolerance 𝜏=10⁻², and let the bound of the common denominator of the polynomial coefficients vector be 100. By use of the rational SOS recovery technique described in Section 4.2, we can obtain an exact Lyapunov function and the corresponding SOS polynomials. Here we only list the Lyapunov function: 𝑉^∼(𝐱)=41𝑥50²¹+9𝑥50²²+33𝑥20²³+107𝑥50²⁴+69𝑥100²⁵+13𝑥10²⁶+⋯+67𝑥100⁴⁴+53𝑥100⁴⁵. (43) Therefore, the locally asymptotic stability is certified.
Furthermore, we will consider the globally asymptotic stability. It suffices to find a Lyapunov function 𝑉(𝐱) with rational coefficients, which satisfies all the conditions in Theorem 3. By solving the associated SOS programming (18), we obtain 𝑉(𝐱)=0.78997𝑥²₁+0.19𝑥²₂+2.2721𝑥²₃+1.4213𝑥²₆+⋯+1.4213𝑥⁴₂+0.71066𝑥⁴₅. (44) By use of the rational SOS recovery technique described in Section 4.2, we then obtain 𝑉^∼(𝐱)=79𝑥100²¹+19𝑥50¹𝑥₂+9𝑥100¹𝑥₅+19𝑥100²²+227𝑥100²³+199𝑥100²⁴+9𝑥100²𝑥₅+71𝑥100²⁵+71𝑥50²⁶+71𝑥50⁴²+71𝑥100⁴⁵, (45) which exactly satisfies the conditions in Theorem 3. Therefore, the globally asymptotic stability of this system is certified.

Table 1 shows the performance of Algorithm 9 on another six examples, for globally asymptotic stability analysis of dynamical systems in the literature. All the computations have been performed on an Intel Core 2 Duo 2.0 GHz processor with 2 GB of memory. In Table 1 Examples  1–3 correspond to [9, Examples 2, 3, 9] and Examples  4–6 correspond to [39, Example 7], [6, Example 22], and [40, page 1341]. For all these examples, let the degree bound of the SOSes be 4, and set 𝜏=10⁻², 𝐷=100. Here 𝑛 and Num denote the number of system variables and the number of decision variables in the LMI problem, respectively; deg(𝑉^∼) denotes the degree of 𝑉^∼(𝐱) obtained from Algorithm 9; Time is that for the entire computation run Algorithm 9 in seconds.

Table 1 

Algorithm performance on benchmarks (globally asymptotic stability).

Example 11 ([41, Example A]). Consider a nonlinear continuous system ̇(𝐱)=𝑥̇𝐱=𝐟¹𝑥̇₂=−𝑥₂𝑥₁+𝑥²₁−1𝑥₂. (46) Firstly, we will certify the locally asymptotic stability of this system. It suffices to find a Lyapunov function 𝑉(𝑥₁,𝑥₂) with rational coefficients, which satisfies all the conditions in Theorem 2. We can set up an associated SOS programming (16) with 𝑔(𝐱)=𝑥²₁+𝑥²₂−0.0001. When 𝑑=2, we obtain 𝑉(𝐱)=1.4957𝑥²₁−0.7279𝑥₁𝑥₂+1.1418𝑥²₂. (47) Let the tolerance 𝜏=10⁻², and, the bound of the common denominator of the polynomial coefficients vector be 100. By use of the rational SOS recovery technique described in Section 4.2, we then obtain a Lyapunov function: 𝑉^∼(𝐱)=−62𝑥85¹𝑥₂+127𝑥85²¹+97𝑥85²². (48) Therefore, the locally asymptotic stability is certified.
We now construct Lyapunov functions to find certified estimates of the ROA. In the associated SOS programming (22) with BMI constraints, we suppose 𝑝(𝑥₁,𝑥₂)=𝑥²₁+𝑥²₂. When 𝑑=2, we obtain 𝑉𝑥₁,𝑥₂=0.6174455𝑥²₁−0.40292𝑥₁𝑥₂+0.43078𝑥²₂,𝛽=1.3402225. (49) Let the tolerance 𝜏=10⁻⁵, and the bound of the common denominator of the polynomial coefficients vector be 10⁵. By use of the rational SOS recovery technique described in Section 4.2, we obtain 𝑉^∼𝑥₁,𝑥₂=6732𝑥10903²¹−4393𝑥10903¹𝑥₂+42271𝑥10903²²_∼,𝛽=6701,5000 (50) which exactly satisfy the conditions in Theorem 6. Therefore, Ω^∼_𝑉=𝑥₁,𝑥₂∈ℝ²∣𝑉^∼𝑥₁,𝑥₂≤1 (51) is a certified estimate of the ROA of the given system.

Table 2 shows the performance of Algorithm 9 on another three examples, for computing verified estimates of ROAs of dynamical systems with the same fixed positive definite polynomial 𝑝(𝐱) given in the literature. All the computations have been performed on an Intel Core 2 Duo 2.0 GHz processor with 2 GB of memory. In Table 2, Examples  1–3 correspond to [22, Examples 1, 2] and [20, page 75]. For all these examples, let the degree bound of the SOSes be 4, 𝜏=10⁻⁴ and 𝐷=10000. Here 𝑛 and Num denote the number of system variables and the number of decision variables in the BMI problem, respectively; deg(𝑉^∼) and 𝛽^∼ denote, respectively, the degree of 𝑉^∼(𝐱) and value of 𝛽 obtained from Algorithm 9, whereas 𝛽 is reported in the literature; time is that for the entire computation run Algorithm 9 in seconds.

Table 2 

Algorithm performance on benchmarks (ROA).

6. Conclusion

In this paper, we present a symbolic-numeric method on asymptotic stability and ROA analysis of nonlinear dynamical systems. A numerical Lyapunov function and an estimate of ROA can be obtained by solving an (bilinear) SOS programming via BMI solver. Then a method based on modified Newton iteration and rational vector recovery techniques is deployed to obtain exact polynomial Lyapunov functions and verified estimates of ROAs with rational coefficients. Some experimental results are given to show the efficiency of our algorithm. For future work, we will consider the problem of stability region analysis of nonpolynomial systems by applying a rigorous polynomial approximate technique to compute an uncertain polynomial system, whose set of trajectories contains that of the given nonpolynomial system.

Acknowledgments

This material is supported in part by the National Natural Science Foundation of China under Grants 91118007 and 61021004 (Wu Yang), the Fundamental Research Funds for the Central Universities under Grant 78210043 (Wu Yang), and the Education Department of Zhejiang Province Project under Grant Y201120383 (Lin).