Abstract and Applied Analysis

Volume 2013 (2013), Article ID 146137, 10 pages

http://dx.doi.org/10.1155/2013/146137

Research Article

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

Min Wu,¹ Zhengfeng Yang,¹ and Wang Lin²

¹Shanghai Key Laboratory of Trustworthy Computing, East China Normal University, Shanghai 200062, China

²College of Mathematics and Information Science, Wenzhou University, Wenzhou Zhejiang 325035, China

Received 14 November 2012; Revised 11 February 2013; Accepted 20 February 2013

Academic Editor: Fabio M. Camilli

Copyright © 2013 Min Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

where is continuous and satisfies the Lipschitz condition. Denote by the solution of (1) corresponding to the initial state , evaluated at time .

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 , there exists such that

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

(iv)globally asymptotically stable if it is stable, and for all ,

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

then the origin is stable. Moreover, if

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

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 .

The ROA of the equilibrium point is a collection of all points such that any trajectory starting at initial state at time 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 . 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 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 , 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:

In general, can be an arbitrary neighborhood of the equilibrium point . However, to simplify calculation, in practice, we can assume , where is chosen to be for , with a given constant , for instance, .

Let us first predetermine a template of Lyapunov functions with the given degree . We assume that

where , with , 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 for all (resp., 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:

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

holds, where for and . According to Stengle's Positivstellensatz, Schmüdgen's Positivstellensatz, or Putinar's Positivstellensatz [32], if there exist SOS polynomials for , such that

then the assertion (14) holds. Therefore, the existence of SOS representations provides a sufficient condition for determining the nonnegativity of over .

Based on the above observation, the problem (11) can be further transformed into the following SOS programming:

where and are SOS in for . Moreover, the degree bound of those unknown SOS polynomials , is exponential in , , , and . In practice, to avoid high computational complexity, we simply set up a truncated SOS programming by fixing a priori (much smaller) degree bound , 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:

By introducing two polynomials , 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:

where are SOS polynomials in for . 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

given by a Lyapunov function . In order to enlarge the computed positively invariant set contained in the ROA, we define a variable sized region

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:

By introducing two polynomials , of the form (13) and based on SOS relaxation, the problem (21) can be further transformed into the following SOS programming:

where and are SOS in for . 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 is SOS if and only if there exists a positive semidefinite matrix such that

where is a monomial vector in of degree . Thus, the SOS programming (16) is equivalent to the following semidefinite program (SDP) problem:

where all the matrices are symmetric and positive semidefinite, and 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:

where matrices , 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 , it suffices to find such that

where , , are SOSes. Suppose that , and and that , and , with parameters. From (1) we have

whose square matrix representation (SMR) [33] is , where

Since all the are SOSes, we have , and . Therefore, the original constraint is translated into a BMI constraint:

Similar to Example 7, the SOS programming (22) can be transformed into a BMI problem of the form

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 , we have

This means that, for any , 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 . 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:

where the coefficients of the polynomials , 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 , 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 . When , we obtain

However, and cannot satisfy the conditions in (21) exactly, because there exists a sample point 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 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 to a nearby rational positive semidefinite matrix by nonnegative truncated -decomposition. In practice, 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 to zeros, we easily get the nearby rational positive semidefinite matrix . We then apply Gauss-Newton iteration to refine , , and simultaneously with respect to a given tolerance .

Now, we will discuss how to recover from the refined , the rational vector , and rational positive semidefinite matrices and that satisfy exactly

Since the equations in (39) are affine in entries of and , , one can define an affine linear hyperplane

Note that the hyperplane (40) can be constructed from a linear system , where consists of the entries of , and, . 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 have full ranks with respect to , or by the rational vector recovery otherwise. Finally, we check whether the matrices are positive semidefinite for . If so, the rational vector and rational positive semidefinite matrices 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): the bound of the common denominator.(iii): the degree bound of the SOSes used to construct the SOS programming. (iv): the given tolerance.

Output

(i)The verified solution of (34) with the 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 .” Otherwise, obtain and .(2)(Compute the verified solution )

(2.1) Convert to a nearby rational positive semidefinite matrix by non-negative truncated -decomposition. (2.2) For the tolerance , apply the modified Newton iteration to refine and , . (2.3) Determine the singularity of and with respect to . Then, for a given common denominator, the rational vector , and rational matrices , can be obtained by orthogonal projection if , are of full rank, or by rational vector recovery method otherwise. (2.4) Check whether the matrices , are positive semidefinite. If so, return , and , . 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

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 . When , we obtain a feasible solution of the associated SDP system. Here we just list one approximate polynomial

