- Email: [email protected]
Domain-of-attraction estimation for uncertain non-polynomial systems
Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Domain-of-attraction estimation for uncertain non-polynomial systems Min Wu a, Zhengfeng Yang a, Wang Lin b,a,⇑ a b
Shanghai Key Laboratory of Trustworthy Computing, East China Normal University, Shanghai 200062, China College of Mathematics and Information Science, Wenzhou University, Zhejiang 325035, China
a r t i c l e
i n f o
Article history: Received 22 May 2013 Received in revised form 3 December 2013 Accepted 4 December 2013 Available online 21 December 2013 Keywords: Non-polynomial systems Domain-of-attraction Polynomial approximation
a b s t r a c t In this paper, we consider the problem of computing estimates of the domain-of-attraction for non-polynomial systems. A polynomial approximation technique, based on multivariate polynomial interpolation and error analysis for remaining functions, is applied to compute an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. The efficiency of the presented method is illustrated by some numerical examples. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Stability for nonlinear control systems plays an important role in control system analysis and design. It will be very useful to know the domain of attraction (DOA) of an equilibrium point, however, this region is usually difficult to find and represent explicitly. Therefore, looking for estimates of the DOA with simple shapes has been a fundamental issue in control system analysis since a long time. Among all the methods, those based on Lyapunov functions are dominant in literature [3,4,6– 8,11,17,22,25,23,10,24,16]. These methods not only yield a Lyapunov function as a stability certificate, but also the corresponding sublevel sets as estimates of the DOA. For polynomial systems, many well-established techniques [6,8,7,11,14,17,22,25,23,10,24,16] are available for computing estimates of DOAs. In [22], a method based on SOS decomposition was presented to find provable DOAs and attractive invariant sets for nonlinear polynomial systems. For odd polynomial systems, [6] employed an LMI-based method to compute the optimal quadratic Lyapunov function for maximizing the volume of the largest estimate of the DOA. To obtain estimates of DOAs of uncertain polynomial systems, the authors of [8] used discretization (in time) to flow invariant sets backwards along the flow of the vector field. In [17], quantifier elimination (QE) method via QEPCAD was also applied to find Lyapunov functions for estimating the DOA. However, these methods cannot be applied directly in practice since most real systems are non-polynomial systems, i.e, their vector fields contain non-polynomial terms. For this kind of systems, only a few approaches have been proposed to deal with the DOA analysis. In [3–5], the author proposed an LMI technique through Taylor expansions as substitution for non-polynomial terms, and this technique can be generalized to compute estimates of DOAs for uncertain non-polynomial systems. In [26], an interval arithmetic approach was proposed. [21,20] suggested a new method, based on quadratic Lyapunov function and the theorem of Ehlich and Zeller. Instead of approximating the original nonlinear model, [19,18] used a maximal Lyapunov function to find an estimation for the DOA of nonlinear hybrid systems where the dynamics is continuous on the boundary of the different regimes of the state space. ⇑ Corresponding author at: College of Mathematics and Information Science, Wenzhou University, Zhejiang 325035, China. Tel.: +86 577 86689216. E-mail address: [email protected] (W. Lin). http://dx.doi.org/10.1016/j.cnsns.2013.12.001 1007-5704/Ó 2013 Elsevier B.V. All rights reserved.
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3045
In this paper, we will consider the problem of stability region analysis of uncertain non-polynomial systems. Through multivariate polynomial interpolation together with the interpolation error analysis, we substitute a non-polynomial system as an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. By computing estimates of the DOA for the resulted uncertain polynomial system, we obtain estimates of the DOA for the original non-polynomial system. Our method is also applicable to the problem of searching for the largest possible guaranteed subset of the DOA via a fixed Lyapunov function. Compared with the classical approximation by Taylor expansions, the error bound obtained using our suggested method is much sharper, which helps to yield a larger estimate of the DOA for a given non-polynomial system. The rest of the paper is organized as follows. In Section 2, some notions related to DOAs are presented. In Section 3, a polynomial approximation method, based on multivariate polynomial interpolation and interpolation error analysis, is proposed to substitute the non-polynomial functions as uncertain polynomials. In Section 4, bilinear SOS programming is applied to estimate DOAs of non-polynomial systems. In Section 5, experiments on some benchmarks are shown to illustrate our suggested method. Section 6 concludes the paper. 2. Problem formulation Consider an autonomous system
x_ ¼ fðxÞ;
ð2:1Þ n
n
where f : D # R ! R is a continuous function defined on an open set D and f satisfies the Lipschitz condition:
kfðxÞ fðyÞk 6 Lkx yk for all x; y 2 D: Denote by /ðt; x0 Þ the solution of (2.1) with the given initial value xð0Þ ¼ x0 . A vector x 2 Rn is an equilibrium point of the system (2.1) if fðxÞ ¼ 0. Since any equilibrium point can be shifted to the origin 0 via a change of variables, we may assume without loss of generality that the equilibrium point of interest occurs at the origin. The equilibrium point 0 of (2.1) is said to be stable, if for any > 0 there exists d such that whenever kx0 k < d we have k/ðt; x0 Þk < for all t > 0; the point 0 is said to be unstable if it is not stable; 0 is asymptotically stable, if, in addition to being stable, there exists d such that limt!1 /ðt; x0 Þ ¼ 0 whenever kx0 k < d; the equilibrium point 0 is globally asymptotically stable, if, in addition to being stable, we have limt!1 /ðt; x0 Þ ¼ 0 for all x0 2 Rn . Globally asymptotic stability is very desirable but is usually difficult to achieve. When the equilibrium point 0 is asymptotically stable, we are interested in determining how far the trajectory of (2.1) can be from 0 and still converge to 0 as t approaches 1. This gives rise to the following definition. Definition 1 (Domain of Attraction). The domain of attraction (DOA) of the equilibrium point 0 for the system (2.1) is defined to be the set fx 2 Rn jlimt!1 /ðt; xÞ ¼ 0g. Usually, no algebraic description for DOAs is available. So researchers are mainly concerned with computing estimates of the DOAs. Many well-established techniques [6,8,7,11,17,22,25,23,10,24,16] are available for computing estimates of DOAs for polynomial (control) systems, i.e., autonomous systems with polynomial vector fields. However, in practice, many autonomous systems often contain non-polynomial terms in their vector fields. Below is an example. Example 1 [12, Example 1.2.1]. Consider the simple pendulum shown in Fig. 1. The motion of the pendulum is described by the following equation
_ ml€h ¼ mg sin h klh; where h denotes the angle subtended by the rod and the vertical axis through the pivot point, l the length of the rod, m the mass of the bob, g the acceleration due to gravity, and k the coefficient of friction. Let us take the state variables as x1 ¼ h, and _ Then the above equation is converted into a non-polynomial system x2 ¼ h.
l θ
mg
Fig. 1. Pendulum.
3046
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
(
x_ 1 ¼ x2 ; x_ 2 ¼ gl sin x1 mk x2 :
For the case of non-polynomial (control) systems, the problem of computing DOAs is still open, and only a few approaches have been proposed to deal with stability region analysis: in [3–5], the authors suggested a way to approximate non-polynomial vector fields by Taylor series expansion at the origin; in [26], an interval arithmetic approach for the estimation of the DOA was proposed; and recently, a method based on the theorem by Ehlich and Zeller was presented in [21,20]. In this paper, we will apply polynomial approximation to transform a non-polynomial system into an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. Therefore, guaranteed subsets of the DOA of the latter system yield those for the original non-polynomial system.
3. Polynomial approximation A key problem in estimating the DOA of a non-polynomial system is how to approximate the involved non-polynomial terms using polynomials, yielding an uncertain polynomial system with the equilibrium 0 being kept. This problem is further reduced to the following problem. Problem 1. Let /ðxÞ : W ! R be a non-polynomial function where W Rn is a bounded subset containing the origin 0. Given d 2 ZP0 , we will find a polynomial pðxÞ with degree d such that the error function r d ðxÞ ¼ /ðxÞ pðxÞ satisfies r d ð0Þ ¼ 0 and the value maxx2W jr d ðxÞj is minimized. The classic method of polynomial approximation is Taylor expansions. Suppose /ðxÞ is a d times continuously differentiable in W. The Taylor expansion of /ðxÞ at the origin 0 is
/ðxÞ ¼
Z 1 Xd1 Da /ð0Þ X d a xa þ x ð1 tÞd1 Da /ðxtÞdt ; j a j¼d jaj¼1 a ! a ! 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} pðxÞ
ð3:1Þ
rd ðxÞ
where a ¼ ða1 ; a2 ; ; an Þ, jaj ¼ a1 þ a2 þ þ an , a! ¼ a1 !a2 ! an !, and Da ¼
@ a1 @ a2 a a @x11 @x22
an
@ @x an . In the above expression, pðxÞ is n
an approximating polynomial of /ðxÞ and the remainder term rd ðxÞ is the error function of this approximation. Clearly, if the size of the region W is small enough, the above Taylor expansion yields a tight bound of rd ðxÞ for all x 2 W. However, when the size of W is large, the associated error bound may be too loose. To obtain a tighter bound, we will apply multivariate polynomial interpolation ([9]) to compute an approximate polynomial pðxÞ of /ðxÞ with a given degree d. Fix the graded lexicographic order in R½x. For the function /ðxÞ, one may find the minimal monomial xc with c :¼ ðc1 ; . . . ; cn Þ 2 ZnP0 , such P –0. Set wðxÞ ¼ /ðxÞ/ð0Þ . Let d 2 ZP0 be such that d > jcj :¼ ni¼1 ci . We construct a mesh M on W with mesh that limx!0 /ðxÞ/ð0Þ xc xc n þ d jcj . Like in [1], the meshes in our paper are spacing s 2 Rþ and mesh points set v ¼ fv 1 ; v 2 ; . . . ; v k g where k ¼ n ~ðxÞ as an approximation either rectangular or simplicial. Then, we apply Lagrange interpolation to construct a polynomial p ~ðv i Þ ¼ wðv i Þ; fori ¼ 1; . . . ; k. Next, we will compute a tight bound of the interof wðxÞ through the interpolation points v, i.e., p ~ðxÞ. Our idea is based on the following lemma. polation error function ~rðxÞ :¼ wðxÞ p Lemma 1 [29, Theorem 3]. Let K Rn be a convex polyhedron, and v 1 ; v 2 ; . . . ; v k and s be the vertices and diameter of K respectively. Suppose that u : K ! R is a continuous and differential function on K, and k ¼ supx2K k 5 uðxÞk. Then for all a1 ; a2 ; . . . ; ak 2 Rþ such that a1 þ a2 þ þ ak ¼ 1, we have
juðxÞ ða1 uðv 1 Þ þ a2 uðv 2 Þ þ þ ak uðv k ÞÞj 6
n ks: nþ1
The following corollary gives an estimated bound of ~r ðxÞ for x 2 W \ M. ~ðxÞ is the Corollary 1. Let s and v :¼ fv 1 ; v 2 ; . . . ; v k g be the mesh spacing and mesh points set of M, respectively. Suppose that p ~ðxÞ is the corresponding error function. Let k ¼ supx2W\M k 5 ~r ðxÞk. interpolation polynomial of wðxÞ through v, and ~r ðxÞ ¼ wðxÞ p Then
j~r ðxÞj 6
n ks for all x 2 M: nþ1
Proof. Clearly, ~rðxÞ is a continuous and differential function on M, and
~rðv 1 Þ ¼ ~r ðv 2 Þ ¼ ¼ ~rðv k Þ ¼ 0:
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3047
Thus, according to Lemma 1, for all a1 ; a2 ; . . . ; av 2 Rþ such that a1 þ a2 þ þ av ¼ 1,
j~rðxÞ ða1~r ðv 1 Þ þ a2~rðv 2 Þ þ þ ak~r ðv k ÞÞj ¼ j~r ðxÞj 6
n ks: nþ1
Therefore, a non-polynomial function /ðxÞ can be relaxed to an uncertain polynomial, as shown in the following theorem. Theorem 1. For a non-polynomial function /ðxÞ with x 2 W, let xc with c 2 ZnP0 be the minimal monomial such that limx!0 /ðxÞ/ð0Þ –0. Let M be a mesh on W with mesh spacing s 2 Rþ and the mesh point set v ¼ fv 1 ; v 2 ; . . . ; v k g, in which xc n þ d jc j ~ðxÞ is the interpolation polynomial of wðxÞ :¼ /ðxÞ/ð0Þ at v with degree 6 d jcj; ~r ðxÞ is the . Suppose that p k¼ xc n associated interpolation error function, and k ¼ supx2W\M k 5 ~rðxÞk. Then for each x 2 W \ M we have /ðxÞ ¼ pðxÞ þ rd ðxÞ, where
~ðxÞ xc pðxÞ ¼ /ð0Þ þ p
and r d ðxÞ ¼ u xc
with juj 6
n ks: nþ1
ð3:2Þ
Clearly, the bound of the error r d ðxÞ in (3.2) depends on the mesh spacing s, which can yield a tighter bound. Furthermore, the bound of r d ðxÞ in (3.2) will converge to zero if d ! 1. Example 2. Consider the function /ðxÞ ¼ cos x with x 2 W ¼ ½1:2; 1:2. We want to compute a polynomial approximation for cos x. Based on Theorem 1, we can obtain an uncertain polynomial with degree 6, where
pðxÞ ¼1 0:5x2 þ 0:0416525x4 0:00134386x6 ; r 6 ðxÞ ¼ux2 ;
0:0000336 6 u 6 0:0000336:
4. Computation of domain of attraction In this section, we will consider an uncertain non-polynomial system of the form:
x_ ¼ fðx; hÞ for all h 2 H Rt ;
ð4:1Þ
where h denotes a vector of uncertainty. Assume that the equilibrium point of interest occurs at the origin 0, i.e, fð0; hÞ ¼ 0 for all h 2 H. Denote by /ðt; x0 ; hÞ the solution of (4.1) for the initial value xð0Þ ¼ x0 and the uncertainty h. The Domain of Attraction (DOA) of the system (4.1) is defined as
n o x 2 Rn jlim /ðt; x; hÞ ¼ 0 for all h 2 H : t!1
Lemma 1 in [22] can be modified a bit to compute guaranteed subsets of the DOA for (4.1) through Lyapunov functions, as described in the following theorem. Theorem 2 [25, Proposition 2.1]. If there exists a continuously differentiable function V : Rn ! R such that
8 XV :¼ fx 2 Rn : VðxÞ 6 1g is bounded; > > > > < Vð0Þ ¼ 0; > VðxÞ > 0; 8x 2 XV n f0g; > > > : VðxÞ _ fðx; hÞ < 0; 8x 2 XV n f0g; 8h 2 H ¼ @V @x
then XV is an invariant subset of the DOA. When the equilibrium 0 is asymptotically stable, the set Xc is clearly a subset of the DOA since every trajectory starting in Xc remains in Xc and approaches 0 as t ! 1. And, if the equilibrium 0 is globally asymptotically stable then the DOA will be the whole space Rn . To enlarge the estimate XV given in Theorem 2, [22] defined a variable sized region
Pb ¼ fx 2 Rn : gðxÞ 6 bg P with gðxÞ a fixed and positive definite polynomial in R½x, for instance, gðxÞ ¼ ni¼1 x2i , and maximize b subject to the constraint P b # XV and the constraints in Theorem 2. Thus, the problem of computing XV can be transformed into the following problem:
max b s:t: XV
is bounded;
9 > > > > > > > > =
Vð0Þ ¼ 0; > VðxÞ > 0; 8x 2 XV n f0g; > > > > _ > VðxÞ ¼ @V fðx; hÞ < 0; 8 x 2 X n f0g; 8 h 2 H > V @x > ; gðxÞ 6 b VðxÞ 6 1:
ð4:2Þ
3048
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
Suppose that the non-polynomial system (4.1) has the following form
x_ i ¼ fi ðx; hÞ ¼ fi0 ðx; hÞ þ
k X
f ij ðx; hÞ/ij ðxÞ;
i ¼ 1; . . . ; n;
ð4:3Þ
j¼1
where fij : Rn Rt ! R are polynomials in x for j ¼ 0; . . . ; k, and /ij : Rn ! R are non-polynomial functions for j ¼ 1; . . . ; k. Using the polynomial approximation technique in Section 3, we can replace each non-polynomial term /ij ðxÞ by an uncertain polynomial pij ðxÞ þ uij xcij with the bound juij j 6 bij . This gives rise to the following uncertain polynomial system:
x_ i ¼ ^f i ðx; hÞ ¼ fi0 ðx; hÞ þ
k X
f ij ðx; hÞðpij ðxÞ þ uij xcij Þ
ð4:4Þ
j¼1
where juij j 6 bij for i ¼ 1; . . . ; n. It is not hard to prove that the set of all trajectories of the system (4.3) is a subset of that of the system (4.4), and, consequently, the DOA of (4.4) is actually a subset of that of the system (4.3). Furthermore, the tighter the bound bij is, the closer the DOA of (4.4) is to the DOA of the original system. Next, we consider how to find an optimal estimate of the DOA for the uncertain system (4.3) through computing the XV in _ ¼ @V (4.2). Remark that the constraint VðxÞ fðx; hÞ < 0 in (4.3) may involve non-polynomial terms due to the existence of @x _ /ij ðxÞ’s. In this situation, replacing the above constraint by VðxÞ ^fðx; hÞ < 0, the problem (4.2) can be relaxed as follow. ¼ @V @x T e ðxÞ is a solution of the following polynomial Theorem 3. Let ^fðx; hÞ ¼ ð^f 1 ðx; hÞ; . . . ; ^f n ðx; hÞÞ with ^f i ðx; hÞ given in (4.4). If V optimization
9 > > > > > > > =
max b V;b
s:t: Vð0Þ ¼ 0; VðxÞ > 0 8x 2 XV n f0g; ^fðx; hÞ < 0 8x 2 X n f0g;
@V @x
V
8h 2 H;
gðxÞ 6 b VðxÞ 6 1;
ð4:5Þ
> > > > 8uij 2 fbij g; > > > ;
e ðxÞ 6 1g is an invariant subset of the DOA for (4.4), and therefore an invariant subset of DOA for (4.3). then Xe :¼ fx 2 Rn : V V Proof. By construction of ^f i ðx; hÞ’s, we have
8x 2 XV n f0g;
8h 2 H ;
_ ~ij 2 ½bij ; bij : VðxÞ 9u ¼
@V ^ fðx; hÞ: @x
Clearly, if the constraints in (4.5) are fulfilled, the conditions in (4.2) also hold. Therefore, Xe is certainly a subset of XV in V (4.2). h Assume that H is a semialgebraic set. For simplicity, we suppose H ¼ fh 2 H : wðhÞ P 0g. By rewriting the third, fourth and fifth constraints into equivalent empty set conditions, the condition (4.5) is transformed as
8 Vð0Þ ¼ 0; > > > < fx 2 Rn : VðxÞ 6 1; x–0; VðxÞ 6 0g ¼ ;; > fx 2 Rn : VðxÞ 6 1; x–0; wðhÞ P 0; @V ^fðxÞ P 0g ¼ ;; > @x > : n fx 2 R : gðxÞ 6 b; VðxÞ P 1; VðxÞ–1g ¼ ;:
ð4:6Þ
8uij 2 fbij g;
As stated in [22], Stengle’s Positivstellensatz [2] be applied directly to solve (4.6). However, from the computational point of view, it is more efficient to replace all the inequations in (4.6) by inequalities of the form f P 0 or f 6 0. This can be done by introducing constants d 2 Rþ and polynomials of the form [22] lðxÞ ¼ Rni¼1 i xm i , where i 2 Rþ and m is assumed to be even. For example, by using d1 and l1 ðxÞ, the second condition in (4.6) can be relaxed as
fx 2 Rn : VðxÞ 1 6 0; VðxÞ l1 ðxÞ þ d1 6 0g ¼ ;: Therefore, the problem (4.6) can be transformed into the following feasibility problem:
9 > > > > > > > =
max b V;b
s:t: Vð0Þ ¼ 0; fx 2 Rn : VðxÞ 1 6 0; fx 2 R : VðxÞ 1 6 0;
VðxÞ l1 ðxÞ þ d1 6 0g ¼ ;; @V ^fðx; hÞ l2 ðxÞ þ d2 6 0;
fx 2 Rn : gðxÞ b 6 0;
VðxÞ 1 d3 P 0g ¼ ;:
n
@x
wðhÞ P 0g ¼ ;;
> > > > ;
> 8uij 2 fbij g > >
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3049
Suppose that R½x is the set of sum of squares (SOS) polynomials in R½x. Since the constraints in the above problem involve no equations and inequations, only a special case of Stengle’s Positivstellensatz is needed, as shown in the following corollary. Corollary 2. Let F ¼ ffi gi¼1;...;r be a set of polynomials in R½x. The semi-algebraic set
fx 2 Rn : fi ðxÞ P 0; i ¼ 1; . . . ; rg is empty if and only if there exist polynomials s0 ; s1 ; . . . ; sl 2 R½x such that
s0 þ
l X sj b j ¼ 0 j¼1
where bj 2 f1t1 frtr : ti 2 ZP0 . Applying Corollary 2, and removing all the crossing products of the involved inequalities, we obtain the following relaxed problem:
max b V;b
s:t: Vð0Þ ¼ 0;
r0 ðxÞ þ r1 ðxÞð1 VðxÞÞ þ r2 ðxÞðVðxÞ þ l1 ðxÞ d1 Þ ¼ 0;
9 > > > > > > > > =
> > > f;^hðxÞ þ l2 ðxÞ d2 ¼ 0; > k0 ðxÞ þ k1 ðxÞð1 VðxÞÞ þ k2 ðxÞ @V > @x > > > ; q0 ðxÞ þ q1 ðxÞðb gðxÞÞ þ q2 ðxÞðVðxÞ 1 d3 Þ ¼ 0;
ð4:7Þ
with ri ðxÞ; ki ðxÞ; qi ðxÞ 2 R½x; i ¼ 0; . . . ; 2 and for any uij 2 fbij g. P Suppose that the Lyapunov function VðxÞ to be computed is a polynomial of degree d and has the form VðxÞ ¼ a ca xa P where ca 2 R, xa ¼ xa11 xann and a ¼ ða1 ; . . . ; an Þ 2 ZnP0 with ni¼1 ai 6 d. The decision variables of the problem (4.7) are b and the coefficients of all the unknown polynomials occurred in (4.7), such as VðxÞ; ri ðxÞ; ki ðxÞ and qi ðxÞ. Clearly, some nonlinear terms which are products of the undetermined coefficients will occur in (4.7), which yields a non-convex bilinear matrix inequalities (BMI) problem. To solve BMI problems, either a Matlab package PENBMI solver [13], which combines the (exterior) penalty and (interior) barrier method with the augmented Lagrangian method, can be applied directly, or an iterative method can be applied by fixing b and the polynomials alternatively, which leads to a sequential convex LMI problem. The reader can refer to [28] for more details. Remark that, the above proposed method is also applicable to computing the largest possible estimate of the DOA for a non-polynomial system (2.1) at 0 through a fixed Lyapunov function VðxÞ. Let Xc :¼ fx 2 Rn : VðxÞ 6 cg. We will compute
Xc ¼ fx 2 Rn : VðxÞ 6 c g; _ where c ¼ supfc 2 R : VðxÞ < 0; 8x 2 Xc n f0gg. Due to the existence of non-polynomial terms, c cannot be computed explicitly. Instead, we will compute the lower and upper bounds of c as follows. Replacing the involved non-polynomial /ij by uncertain polynomials, computing the lower bound cd of c can be relaxed as the following problem:
cd ¼ max c c;x
s:t:
@V @x
^fðxÞ < 0;
ð4:8Þ
8x 2 Xc n f0g;
8uij 2 fbij g;
where d is the degree of the interpolation polynomials. Clearly, cd will converge to c when d tends to 1. Next, we will search for a tight upper bound td of c . To achieve this, let us look for td such that, for each uij 2 fbij g, the constant semi-algebraic system
VðxÞ td 6 0 ^ x–0 ^
@V ^ fðxÞ P 0 @x
has real solutions, which implies that Xtd is not an estimate of the DOA. Based on bisection, packages RegularChains, DISCOVERER [27] and RAGLib [15]. 5. Experiments Let us present some examples of DOA analysis of non-polynomial systems. Example 3 [4, Example 1]. Consider a non-polynomial system
x_ 1 ¼ x1 þ x2 þ 0:5ðex1 1Þ; x_ 2 ¼ x1 x2 þ x1 x2 þ x1 cos x1 :
ð4:9Þ
td can be computed by Maple
3050
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
To estimate the DOA of this system, we need to approximate the occurred non-polynomial terms ex1 and cos x1 by uncertain polynomials. Based on the technique in Section 3, we obtain
8 x e 1 ¼ 1 þ ð1:0000004 þ u1 Þx1 þ þ 0:0014482244x61 ; > > > > 2 4 6 > > < cos x1 ¼ 1 ð0:5 þ u2 Þx1 þ 0:041669352x1 0:0013878601x1 ; 0:6 6 x1 6 0:6; > > > > 3 106 6 u1 6 3 106 ; > > : 1:2 106 6 u2 6 1:2 106 ; and the associated uncertain polynomial system. We first consider a fixed Lyapunov function Vðx1 ; x2 Þ ¼ x21 þ x22 . For the given degree 6 of interpolation polynomials, after solving the corresponding SOS programming (4.8), we obtain the lower bound c6 ¼ 0:321064 of c , which is an improvement over the results in [4] with the lower bound 0.3210. Furthermore, by solving the problem (4.9) we obtain a tight upper bound t6 ¼ 0:3216 of c . Next, we estimate the DOA with variable Lyapunov functions. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . For deg V ¼ 2, solving the SOS programming (4.7) with BMI constraints yields
Vðx1 ; x2 Þ ¼ 0:56678683x21 þ 0:23598133x1 x2 þ 0:92086339x22 ; b ¼ 1:0453916; which is an improvement over the result from [4] where b ¼ 1:0404. Similarly, for deg V ¼ 4 solving the SOS programming (4.7) with BMI constraints yields b ¼ 1:4001306 and
Vðx1 ; x2 Þ ¼ 0:068693712x21 þ þ 0:27723966x41 þ 0:2167918x42 ; which is an improvement over the result from [4] where b ¼ 1:2769. Therefore, XV is an estimate of the DOA of the given system. Fig. 2 shows the results obtained with Lyapunov functions of degrees 2 and 4. h Example 4 [4, Example 2]. Consider a non-polynomial system
x_ 1 ¼ x2 ; x_ 2 ¼ 0:2x2 þ 0:81 sin x1 cos x1 sin x1 :
Using the technique in Section 3, we obtain the approximations of the non-polynomial terms sin x1 and cos x1 as follows
8 sin 2x1 ¼ ð2 þ u1 Þx1 1:3333091x31 þ 0:26625372x51 0:023889715x71 ; > > > > 3 5 7 > > < sin x1 ¼ ð1 þ u2 Þx1 0:16666643x1 þ 0:008331760x1 0:00019471928x1 ; 0:84 6 x1 6 0:84; > > > > 5:4 105 6 u1 6 5:4 105 ; > > : 1 107 6 u2 6 1 107 ;
and the associated uncertain polynomial system. We first fix the Lyapunov function Vðx1 ; x2 Þ ¼ x21 þ x1 x2 þ 4x22 . Let the degree of interpolation polynomials be 7. Solving the corresponding SOS programming (4.8), we obtain the results for the lower bound c7 ¼ 0:69922 of c , which is an improvement over the result from [4] where the lower bound was 0.6990. Furthermore, by solving the problem (4.9) we obtain a tighter upper bound t7 ¼ 0:6998 of c .
(a)
(b) Fig. 2. Results of Example 4.
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3051
(b)
(a)
Fig. 3. Results of Example 5.
We then estimate the DOA with variable Lyapunov functions. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . When deg V ¼ 2, by solving the SOS programming (4.7) with BMI constraints, we obtain
Vðx1 ; x2 Þ ¼ 1:01636667x21 þ 0:84993333x1 x2 þ 3:40233333x22 ; b ¼ 0:287706; which is an improvement over the result from [4] where b ¼ 0:2809. Similarly, when deg V ¼ 4, by solving the SOS programming (4.7) with BMI constraints, we obtain b ¼ 1:92156 and
Vðx1 ; x2 Þ ¼ 0:053849691x21 þ þ 0:11144243x41 þ 0:058855229x42 ; which is an improvement over the result from [4] where b ¼ 1:1236. Therefore, Xe is an estimate of the DOA of the given V system. Fig. 3 shows the results obtained with Lyapunov functions of degrees 2 and 4. h Example 5 [12, Example 2.2]. Consider an uncertain non-polynomial system
x_ 1 ¼ x2 ; x_ 2 ¼ hx2 10 sin x1 ;
for 0:2 6 h 6 1. Based on the technique in Section 3, we obtain an approximation of the non-polynomial term sin x1 as follows
8 3 5 7 > < sin x1 ¼ ð1 þ u1 Þx1 0:1666426901x1 þ 0:008282118073x1 0:0001751721223x1 ; 2:4 6 x1 6 2:4; > : 4:365606 104 6 u1 6 4:365606 104 ; and the associated uncertain polynomial system. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . For deg V ¼ 4, solving the SOS programming (4.7) with BMI constraints yields b ¼ 0:66552836 and
Vðx1 ; x2 Þ ¼ 1:1629845x21 þ þ 0:51014802x41 þ 0:010528x42 : Then XV is an estimate of the DOA of the given system.
h
6. Conclusion In this paper, we present a method on stability region analysis of non-polynomial systems via Lyapunov functions. A polynomial approximation technique, based on multivariate polynomial interpolation and error analysis, is applied to compute an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. To estimate DOA of the uncertain polynomial system, we apply Positivstellensatz to transform polynomial optimization problem into the corresponding (bilinear) sum of squares programming, which can be solved using the PENBMI solver or iterative method. The efficiency of the presented method is illustrated by some numerical examples. Acknowledgments We would like to thank anonymous reviewers for the very valuable comments. This material is supported in part by the National Natural Science Foundation of China under Grants 11371143, 91118007, and 61021004, the Innovation Program of
3052
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
Shanghai Municipal Education Commission under Grant 14ZZ046, the National High Technology Research and Development Program of China under Grant 2011AA010101, and the Natural Science Foundation of Zhejiang Province under Grant Q13F020041. References [1] Asarin E, Dang T, Girard A. Reachability analysis of nonlinear systems using conservative approximation. In: Proceedings of the 6th international conference on hybrid systems: computation and control. Springer-Verlag; 2003. p. 20–35. [2] Bochnak J, Coste M, Roy M. Real algebraic geometry. Springer-Verlag; 1998. [3] Chesi G. Domain of attraction: estimates for non-polynomial systems via LMIs. In: Proc. 16th IFAC world congress on automatic, control. IFAC; 2005. [4] Chesi G. Estimating the domain of attraction for non-polynomial systems via LMI optimizations. Automatica 2009;45(6):1536–41. [5] Chesi G. Domain of attraction: analysis and control via SOS programming. Springer; 2011. [6] Chesi G, Garulli A, Tesi A, Vicino A. LMI-based computation of optimal quadratic Lyapunov functions for odd polynomial systems. Int J Robust Nonlinear Control 2005;15(1):35–49. [7] Chiang H, Thorp J. Stability regions of nonlinear dynamical systems: a constructive methodology. IEEE Trans Autom Control 1989;34(12):1229–41. [8] Cruck E, Moitie R, Seube N. Estimation of basins of attraction for uncertain systems with affine and Lipschitz dynamics. Dyn Control 2001;11(3):211–27. [9] Gasca M, Sauer T. On the history of multivariate polynomial interpolation. J Comput Appl Math 2000;122(1):23–35. [10] Hachicho O, Tibken B. Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments. Proceedings of the 41st IEEE conference on decision and control, vol. 3. IEEE; 2002. p. 3150–5. [11] Jarvis-Wloszek Z. Lyapunov Based Analysis and Controller Synthesis for Polynomial Systems Using Sum-Of-Squares Optimization. PhD thesis, University of California, 2003. [12] Khalil H. Nonlinear systems. 3rd ed. New Jewsey: Prentice hall; 2002. [13] Kocˇvara, M., and Stingl, M. PENBMI User’s Guide (Version 2.0). http://www.penopt.com, 2005 [14] Liu Y, Pang G. The basin of attraction of the liu system. Commun Nonlinear Sci Numer Simul 2011;16(4):2065–71. [15] Mohab, S.E.D. Raglib (Real Algebraic Library Maple package)., 2003 [16] Prajna S, Parrilo P, Rantzer A. Nonlinear control synthesis by convex optimization. IEEE Trans Autom Control 2004;49(2):310–4. [17] Prakash S, Vanualailai J, Soma T. Obtaining approximate region of asymptotic stability by computer algebra: a case study. South Pac J Nat Appl Sci 2002;20(1):56–61. [18] Rozgonyi S, Hangos K. Domain of attraction analysis of a controlled hybrid reactor model. Ann Nucl Eng 2011;38(5):969–75. [19] Rozgonyi S, Hangos KM, Szederkényi G. Determining the domain of attraction of hybrid nonlinear systems using maximal lyapunov functions. Kybernetika 2010;46(1):19–37. [20] Saleme A, Tibken B. A new method to estimate a guaranteed subset of the domain of attraction for non-polynomial systems. American control conference. IEEE; 2012. p. 2577–2582. [21] Saleme A, Tibken B, Warthenpfuhl S, Selbach C. Estimation of the domain of attraction for non-polynomial systems: a novel method. In: Proceedings of the 18th IFAC world congress. Milano, Italy: IFAC; 2011. p. 10976–10981. [22] Tan W, Packard A. Stability region analysis using sum of squares programming. American control conference. IEEE; 2006. p. 2297–2302. [23] Tibken B. Estimation of the domain of attraction for polynomial systems via LMIs. Proceedings of the 39th IEEE conference on decision and control, vol. 4. IEEE; 2000. p. 3860–4. [24] Tibken B, Dilaver K. Computation of subsets of the domain of attraction for polynomial systems. Proceedings of the 41st IEEE conference on decision and control, vol. 3. IEEE; 2002. p. 2651–6. [25] Topcu U, Packard AK, Seiler P, Balas GJ. Robust region-of-attraction estimation. IEEE Trans Autom Control 2010;55(1):137–42. [26] Warthenpfuhl S, Tibken B, Mayer S. An interval arithmetic approach for the estimation of the domain of attraction. In: International symposium on computer-aided control system design. IEEE; 2010. p. 1999–2004. [27] Xia B. DISCOVERER: a tool for solving semi-algebraic systems. ACM Commun Comput Algebra 2007;41(3):102–3. [28] Yang Z, Wu M, Lin W. Exact safety verification of hybrid systems based on bilinear SOS representation. http://arxiv.org/abs/1201.4219, 2012. [29] Zeng Z, Zhang J. A mechanical proof to a geometric inequality of Zirakzadeh through rectangular partition of polyhedra (in Chinese). J Syst Sci Math Sci 2010;30(11):1430–58.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Domain-of-attraction estimation for uncertain non-polynomial systems Min Wu a, Zhengfeng Yang a, Wang Lin b,a,⇑ a b
Shanghai Key Laboratory of Trustworthy Computing, East China Normal University, Shanghai 200062, China College of Mathematics and Information Science, Wenzhou University, Zhejiang 325035, China
a r t i c l e
i n f o
Article history: Received 22 May 2013 Received in revised form 3 December 2013 Accepted 4 December 2013 Available online 21 December 2013 Keywords: Non-polynomial systems Domain-of-attraction Polynomial approximation
a b s t r a c t In this paper, we consider the problem of computing estimates of the domain-of-attraction for non-polynomial systems. A polynomial approximation technique, based on multivariate polynomial interpolation and error analysis for remaining functions, is applied to compute an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. The efficiency of the presented method is illustrated by some numerical examples. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Stability for nonlinear control systems plays an important role in control system analysis and design. It will be very useful to know the domain of attraction (DOA) of an equilibrium point, however, this region is usually difficult to find and represent explicitly. Therefore, looking for estimates of the DOA with simple shapes has been a fundamental issue in control system analysis since a long time. Among all the methods, those based on Lyapunov functions are dominant in literature [3,4,6– 8,11,17,22,25,23,10,24,16]. These methods not only yield a Lyapunov function as a stability certificate, but also the corresponding sublevel sets as estimates of the DOA. For polynomial systems, many well-established techniques [6,8,7,11,14,17,22,25,23,10,24,16] are available for computing estimates of DOAs. In [22], a method based on SOS decomposition was presented to find provable DOAs and attractive invariant sets for nonlinear polynomial systems. For odd polynomial systems, [6] employed an LMI-based method to compute the optimal quadratic Lyapunov function for maximizing the volume of the largest estimate of the DOA. To obtain estimates of DOAs of uncertain polynomial systems, the authors of [8] used discretization (in time) to flow invariant sets backwards along the flow of the vector field. In [17], quantifier elimination (QE) method via QEPCAD was also applied to find Lyapunov functions for estimating the DOA. However, these methods cannot be applied directly in practice since most real systems are non-polynomial systems, i.e, their vector fields contain non-polynomial terms. For this kind of systems, only a few approaches have been proposed to deal with the DOA analysis. In [3–5], the author proposed an LMI technique through Taylor expansions as substitution for non-polynomial terms, and this technique can be generalized to compute estimates of DOAs for uncertain non-polynomial systems. In [26], an interval arithmetic approach was proposed. [21,20] suggested a new method, based on quadratic Lyapunov function and the theorem of Ehlich and Zeller. Instead of approximating the original nonlinear model, [19,18] used a maximal Lyapunov function to find an estimation for the DOA of nonlinear hybrid systems where the dynamics is continuous on the boundary of the different regimes of the state space. ⇑ Corresponding author at: College of Mathematics and Information Science, Wenzhou University, Zhejiang 325035, China. Tel.: +86 577 86689216. E-mail address: [email protected] (W. Lin). http://dx.doi.org/10.1016/j.cnsns.2013.12.001 1007-5704/Ó 2013 Elsevier B.V. All rights reserved.
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3045
In this paper, we will consider the problem of stability region analysis of uncertain non-polynomial systems. Through multivariate polynomial interpolation together with the interpolation error analysis, we substitute a non-polynomial system as an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. By computing estimates of the DOA for the resulted uncertain polynomial system, we obtain estimates of the DOA for the original non-polynomial system. Our method is also applicable to the problem of searching for the largest possible guaranteed subset of the DOA via a fixed Lyapunov function. Compared with the classical approximation by Taylor expansions, the error bound obtained using our suggested method is much sharper, which helps to yield a larger estimate of the DOA for a given non-polynomial system. The rest of the paper is organized as follows. In Section 2, some notions related to DOAs are presented. In Section 3, a polynomial approximation method, based on multivariate polynomial interpolation and interpolation error analysis, is proposed to substitute the non-polynomial functions as uncertain polynomials. In Section 4, bilinear SOS programming is applied to estimate DOAs of non-polynomial systems. In Section 5, experiments on some benchmarks are shown to illustrate our suggested method. Section 6 concludes the paper. 2. Problem formulation Consider an autonomous system
x_ ¼ fðxÞ;
ð2:1Þ n
n
where f : D # R ! R is a continuous function defined on an open set D and f satisfies the Lipschitz condition:
kfðxÞ fðyÞk 6 Lkx yk for all x; y 2 D: Denote by /ðt; x0 Þ the solution of (2.1) with the given initial value xð0Þ ¼ x0 . A vector x 2 Rn is an equilibrium point of the system (2.1) if fðxÞ ¼ 0. Since any equilibrium point can be shifted to the origin 0 via a change of variables, we may assume without loss of generality that the equilibrium point of interest occurs at the origin. The equilibrium point 0 of (2.1) is said to be stable, if for any > 0 there exists d such that whenever kx0 k < d we have k/ðt; x0 Þk < for all t > 0; the point 0 is said to be unstable if it is not stable; 0 is asymptotically stable, if, in addition to being stable, there exists d such that limt!1 /ðt; x0 Þ ¼ 0 whenever kx0 k < d; the equilibrium point 0 is globally asymptotically stable, if, in addition to being stable, we have limt!1 /ðt; x0 Þ ¼ 0 for all x0 2 Rn . Globally asymptotic stability is very desirable but is usually difficult to achieve. When the equilibrium point 0 is asymptotically stable, we are interested in determining how far the trajectory of (2.1) can be from 0 and still converge to 0 as t approaches 1. This gives rise to the following definition. Definition 1 (Domain of Attraction). The domain of attraction (DOA) of the equilibrium point 0 for the system (2.1) is defined to be the set fx 2 Rn jlimt!1 /ðt; xÞ ¼ 0g. Usually, no algebraic description for DOAs is available. So researchers are mainly concerned with computing estimates of the DOAs. Many well-established techniques [6,8,7,11,17,22,25,23,10,24,16] are available for computing estimates of DOAs for polynomial (control) systems, i.e., autonomous systems with polynomial vector fields. However, in practice, many autonomous systems often contain non-polynomial terms in their vector fields. Below is an example. Example 1 [12, Example 1.2.1]. Consider the simple pendulum shown in Fig. 1. The motion of the pendulum is described by the following equation
_ ml€h ¼ mg sin h klh; where h denotes the angle subtended by the rod and the vertical axis through the pivot point, l the length of the rod, m the mass of the bob, g the acceleration due to gravity, and k the coefficient of friction. Let us take the state variables as x1 ¼ h, and _ Then the above equation is converted into a non-polynomial system x2 ¼ h.
l θ
mg
Fig. 1. Pendulum.
3046
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
(
x_ 1 ¼ x2 ; x_ 2 ¼ gl sin x1 mk x2 :
For the case of non-polynomial (control) systems, the problem of computing DOAs is still open, and only a few approaches have been proposed to deal with stability region analysis: in [3–5], the authors suggested a way to approximate non-polynomial vector fields by Taylor series expansion at the origin; in [26], an interval arithmetic approach for the estimation of the DOA was proposed; and recently, a method based on the theorem by Ehlich and Zeller was presented in [21,20]. In this paper, we will apply polynomial approximation to transform a non-polynomial system into an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. Therefore, guaranteed subsets of the DOA of the latter system yield those for the original non-polynomial system.
3. Polynomial approximation A key problem in estimating the DOA of a non-polynomial system is how to approximate the involved non-polynomial terms using polynomials, yielding an uncertain polynomial system with the equilibrium 0 being kept. This problem is further reduced to the following problem. Problem 1. Let /ðxÞ : W ! R be a non-polynomial function where W Rn is a bounded subset containing the origin 0. Given d 2 ZP0 , we will find a polynomial pðxÞ with degree d such that the error function r d ðxÞ ¼ /ðxÞ pðxÞ satisfies r d ð0Þ ¼ 0 and the value maxx2W jr d ðxÞj is minimized. The classic method of polynomial approximation is Taylor expansions. Suppose /ðxÞ is a d times continuously differentiable in W. The Taylor expansion of /ðxÞ at the origin 0 is
/ðxÞ ¼
Z 1 Xd1 Da /ð0Þ X d a xa þ x ð1 tÞd1 Da /ðxtÞdt ; j a j¼d jaj¼1 a ! a ! 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} pðxÞ
ð3:1Þ
rd ðxÞ
where a ¼ ða1 ; a2 ; ; an Þ, jaj ¼ a1 þ a2 þ þ an , a! ¼ a1 !a2 ! an !, and Da ¼
@ a1 @ a2 a a @x11 @x22
an
@ @x an . In the above expression, pðxÞ is n
an approximating polynomial of /ðxÞ and the remainder term rd ðxÞ is the error function of this approximation. Clearly, if the size of the region W is small enough, the above Taylor expansion yields a tight bound of rd ðxÞ for all x 2 W. However, when the size of W is large, the associated error bound may be too loose. To obtain a tighter bound, we will apply multivariate polynomial interpolation ([9]) to compute an approximate polynomial pðxÞ of /ðxÞ with a given degree d. Fix the graded lexicographic order in R½x. For the function /ðxÞ, one may find the minimal monomial xc with c :¼ ðc1 ; . . . ; cn Þ 2 ZnP0 , such P –0. Set wðxÞ ¼ /ðxÞ/ð0Þ . Let d 2 ZP0 be such that d > jcj :¼ ni¼1 ci . We construct a mesh M on W with mesh that limx!0 /ðxÞ/ð0Þ xc xc n þ d jcj . Like in [1], the meshes in our paper are spacing s 2 Rþ and mesh points set v ¼ fv 1 ; v 2 ; . . . ; v k g where k ¼ n ~ðxÞ as an approximation either rectangular or simplicial. Then, we apply Lagrange interpolation to construct a polynomial p ~ðv i Þ ¼ wðv i Þ; fori ¼ 1; . . . ; k. Next, we will compute a tight bound of the interof wðxÞ through the interpolation points v, i.e., p ~ðxÞ. Our idea is based on the following lemma. polation error function ~rðxÞ :¼ wðxÞ p Lemma 1 [29, Theorem 3]. Let K Rn be a convex polyhedron, and v 1 ; v 2 ; . . . ; v k and s be the vertices and diameter of K respectively. Suppose that u : K ! R is a continuous and differential function on K, and k ¼ supx2K k 5 uðxÞk. Then for all a1 ; a2 ; . . . ; ak 2 Rþ such that a1 þ a2 þ þ ak ¼ 1, we have
juðxÞ ða1 uðv 1 Þ þ a2 uðv 2 Þ þ þ ak uðv k ÞÞj 6
n ks: nþ1
The following corollary gives an estimated bound of ~r ðxÞ for x 2 W \ M. ~ðxÞ is the Corollary 1. Let s and v :¼ fv 1 ; v 2 ; . . . ; v k g be the mesh spacing and mesh points set of M, respectively. Suppose that p ~ðxÞ is the corresponding error function. Let k ¼ supx2W\M k 5 ~r ðxÞk. interpolation polynomial of wðxÞ through v, and ~r ðxÞ ¼ wðxÞ p Then
j~r ðxÞj 6
n ks for all x 2 M: nþ1
Proof. Clearly, ~rðxÞ is a continuous and differential function on M, and
~rðv 1 Þ ¼ ~r ðv 2 Þ ¼ ¼ ~rðv k Þ ¼ 0:
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3047
Thus, according to Lemma 1, for all a1 ; a2 ; . . . ; av 2 Rþ such that a1 þ a2 þ þ av ¼ 1,
j~rðxÞ ða1~r ðv 1 Þ þ a2~rðv 2 Þ þ þ ak~r ðv k ÞÞj ¼ j~r ðxÞj 6
n ks: nþ1
Therefore, a non-polynomial function /ðxÞ can be relaxed to an uncertain polynomial, as shown in the following theorem. Theorem 1. For a non-polynomial function /ðxÞ with x 2 W, let xc with c 2 ZnP0 be the minimal monomial such that limx!0 /ðxÞ/ð0Þ –0. Let M be a mesh on W with mesh spacing s 2 Rþ and the mesh point set v ¼ fv 1 ; v 2 ; . . . ; v k g, in which xc n þ d jc j ~ðxÞ is the interpolation polynomial of wðxÞ :¼ /ðxÞ/ð0Þ at v with degree 6 d jcj; ~r ðxÞ is the . Suppose that p k¼ xc n associated interpolation error function, and k ¼ supx2W\M k 5 ~rðxÞk. Then for each x 2 W \ M we have /ðxÞ ¼ pðxÞ þ rd ðxÞ, where
~ðxÞ xc pðxÞ ¼ /ð0Þ þ p
and r d ðxÞ ¼ u xc
with juj 6
n ks: nþ1
ð3:2Þ
Clearly, the bound of the error r d ðxÞ in (3.2) depends on the mesh spacing s, which can yield a tighter bound. Furthermore, the bound of r d ðxÞ in (3.2) will converge to zero if d ! 1. Example 2. Consider the function /ðxÞ ¼ cos x with x 2 W ¼ ½1:2; 1:2. We want to compute a polynomial approximation for cos x. Based on Theorem 1, we can obtain an uncertain polynomial with degree 6, where
pðxÞ ¼1 0:5x2 þ 0:0416525x4 0:00134386x6 ; r 6 ðxÞ ¼ux2 ;
0:0000336 6 u 6 0:0000336:
4. Computation of domain of attraction In this section, we will consider an uncertain non-polynomial system of the form:
x_ ¼ fðx; hÞ for all h 2 H Rt ;
ð4:1Þ
where h denotes a vector of uncertainty. Assume that the equilibrium point of interest occurs at the origin 0, i.e, fð0; hÞ ¼ 0 for all h 2 H. Denote by /ðt; x0 ; hÞ the solution of (4.1) for the initial value xð0Þ ¼ x0 and the uncertainty h. The Domain of Attraction (DOA) of the system (4.1) is defined as
n o x 2 Rn jlim /ðt; x; hÞ ¼ 0 for all h 2 H : t!1
Lemma 1 in [22] can be modified a bit to compute guaranteed subsets of the DOA for (4.1) through Lyapunov functions, as described in the following theorem. Theorem 2 [25, Proposition 2.1]. If there exists a continuously differentiable function V : Rn ! R such that
8 XV :¼ fx 2 Rn : VðxÞ 6 1g is bounded; > > > > < Vð0Þ ¼ 0; > VðxÞ > 0; 8x 2 XV n f0g; > > > : VðxÞ _ fðx; hÞ < 0; 8x 2 XV n f0g; 8h 2 H ¼ @V @x
then XV is an invariant subset of the DOA. When the equilibrium 0 is asymptotically stable, the set Xc is clearly a subset of the DOA since every trajectory starting in Xc remains in Xc and approaches 0 as t ! 1. And, if the equilibrium 0 is globally asymptotically stable then the DOA will be the whole space Rn . To enlarge the estimate XV given in Theorem 2, [22] defined a variable sized region
Pb ¼ fx 2 Rn : gðxÞ 6 bg P with gðxÞ a fixed and positive definite polynomial in R½x, for instance, gðxÞ ¼ ni¼1 x2i , and maximize b subject to the constraint P b # XV and the constraints in Theorem 2. Thus, the problem of computing XV can be transformed into the following problem:
max b s:t: XV
is bounded;
9 > > > > > > > > =
Vð0Þ ¼ 0; > VðxÞ > 0; 8x 2 XV n f0g; > > > > _ > VðxÞ ¼ @V fðx; hÞ < 0; 8 x 2 X n f0g; 8 h 2 H > V @x > ; gðxÞ 6 b VðxÞ 6 1:
ð4:2Þ
3048
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
Suppose that the non-polynomial system (4.1) has the following form
x_ i ¼ fi ðx; hÞ ¼ fi0 ðx; hÞ þ
k X
f ij ðx; hÞ/ij ðxÞ;
i ¼ 1; . . . ; n;
ð4:3Þ
j¼1
where fij : Rn Rt ! R are polynomials in x for j ¼ 0; . . . ; k, and /ij : Rn ! R are non-polynomial functions for j ¼ 1; . . . ; k. Using the polynomial approximation technique in Section 3, we can replace each non-polynomial term /ij ðxÞ by an uncertain polynomial pij ðxÞ þ uij xcij with the bound juij j 6 bij . This gives rise to the following uncertain polynomial system:
x_ i ¼ ^f i ðx; hÞ ¼ fi0 ðx; hÞ þ
k X
f ij ðx; hÞðpij ðxÞ þ uij xcij Þ
ð4:4Þ
j¼1
where juij j 6 bij for i ¼ 1; . . . ; n. It is not hard to prove that the set of all trajectories of the system (4.3) is a subset of that of the system (4.4), and, consequently, the DOA of (4.4) is actually a subset of that of the system (4.3). Furthermore, the tighter the bound bij is, the closer the DOA of (4.4) is to the DOA of the original system. Next, we consider how to find an optimal estimate of the DOA for the uncertain system (4.3) through computing the XV in _ ¼ @V (4.2). Remark that the constraint VðxÞ fðx; hÞ < 0 in (4.3) may involve non-polynomial terms due to the existence of @x _ /ij ðxÞ’s. In this situation, replacing the above constraint by VðxÞ ^fðx; hÞ < 0, the problem (4.2) can be relaxed as follow. ¼ @V @x T e ðxÞ is a solution of the following polynomial Theorem 3. Let ^fðx; hÞ ¼ ð^f 1 ðx; hÞ; . . . ; ^f n ðx; hÞÞ with ^f i ðx; hÞ given in (4.4). If V optimization
9 > > > > > > > =
max b V;b
s:t: Vð0Þ ¼ 0; VðxÞ > 0 8x 2 XV n f0g; ^fðx; hÞ < 0 8x 2 X n f0g;
@V @x
V
8h 2 H;
gðxÞ 6 b VðxÞ 6 1;
ð4:5Þ
> > > > 8uij 2 fbij g; > > > ;
e ðxÞ 6 1g is an invariant subset of the DOA for (4.4), and therefore an invariant subset of DOA for (4.3). then Xe :¼ fx 2 Rn : V V Proof. By construction of ^f i ðx; hÞ’s, we have
8x 2 XV n f0g;
8h 2 H ;
_ ~ij 2 ½bij ; bij : VðxÞ 9u ¼
@V ^ fðx; hÞ: @x
Clearly, if the constraints in (4.5) are fulfilled, the conditions in (4.2) also hold. Therefore, Xe is certainly a subset of XV in V (4.2). h Assume that H is a semialgebraic set. For simplicity, we suppose H ¼ fh 2 H : wðhÞ P 0g. By rewriting the third, fourth and fifth constraints into equivalent empty set conditions, the condition (4.5) is transformed as
8 Vð0Þ ¼ 0; > > > < fx 2 Rn : VðxÞ 6 1; x–0; VðxÞ 6 0g ¼ ;; > fx 2 Rn : VðxÞ 6 1; x–0; wðhÞ P 0; @V ^fðxÞ P 0g ¼ ;; > @x > : n fx 2 R : gðxÞ 6 b; VðxÞ P 1; VðxÞ–1g ¼ ;:
ð4:6Þ
8uij 2 fbij g;
As stated in [22], Stengle’s Positivstellensatz [2] be applied directly to solve (4.6). However, from the computational point of view, it is more efficient to replace all the inequations in (4.6) by inequalities of the form f P 0 or f 6 0. This can be done by introducing constants d 2 Rþ and polynomials of the form [22] lðxÞ ¼ Rni¼1 i xm i , where i 2 Rþ and m is assumed to be even. For example, by using d1 and l1 ðxÞ, the second condition in (4.6) can be relaxed as
fx 2 Rn : VðxÞ 1 6 0; VðxÞ l1 ðxÞ þ d1 6 0g ¼ ;: Therefore, the problem (4.6) can be transformed into the following feasibility problem:
9 > > > > > > > =
max b V;b
s:t: Vð0Þ ¼ 0; fx 2 Rn : VðxÞ 1 6 0; fx 2 R : VðxÞ 1 6 0;
VðxÞ l1 ðxÞ þ d1 6 0g ¼ ;; @V ^fðx; hÞ l2 ðxÞ þ d2 6 0;
fx 2 Rn : gðxÞ b 6 0;
VðxÞ 1 d3 P 0g ¼ ;:
n
@x
wðhÞ P 0g ¼ ;;
> > > > ;
> 8uij 2 fbij g > >
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3049
Suppose that R½x is the set of sum of squares (SOS) polynomials in R½x. Since the constraints in the above problem involve no equations and inequations, only a special case of Stengle’s Positivstellensatz is needed, as shown in the following corollary. Corollary 2. Let F ¼ ffi gi¼1;...;r be a set of polynomials in R½x. The semi-algebraic set
fx 2 Rn : fi ðxÞ P 0; i ¼ 1; . . . ; rg is empty if and only if there exist polynomials s0 ; s1 ; . . . ; sl 2 R½x such that
s0 þ
l X sj b j ¼ 0 j¼1
where bj 2 f1t1 frtr : ti 2 ZP0 . Applying Corollary 2, and removing all the crossing products of the involved inequalities, we obtain the following relaxed problem:
max b V;b
s:t: Vð0Þ ¼ 0;
r0 ðxÞ þ r1 ðxÞð1 VðxÞÞ þ r2 ðxÞðVðxÞ þ l1 ðxÞ d1 Þ ¼ 0;
9 > > > > > > > > =
> > > f;^hðxÞ þ l2 ðxÞ d2 ¼ 0; > k0 ðxÞ þ k1 ðxÞð1 VðxÞÞ þ k2 ðxÞ @V > @x > > > ; q0 ðxÞ þ q1 ðxÞðb gðxÞÞ þ q2 ðxÞðVðxÞ 1 d3 Þ ¼ 0;
ð4:7Þ
with ri ðxÞ; ki ðxÞ; qi ðxÞ 2 R½x; i ¼ 0; . . . ; 2 and for any uij 2 fbij g. P Suppose that the Lyapunov function VðxÞ to be computed is a polynomial of degree d and has the form VðxÞ ¼ a ca xa P where ca 2 R, xa ¼ xa11 xann and a ¼ ða1 ; . . . ; an Þ 2 ZnP0 with ni¼1 ai 6 d. The decision variables of the problem (4.7) are b and the coefficients of all the unknown polynomials occurred in (4.7), such as VðxÞ; ri ðxÞ; ki ðxÞ and qi ðxÞ. Clearly, some nonlinear terms which are products of the undetermined coefficients will occur in (4.7), which yields a non-convex bilinear matrix inequalities (BMI) problem. To solve BMI problems, either a Matlab package PENBMI solver [13], which combines the (exterior) penalty and (interior) barrier method with the augmented Lagrangian method, can be applied directly, or an iterative method can be applied by fixing b and the polynomials alternatively, which leads to a sequential convex LMI problem. The reader can refer to [28] for more details. Remark that, the above proposed method is also applicable to computing the largest possible estimate of the DOA for a non-polynomial system (2.1) at 0 through a fixed Lyapunov function VðxÞ. Let Xc :¼ fx 2 Rn : VðxÞ 6 cg. We will compute
Xc ¼ fx 2 Rn : VðxÞ 6 c g; _ where c ¼ supfc 2 R : VðxÞ < 0; 8x 2 Xc n f0gg. Due to the existence of non-polynomial terms, c cannot be computed explicitly. Instead, we will compute the lower and upper bounds of c as follows. Replacing the involved non-polynomial /ij by uncertain polynomials, computing the lower bound cd of c can be relaxed as the following problem:
cd ¼ max c c;x
s:t:
@V @x
^fðxÞ < 0;
ð4:8Þ
8x 2 Xc n f0g;
8uij 2 fbij g;
where d is the degree of the interpolation polynomials. Clearly, cd will converge to c when d tends to 1. Next, we will search for a tight upper bound td of c . To achieve this, let us look for td such that, for each uij 2 fbij g, the constant semi-algebraic system
VðxÞ td 6 0 ^ x–0 ^
@V ^ fðxÞ P 0 @x
has real solutions, which implies that Xtd is not an estimate of the DOA. Based on bisection, packages RegularChains, DISCOVERER [27] and RAGLib [15]. 5. Experiments Let us present some examples of DOA analysis of non-polynomial systems. Example 3 [4, Example 1]. Consider a non-polynomial system
x_ 1 ¼ x1 þ x2 þ 0:5ðex1 1Þ; x_ 2 ¼ x1 x2 þ x1 x2 þ x1 cos x1 :
ð4:9Þ
td can be computed by Maple
3050
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
To estimate the DOA of this system, we need to approximate the occurred non-polynomial terms ex1 and cos x1 by uncertain polynomials. Based on the technique in Section 3, we obtain
8 x e 1 ¼ 1 þ ð1:0000004 þ u1 Þx1 þ þ 0:0014482244x61 ; > > > > 2 4 6 > > < cos x1 ¼ 1 ð0:5 þ u2 Þx1 þ 0:041669352x1 0:0013878601x1 ; 0:6 6 x1 6 0:6; > > > > 3 106 6 u1 6 3 106 ; > > : 1:2 106 6 u2 6 1:2 106 ; and the associated uncertain polynomial system. We first consider a fixed Lyapunov function Vðx1 ; x2 Þ ¼ x21 þ x22 . For the given degree 6 of interpolation polynomials, after solving the corresponding SOS programming (4.8), we obtain the lower bound c6 ¼ 0:321064 of c , which is an improvement over the results in [4] with the lower bound 0.3210. Furthermore, by solving the problem (4.9) we obtain a tight upper bound t6 ¼ 0:3216 of c . Next, we estimate the DOA with variable Lyapunov functions. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . For deg V ¼ 2, solving the SOS programming (4.7) with BMI constraints yields
Vðx1 ; x2 Þ ¼ 0:56678683x21 þ 0:23598133x1 x2 þ 0:92086339x22 ; b ¼ 1:0453916; which is an improvement over the result from [4] where b ¼ 1:0404. Similarly, for deg V ¼ 4 solving the SOS programming (4.7) with BMI constraints yields b ¼ 1:4001306 and
Vðx1 ; x2 Þ ¼ 0:068693712x21 þ þ 0:27723966x41 þ 0:2167918x42 ; which is an improvement over the result from [4] where b ¼ 1:2769. Therefore, XV is an estimate of the DOA of the given system. Fig. 2 shows the results obtained with Lyapunov functions of degrees 2 and 4. h Example 4 [4, Example 2]. Consider a non-polynomial system
x_ 1 ¼ x2 ; x_ 2 ¼ 0:2x2 þ 0:81 sin x1 cos x1 sin x1 :
Using the technique in Section 3, we obtain the approximations of the non-polynomial terms sin x1 and cos x1 as follows
8 sin 2x1 ¼ ð2 þ u1 Þx1 1:3333091x31 þ 0:26625372x51 0:023889715x71 ; > > > > 3 5 7 > > < sin x1 ¼ ð1 þ u2 Þx1 0:16666643x1 þ 0:008331760x1 0:00019471928x1 ; 0:84 6 x1 6 0:84; > > > > 5:4 105 6 u1 6 5:4 105 ; > > : 1 107 6 u2 6 1 107 ;
and the associated uncertain polynomial system. We first fix the Lyapunov function Vðx1 ; x2 Þ ¼ x21 þ x1 x2 þ 4x22 . Let the degree of interpolation polynomials be 7. Solving the corresponding SOS programming (4.8), we obtain the results for the lower bound c7 ¼ 0:69922 of c , which is an improvement over the result from [4] where the lower bound was 0.6990. Furthermore, by solving the problem (4.9) we obtain a tighter upper bound t7 ¼ 0:6998 of c .
(a)
(b) Fig. 2. Results of Example 4.
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
3051
(b)
(a)
Fig. 3. Results of Example 5.
We then estimate the DOA with variable Lyapunov functions. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . When deg V ¼ 2, by solving the SOS programming (4.7) with BMI constraints, we obtain
Vðx1 ; x2 Þ ¼ 1:01636667x21 þ 0:84993333x1 x2 þ 3:40233333x22 ; b ¼ 0:287706; which is an improvement over the result from [4] where b ¼ 0:2809. Similarly, when deg V ¼ 4, by solving the SOS programming (4.7) with BMI constraints, we obtain b ¼ 1:92156 and
Vðx1 ; x2 Þ ¼ 0:053849691x21 þ þ 0:11144243x41 þ 0:058855229x42 ; which is an improvement over the result from [4] where b ¼ 1:1236. Therefore, Xe is an estimate of the DOA of the given V system. Fig. 3 shows the results obtained with Lyapunov functions of degrees 2 and 4. h Example 5 [12, Example 2.2]. Consider an uncertain non-polynomial system
x_ 1 ¼ x2 ; x_ 2 ¼ hx2 10 sin x1 ;
for 0:2 6 h 6 1. Based on the technique in Section 3, we obtain an approximation of the non-polynomial term sin x1 as follows
8 3 5 7 > < sin x1 ¼ ð1 þ u1 Þx1 0:1666426901x1 þ 0:008282118073x1 0:0001751721223x1 ; 2:4 6 x1 6 2:4; > : 4:365606 104 6 u1 6 4:365606 104 ; and the associated uncertain polynomial system. Suppose gðx1 ; x2 Þ ¼ x21 þ x22 . For deg V ¼ 4, solving the SOS programming (4.7) with BMI constraints yields b ¼ 0:66552836 and
Vðx1 ; x2 Þ ¼ 1:1629845x21 þ þ 0:51014802x41 þ 0:010528x42 : Then XV is an estimate of the DOA of the given system.
h
6. Conclusion In this paper, we present a method on stability region analysis of non-polynomial systems via Lyapunov functions. A polynomial approximation technique, based on multivariate polynomial interpolation and error analysis, is applied to compute an uncertain polynomial system, whose set of trajectories contains that of the original non-polynomial system. To estimate DOA of the uncertain polynomial system, we apply Positivstellensatz to transform polynomial optimization problem into the corresponding (bilinear) sum of squares programming, which can be solved using the PENBMI solver or iterative method. The efficiency of the presented method is illustrated by some numerical examples. Acknowledgments We would like to thank anonymous reviewers for the very valuable comments. This material is supported in part by the National Natural Science Foundation of China under Grants 11371143, 91118007, and 61021004, the Innovation Program of
3052
M. Wu et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 3044–3052
Shanghai Municipal Education Commission under Grant 14ZZ046, the National High Technology Research and Development Program of China under Grant 2011AA010101, and the Natural Science Foundation of Zhejiang Province under Grant Q13F020041. References [1] Asarin E, Dang T, Girard A. Reachability analysis of nonlinear systems using conservative approximation. In: Proceedings of the 6th international conference on hybrid systems: computation and control. Springer-Verlag; 2003. p. 20–35. [2] Bochnak J, Coste M, Roy M. Real algebraic geometry. Springer-Verlag; 1998. [3] Chesi G. Domain of attraction: estimates for non-polynomial systems via LMIs. In: Proc. 16th IFAC world congress on automatic, control. IFAC; 2005. [4] Chesi G. Estimating the domain of attraction for non-polynomial systems via LMI optimizations. Automatica 2009;45(6):1536–41. [5] Chesi G. Domain of attraction: analysis and control via SOS programming. Springer; 2011. [6] Chesi G, Garulli A, Tesi A, Vicino A. LMI-based computation of optimal quadratic Lyapunov functions for odd polynomial systems. Int J Robust Nonlinear Control 2005;15(1):35–49. [7] Chiang H, Thorp J. Stability regions of nonlinear dynamical systems: a constructive methodology. IEEE Trans Autom Control 1989;34(12):1229–41. [8] Cruck E, Moitie R, Seube N. Estimation of basins of attraction for uncertain systems with affine and Lipschitz dynamics. Dyn Control 2001;11(3):211–27. [9] Gasca M, Sauer T. On the history of multivariate polynomial interpolation. J Comput Appl Math 2000;122(1):23–35. [10] Hachicho O, Tibken B. Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments. Proceedings of the 41st IEEE conference on decision and control, vol. 3. IEEE; 2002. p. 3150–5. [11] Jarvis-Wloszek Z. Lyapunov Based Analysis and Controller Synthesis for Polynomial Systems Using Sum-Of-Squares Optimization. PhD thesis, University of California, 2003. [12] Khalil H. Nonlinear systems. 3rd ed. New Jewsey: Prentice hall; 2002. [13] Kocˇvara, M., and Stingl, M. PENBMI User’s Guide (Version 2.0). http://www.penopt.com, 2005 [14] Liu Y, Pang G. The basin of attraction of the liu system. Commun Nonlinear Sci Numer Simul 2011;16(4):2065–71. [15] Mohab, S.E.D. Raglib (Real Algebraic Library Maple package).