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

and Applied Analysis 3 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 t approaches∞. This gives rise to the following definition. Definition 4 (region of attraction). The region of attraction (ROA) Ω of the equilibrium point 0 is defined as Ω = {x ∈ Rn | lim t→∞ φ(t; x) = 0}. The ROA of the equilibrium point 0 is a collection of all points x such that any trajectory starting at initial state x 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 setM ⊂ Rn is called a positively invariant set of the system (1), if x 0 ∈ M implies φ(t; x 0 ) ∈ M for all t ≥ 0. Namely, if a solution belongs to a positively invariant setM at some time instant, then it belongs toM 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 D ⊂ Rn be a domain containing the equilibrium point 0 of (1). If there exists a continuously differentiable function V : D → R satisfying (5) and (7), then any region Ω c = {x ∈ Rn | V(x) ≤ c} with c ≥ 0 such that Ω c ⊆ D is a positively invariant set contained in the ROA of the equilibrium 0. Theorem 6 describes an approach to compute estimates of the ROA through Lyapunov functions.More specifically, in the case of asymptotic stability, ifΩ c = {x ∈ Rn | V(x) ≤ c} is bounded and contained inD, then every trajectory starting in Ω c remains inΩ c and approaches the origin as t → ∞.Thus, Ω c is an estimate of the ROA. Remark that when the origin is globally asymptotically stable, the region of attraction is the whole space R. 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 f i (x) ∈ R[x] for 1 ≤ i ≤ n, where R[x] denotes the ring of real polynomials in the variables x. 3.1. Asymptotic Stability. Firstly, we consider the asymptotic stability of system (1). As shown in Theorem 2, the existence of a Lyapunov function V(x) which satisfies the conditions (5) and (7) is a certificate for asymptotical stability of the equilibrium point 0, and the problem of computing such a V(x) can be transformed into the following problem: find V (x) ∈ R [x] s.t. V (0) = 0, V (x) > 0 ∀x ∈ D \ {0} ,


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][2][3][4][5][6][7][8][9][10][11][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 wellestablished techniques for computation of ROAs [5,[16][17][18][19][20][21][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 MAT-LAB, 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.

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 where f : R  → R  is continuous and f satisfies the Lipschitz condition.Denote by (; x 0 ) the solution of (1) corresponding to the initial state x 0 = x(0), evaluated at time  > 0.
A vector x ∈ R  is an equilibrium point of the system (1) if f(x) = 0. Since any equilibrium point can be shifted to the origin 0 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 (; x 0 ) where the initial state x 0 is not at the equilibrium 0 but is "close" to it.
Definition 1 (Lyapunov stability).The equilibrium point 0 of ( 1) is (i) stable, if for any  > 0, there exists  = () such that     x 0     <  ⇒      (; x 0 )     < , ∀ > 0, here ‖ ⋅ ‖ denotes any norm defined on R  ; (ii) unstable, if it is not stable, (iii) asymptotically stable, if it is stable, and  can be chosen such that (iv) globally asymptotically stable if it is stable, and for all Intuitively, the equilibrium point 0 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 0 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 0 is globally asymptotically stable if it is asymptotically stable for all initial conditions x 0 ∈ R  .The stability is an important property in practice, because it means arbitrarily small perturbations of the initial state about the equilibrium point 0 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  ⊂ R  be a domain containing the equilibrium point 0 of (1).If there exists a continuously differentiable function  :  → R such that then the origin is stable.Moreover, if then the origin is asymptotically stable.
A function (x) 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  : R  → R such that then the origin is globally asymptotically stable.
Abstract and Applied Analysis 3 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 0 is defined as The ROA of the equilibrium point 0 is a collection of all points x such that any trajectory starting at initial state x 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  ⊂ R  is called a positively invariant set of the system (1), if x 0 ∈  implies (; x 0 ) ∈  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 describes an approach to compute estimates of the ROA through Lyapunov functions.More specifically, in the case of asymptotic stability, if Ω  = {x ∈ R  | (x) ≤ } 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 R  .

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   (x) ∈ R[x] for 1 ≤  ≤ , where R[x] denotes the ring of real polynomials in the variables x.

Asymptotic Stability.
Firstly, we consider the asymptotic stability of system (1).As shown in Theorem 2, the existence of a Lyapunov function (x) which satisfies the conditions ( 5) and ( 7) is a certificate for asymptotical stability of the equilibrium point 0, and the problem of computing such a (x) can be transformed into the following problem: In general,  can be an arbitrary neighborhood of the equilibrium point 0. However, to simplify calculation, in practice, we can assume  = {x ∈ R  : (x) ≤ 0}, where Let us first predetermine a template of Lyapunov functions with the given degree .We assume that where =1   ≤ , and   ∈ R are parameters.We can rewrite (x) = c  ⋅ (x), where (x) is the (column) vector of all terms in  1 , . . .,   with total degree ≤ , and c ∈ R ] , with ] = ( +  ), is the coefficient vector of (x).In the sequel, we write (x) as (x, c) for clarity.
For a polynomial (x) ∈ R[x], we say that (x) is positive definite (resp., positive semidefinite), if (x) > 0 for all x ∈ R  \ {0} (resp., (x) ≥ 0 for all x ∈ R  ).Observe that, if (x) is a sum of squares (SOS), then (x) is globally nonnegative.And, to ensure positive definiteness of (x), we can use a polynomial (x) of the form: where   ∈ R + and  is assumed to be even.Clearly, the condition that (x) − (x) is an SOS polynomial guarantees the positive definiteness of (x).Therefore, to ensure positive definiteness of the second and third constraints in (11), we can use two polynomials  1 (x),  2 (x) 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 holds, where According to Stengle's Positivstellensatz, Schmüdgen's Positivstellensatz, or Putinar's Positivstellensatz [32], if there exist SOS polynomials   ∈ R[x] for  = 0, . . ., , such that then the assertion ( 14) holds.Therefore, the existence of SOS representations provides a sufficient condition for determining the nonnegativity of (x) over {x ∈ R  : ∧  =1   (x) ≥ 0}.Based on the above observation, the problem (11) can be further transformed into the following SOS programming: where   (x) and   (x) are SOS in R[x] for  = 1, 2.Moreover, the degree bound of those unknown SOS polynomials   ,   is exponential in , deg(), deg(f), 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  ∈ Z + , of the unknown SOS polynomials.Consequently, the existence of a solution c of ( 16) can guarantee the asymptotic stability of the given system.
Similarly, the problem of computing the Lyapunov function (x) for globally asymptotic stability of system (1), that satisfies the conditions in Theorem 3 can be rewritten as the following problem: By introducing two polynomials   (x),  = 1, 2 of the form (13), the condition, that (x, c) −  1 (x) is SOS, guarantees the radially unboundedness of (x, c).Consequently, the problem (17) can be further transformed into the following SOS programming: where   (x) are SOS polynomials in R[x] for  = 1, 2. The decision variables are the coefficients of all the polynomials appearing in (18), such as (x, c),   (x), and   (x).

Region of Attraction.
Suppose that the equilibrium point 0 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 0 is globally asymptotically stable, the ROA becomes the whole space R  .
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 given by a Lyapunov function (x).In order to enlarge the computed positively invariant set contained in the ROA, we define a variable sized region where (x) ∈ R[x] 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: By introducing two polynomials   (x),  = 1, 2 of the form (13) and based on SOS relaxation, the problem ( 21) can be further transformed into the following SOS programming: where   (x) and   (x) are SOS in R[x] for 1 ≤  ≤ 3.
In (22), the decision variables are  and the coefficients of all the polynomials appearing in (22), such as (x, c),   (x), and   (x).Since  and the coefficients of (x, c) and   (x) 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.

Exact Certificate of Sum of Squares Decomposition
function (x) 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 (x) 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 (x) of degree 2 is SOS if and only if there exists a positive semidefinite matrix  such that where (x) is a monomial vector in x of degree .Thus, the SOS programming ( 16) is equivalent to the following semidefinite program (SDP) problem: inf Trace ( [] +  [] ) where all the matrices  [] ,  [] are symmetric and positive semidefinite, and ∑ 2 =1 Trace( [] +  [] ) 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.
Many Matlab packages of SDP solvers, such as SOS-TOOLS [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.
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 This means that, for any  > 0, the problem ( 30) is equivalent to the following augmented problem: The associated Lagrangian of ( 32) can be viewed as a (generalized) augmented Lagrangian of (30): 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 B. For more details please refer to [35].Alternatively, observing (30), B(u, k) involves no crossing products like     and V  V  .Taking this special form into account, an iterative method can be applied by fixing u and k 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 (u, k) efficiently in practice.

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: where the coefficients of the polynomials   (x, c), 1 ≤  ≤ 3 are affine in c.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 c and the numerical positive semidefinite matrices  [] ,  = 1, 2, 3 may not satisfy the conditions in (34) exactly, that is, as illustrated by the following example.
Example 8. Consider the following nonlinear system: We want to find a certified estimate of the ROA.In the associated SOS programming (22) with BMI constraints, we suppose However, ( 1 ,  2 ) 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, 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 c 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 c and the associated positive semidefinite matrices simultaneously.The main ideal is as follows.
We first convert  [3] to a nearby rational positive semidefinite matrix W[3] by nonnegative truncated PLDL T P Tdecomposition.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 W [3] .We then apply Gauss-Newton iteration to refine c,  [1] , and  [2]  simultaneously with respect to a given tolerance .Now, we will discuss how to recover from the refined c,  [1] ,  [2] , the rational vector c, and rational positive semidefinite matrices W [1] and W [2] that satisfy exactly Since the equations in (39) are affine in entries of c and W1 , W2 , one can define an affine linear hyperplane L = {c,  [1] ,  [2] |  1 (x, c) − m 1 (x)  ⋅  [1] ⋅ m 1 (x) = 0, Note that the hyperplane (40) can be constructed from a linear system y = b, where y consists of the entries of c,  [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 W[] are positive semidefinite for  = 1, 2. If so, the rational vector c and rational positive semidefinite matrices W[] 1 ≤  ≤ 3 can satisfy the conditions of problem (34) exactly.For more details the reader can refer to [38].

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)  ∈ Z >0 : the bound of the common denominator.
(iii)  ∈ Z ≥0 : the degree bound 2 of the SOSes used to construct the SOS programming.
(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 c and  [] ⪸ 0, 1 ≤  ≤ 3.
(2) (Compute the verified solution c) (2.1) Convert  [3] to a nearby rational positive semidefinite matrix W[3] by non-negative truncated PLDL T P T -decomposition.(2.2) For the tolerance , apply the modified Newton iteration to refine c and  [1] ,  [2] .(2.3) Determine the singularity of  [1] and  [2]  with respect to .Then, for a given common denominator, the rational vector c, and rational matrices W [1] , W [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 W [1] , W [2] are positive semidefinite.If so, return c, and Otherwise, return "we can't find the solutions of (34) with the given degree bound".

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 Firstly, we will certify the locally asymptotic stability of this system.According to Theorem 2, we need to find a Lyapunov function (x), which satisfies all the conditions in Theorem 2.
We can set up an associated SOS programming (16) with (x) = Therefore, the locally asymptotic stability is certified.Furthermore, we will consider the globally asymptotic stability.It suffices to find a Lyapunov function (x) with rational coefficients, which satisfies all the conditions in Theorem 3. By solving the associated SOS programming (18), we obtain  (x) = 0.78997 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 Firstly, we will certify the locally asymptotic stability of this system.It suffices to find a Lyapunov function ( 1 ,  2 ) with rational coefficients, which satisfies all the conditions in Theorem 2. We can set up an associated SOS programming (16) with (x) =  2 1 +  2 2 − 0.0001.When  = 2, we obtain 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 ( 1 ,  2 ) =  2  1 +  2 2 .When  = 2, we obtain  ( 1 ,  2 ) = 0.6174455 2 1 − 0.40292 1  2 + 0.43078 2 2 ,  = 1.3402225. (49) Let the tolerance  = 10 −5 , and the bound of the common denominator of the polynomial coefficients vector be 10 which exactly satisfy the conditions in Theorem 6.Therefore, 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 (x) 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 −4 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 Ṽ(x) 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.

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.
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 −2 ,  = 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 Ṽ(x) obtained from Algorithm 9; Time is that for the entire computation run Algorithm 9 in seconds.Example 11 ([41, Example A]).Consider a nonlinear continuous system ẋ

Table 2 :
Algorithm performance on benchmarks (ROA).Let the tolerance  = 10 −2 , 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: 5. By use of the rational SOS recovery technique described in Section 4.2, we obtain