Let the tolerance , 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:

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

By use of the rational SOS recovery technique described in Section 4.2, we then obtain

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 , . Here and Num denote the number of system variables and the number of decision variables in the LMI problem, respectively; 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

Num
deg
Time

(s)

1
2
11

1.391

2
2
11

1.418

3
3
18

1.674

4
4
81

2.951

5
2
24

1.832

6
2
24

1.830

Example 11 ([41, Example A]). Consider a nonlinear continuous system

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 . When , we obtain

Let the tolerance , 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:

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 , we obtain

Let the tolerance , and the bound of the common denominator of the polynomial coefficients vector be . By use of the rational SOS recovery technique described in Section 4.2, we obtain

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 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, and . Here and Num denote the number of system variables and the number of decision variables in the BMI problem, respectively; 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).

Example

Num
deg

β
Time (s)

1
2
28

19.417

2
3
46

26.635

3
2
28

9.915

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).

References

A. Polański, “Lyapunov function construction by linear programming,” Transactions on Automatic Control, vol. 42, no. 7, pp. 1013–1016, 1997.

T. A. Johansen, “Computation of Lyapunov functions for smooth nonlinear systems using convex optimization,” Automatica, vol. 36, no. 11, pp. 1617–1626, 2000.

L. Grüne and F. Wirth, “Computing control Lyapunov functions via a Zubov type algorithm,” in Proceedings of the 39th IEEE Confernce on Decision and Control, vol. 3, pp. 2129–2134, IEEE, December 2000.

R. C. L. F. Oliveira and P. L. D. Peres, “LMI conditions for robust stability analysis based on polynomially parameter-dependent Lyapunov functions,” Systems & Control Letters, vol. 55, no. 1, pp. 52–61, 2006.

G. Chesi, A. Garulli, A. Tesi, and A. Vicino, “LMI-based computation of optimal quadratic Lyapunov functions for odd polynomial systems,” International Journal of Robust and Nonlinear Control, vol. 15, no. 1, pp. 35–49, 2005.

J. Liu, N. Zhan, and H. Zhao, “Automatically discovering relaxed Lyapunov functions for polynomial dynamical systems,” Mathematics in Computer Science, vol. 6, no. 4, pp. 395–408, 2012.

T. Nguyen, Y. Mori, T. Mori, and Y. Kuroe, “QE approach to common Lyapunov function problem,” Journal of Japan Society for Symbolic and Algebraic Computation, vol. 10, no. 1, pp. 52–62, 2003.

A. Papachristodoulou and S. Prajna, “On the construction of Lyapunov functions using the sum of squares decomposition,” in Proceedings of the 41st IEEE Conference on Decision and Control, vol. 3, pp. 3482–3487, IEEE, December 2002.

Z. She, B. Xia, R. Xiao, and Z. Zheng, “A semi-algebraic approach for asymptotic stability analysis,” Nonlinear Analysis: Hybrid Systems, vol. 3, no. 4, pp. 588–596, 2009.

B. Grosman and D. R. Lewin, “Lyapunov-based stability analysis automated by genetic programming,” Automatica, vol. 45, no. 1, pp. 252–256, 2009.

K. Forsman, “Construction of Lyapunov functions using Grobner bases,” in Proceedings of the 30th IEEE Conference on Decision and Control, pp. 798–799, IEEE, December 1991.

Z. She, B. Xue, and Z. Zheng, “Algebraic analysis on asymptotic stability of continuous dynamical systems,” in Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation, pp. 313–320, ACM, 2011.

S. Prajna, P. A. Parrilo, and A. Rantzer, “Nonlinear control synthesis by convex optimization,” Transactions on Automatic Control, vol. 49, no. 2, pp. 310–314, 2004.

B. Xia, “DISCOVERER: a tool for solving semi-algebraic systems,” ACM Communications in Computer Algebra, vol. 41, no. 3, pp. 102–103, 2007.

P. Giesl, “Construction of a local and global Lyapunov function using radial basis functions,” IMA Journal of Applied Mathematics, vol. 73, no. 5, pp. 782–802, 2008.

H.-D. Chiang and J. S. Thorp, “Stability regions of nonlinear dynamical systems: a constructive methodology,” Transactions on Automatic Control, vol. 34, no. 12, pp. 1229–1241, 1989.

E. J. Davison and E. M. Kurak, “A computational method for determining quadratic Lyapunov functions for non-linear systems,” vol. 7, no. 5, pp. 627–636, 1971.

E. Cruck, R. Moitie, and N. Seube, “Estimation of basins of attraction for uncertain systems with affine and Lipschitz dynamics,” Dynamics and Control, vol. 11, no. 3, pp. 211–227, 2001.

R. Genesio, M. Tartaglia, and A. Vicino, “On the estimation of asymptotic stability regions: state of the art and new proposals,” Transactions on Automatic Control, vol. 30, no. 8, pp. 747–755, 1985.

Z. Jarvis-Wloszek, Lyapunov based analysis and controller synthesis for polynomial systems using sum-of-squares optimization [Ph.D. thesis], University of California, Berkeley, Calif, USA, 2003.

S. Prakash, J. Vanualailai, and T. Soma, “Obtaining approximate region of asymptotic stability by computer algebra: a case study,” The South Pacific Journal of Natural and Applied Sciences, vol. 20, no. 1, pp. 56–61, 2002.

T. Weehong and A. Packard, “Stability region analysis using sum of squares programming,” in Proceedings of the American Control Conference (ACC '06), pp. 2297–2302, IEEE, June 2006.

P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic problems,” Mathematical Programming B, vol. 96, no. 2, pp. 293–320, 2003.

S. Prajna, A. Papachristodoulou, and P. Parrilo, 2002, SOSTOOLS: Sum of squares optimization toolbox for MATLAB, http://www.cds.caltech.edu/sostools

J. Löfberg, “YALMIP: a toolbox for modeling and optimization in MATLAB,” in Proceedings of the 2004 IEEE International Symposium on Computer Aided Control System Design, pp. 284–289, September 2004.

J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,” Optimization Methods and Software, vol. 11, no. 1–4, pp. 625–653, 1999.

M. Kocvara and M. Stingl, PENBMI User's Guide (version 2. 0), http://www.penopt.com, 2005

E. Kaltofen, B. Li, Z. Yang, and L. Zhi, “Exact certification of global optimality of approximate factorizations via rationalizing sums-of-squares with floating point scalars,” in Proceedings of the 21th International Symposium on Symbolic and Algebraic Computation (ISSAC '08), pp. 155–164, ACM, Hagenberg, Austria, July 2008.

E. L. Kaltofen, B. Li, Z. Yang, and L. Zhi, “Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients,” Journal of Symbolic Computation, vol. 47, no. 1, pp. 1–15, 2012.

A. Platzer and E. M. Clarke, “The image computation problem in hybrid systems model checking,” in Proceedings of the 10th International Conference on Hybrid Systems: Computation and Control (HSCC '07), vol. 4416, pp. 473–486, Springer, April 2007.

H. Khalil, Nonlinear Systems, Prentice hall, Upper Saddle River, NJ, USA, 3rd edition, 2002.

J. Bochnak, M. Coste, and M. Roy, Real Algebraic Geometry, Springer, Berlin, Germany, 1998.

G. Chesi, “LMI techniques for optimization over polynomials in control: a survey,” Transactions on Automatic Control, vol. 55, no. 11, pp. 2500–2510, 2010.

F. Leibfritz and E. M. E. Mostafa, “An interior point constrained trust region method for a special class of nonlinear semidefinite programming problems,” SIAM Journal on Optimization, vol. 12, no. 4, pp. 1048–1074, 2002.

M. Kočvara and M. Stingl, “Pennon: a code for convex nonlinear and semidefinite programming,” Optimization Methods & Software, vol. 18, no. 3, pp. 317–333, 2003.

M. Wu and Z. Yang, “Generating invariants of hybrid systems via sums-of-squares of polynomials with rational coefficients,” in Proceedings of the International Workshop on Symbolic-Numeric Computation, pp. 104–111, ACM Press, 2011.

W. Lin, M. Wu, Z. Yang, and Z. Zeng, “Exact safety verification of hybrid systems using sums-of-squares representation,” Science China Information Sciences, 15 pages, 2012.

Z. Yang, M. Wu, and W. Lin, “Exact verification of hybrid systems based on bilinear SOS representation,” 19 pages, 2012, http://arxiv.org/abs/1201.4219

A. Papachristodoulou and S. Prajna, “A tutorial on sum of squares techniques for systems analysis,” in Proceedings of the American Control Conference (ACC '05), pp. 2686–2700, IEEE, June 2005.

J. Löfberg, “Modeling and solving uncertain optimization problems in YALMIP,” in Proceedings of the 17th World Congress, International Federation of Automatic Control (IFAC '08), pp. 1337–1341, July 2008.

U. Topcu, A. Packard, and P. Seiler, “Local stability analysis using simulations and sum-of-squares programming,” Automatica, vol. 44, no. 10, pp. 2669–2675, 2008.