# Padé-like approximation and its application in domain of attraction estimation

Padé-like approximation and its application in domain of attraction estimation Abstract In this paper, a rational approximation of multivariable functions is proposed. Since the proposed approach, hereafter called Padé-like approximation, uses the coefficients of Taylor series expansion at a point, it can be considered as a novel extension of conventional Padé approximation for multivariable functions. This extension is free from some of previous extensions imperfections. Similar to the conventional Padé approximation for one-variable functions, the outcome of the proposed method can reveal the points where the multi-variable function is non-analytic. The proposed method, whose accuracy is analytically guaranteed, could find a wide range of applications wherever an accurate approximation is needed. In the line of estimating domain of attraction in non-linear systems, Padé-like approximation of an optimal Lyapunov function (OLF) is obtained from the corresponding partial differential equation. Using the approximate OLF, an estimation of domain of attraction is obtained which is less conservative compared to those based on the previous techniques. 1. Introduction For one-variable functions, the Padé approximation is a well-known method and is considered as the best rational function approximation (see Baker, 1975). Padé approximation could be valid even in the region where the Taylor series does not converge (see Guillaume & Huard, 2000). Due to the rationality of approximation, it can be used for approximating the non-analytic points of a function. Also, construction by Taylor series coefficients is another great property of Padé approximation. That is because, in some cases, Taylor series coefficients of an even non-explicitly known function may be obtained from the equation (possibly a partial differential equation, PDE) satisfying it. It should be noted that approximate solution of an equation is usually different from the expansion of the exact solution. The aforementioned property of the Padé approximation, however, causes no difference between the equation and the function approximations. It is important to note that for multi-variable functions, the matter is not so simple, which is emergence of the fact that the extension of Padé approximation to multi-variable functions is not straightforward. There are some extensions of Padé approximation for multi-variable functions such as those posed in Jones (1976), Cuyt (1983), Guillaume (1997, 1998), Guillaume et al. (1998), Kida (1990), Guillaume & Huard (2000). With regard to the estimation of domain of attraction, none of the mentioned methods are sufficient. Methods like the one in Jones (1976) need some symmetry property in coefficients of Taylor series and do not work for some functions, especially those vanish at the origin. On the other hand, methods like the one posed in Cuyt (1983) give a high-order polynomial, even for a low order approximation, and in the meantime, the origin is always one of the roots of approximation denominator. The recent fact makes this method useless for those purposes which need the study of non-analyticity property of functions. This matter becomes important for our application, i.e. the computation of $$D_A$$ (domain of attraction), which may be related to non-analytic points of a multi-variable function (see Balint et al., 2009). Guillaume (1997, 1998) has proposed a nested Padé approximation technique which suffers from non-homogeneity of approximation degree for each variable. In the meantime, the denominator is not a polynomial and furthermore, it could cause some computational costs. Guillaume et al. (1998) have used the least square technique for obtaining the coefficients of numerator and denominator of the rational approximation. Through this method, Taylor series expansion coefficients of the difference of a function and its approximation become marginal, but not zero. These circumstances convinced the authors to develop a novel rational approximation and use its denominators roots as an estimation of domain of attraction. By computation of $$D_A$$, we mean the solution of the problem of finding initial conditions which converge to an asymptotic stable equilibrium point. A linear model is usually an approximation of a real non-linear system about its operating point. Hence, the validity of all subsequent analyses and designs is local. As a result, after stabilization, one has to determine the corresponding domain of attraction (see Yazdanpanah et al., 1999). From control engineering point of view, the importance of finding the $$D_A$$ seems obvious in analysis and synthesis. From application point of view, computation of $$D_A$$ has found numerous applications in different fields of science and technology. Power plant systems, chemical reactors, aircrafts and robotic systems are some typical examples of $$D_A$$ applications in engineering. Biological systems, economy and ecology are some other applied fields where $$D_A$$ plays a crucial role in the ground of behavioural analysis (see Genesio et al., 1985, Yazdanpanah et al., 1998). Moreover, enlargement of the closed loop $$D_A$$ has been a goal for some controller design strategies (see Bakhtiari & Yazdanpanah, 2005, Benzaouia, 2005, Hashemzadeh & Yazdanpanah, 2006). $$D_A$$ has rich and challenging theoretical aspects as it has a practically important and beneficial insight. It is clear that $$D_A$$ is a complicated set and, except very special cases, it has no analytical and simple representation. Therefore, massive research effort has been done on the $$D_A$$ estimation during past decades. $$D_A$$ estimation is extensively considered in the literature and many different methods have been applied for it. These methods can be divided into two general categories: Lyapunov-based methods and non-Lyapunov-based methods (see Genesio et al., 1985). This classification is very common and is almost used by every paper in the field. Lyapunov-based methods try to find a Lyapunov function based on which an invariant set is proposed as the estimation of $$D_A$$. Trajectory reversing method (see Genesio, 1984) and topological-based method (see Chiang et al., 1988) can be mentioned as some non-Lyapunov-based methods. Genesio et al. (1985) have reviewed the relevant works published till 1985. Papers related to the subject have reviewed briefly these methods in their introduction sections. Zubov theorem and the estimation methods based on that are early approaches in this field (see Genesio et al., 1985, Melchor-Aguilar & Niculescu, 2007 and Balint et al., 2009 for a class of time-delay systems). Vannelli & Vidyasagar (1985) introduced maximal Lyapunov function and proposed an estimation method which is categorized under Lyapunov-based methods. This method is one of the few ones which uses a rational Lyapunov function. Hachicho (2007) has extended this method through an optimization technique with LMI (linear matrix inequality) constraints. An important branch of Lyapunov-based method which is concentrated in the past decade is optimization-based methods. These methods try to formulate the finding of a Lyapunov function which maximizes $$D_A$$ estimation—towards its exact domain—as an optimization problem. This optimization problem, in general, is a non-convex problem (see Tibken, 2000). These methods are applied to odd polynomial systems (see Chesi & Garulli, 2005), polynomial systems (see Tibken, 2000, Chesi, 2004, Valmórbida, 2009) and some non-polynomial systems (see Chesi, 2009). Square programming is also a new and powerful optimization-based method that is used for $$D_A$$ estimation (see Topcu & Packard, 2007, Tan & Packard, 2008, Topcu et al., 2010). Although $$D_A$$ estimation is extensively studied, almost all of the obtained results are conservative. For example, the results regarding unbounded $$D_A$$s are not satisfactory. Consequently, $$D_A$$ estimation is still an open problem for research. Balint et al. have presented one of the interesting theorems in this area which presents the conditions for exact determination of $$D_A$$ (see Balint, 1985, Kaslik et al., 2003, Balint et al., 2009). Simply, Balint et al. (2009) states that exact $$D_A$$ for a system is coincides with the domain of analyticity of a special Lyapunov function. This Lyapunov function which is introduced by a PDE is called the optimal Lyapunov function (OLF). It is worth noting that such an OLF is different from the similarly called function in the Zubov theore. Like many other PDEs, obtaining the closed-form solution of the OLF—but in a very special cases—generally is not possible. Furthermore, for this type of PDEs, solving PDE by characteristic method (Evans, 1998)—a known analytical way—is equivalent to solving the dynamical equations of the system. Hence, finding the approximate solution for this PDE becomes significant. Rational polynomials have the potential of revealing the possible non-analytic nature of the functions. Consequently, it seems that using a rational function as an approximation of the OLF can approximate non-analytical points of OLF based on which an estimation of $$D_A$$ may be obtained. This fact is the linkage key of Padé-like approximation and estimation of $$D_A$$. In the line of aforementioned justification, here, first, a novel rational approximation is developed for multi-variable functions. This approximation, like conventional Padé approximation, uses Taylor series coefficients and preserves the great property of conventional Padé approximation discussed above. This method is quite different from previous ones and resolves some of their deficiencies. In this paper, this method is used for the purpose of estimation of $$D_A$$, but like other approximators, has a wide potential applications. The convergence of this method is shown mathematically in Theorem 3.1. In the next step, this approximation method is applied to approximate an OLF. The roots of denominator of this rational approximation are near to the non-analytical points of OLF. Thus, the denominator of this approximation is used to present an estimation of $$D_A$$. Applying this method to some known sample systems shows its power, especially for unbounded $$D_A$$s. Some applications of $$D_A$$ estimation need a guarantee for the estimated region to be inside the exact domain. One possible way to present this guarantee is to prove that the roots of approximation denominator converge to the non-analytical points of the OLF. This paper does not pursue the mentioned line of proof because of the following reason. Usually, convergence theorems of approximation methods discuss the differences between approximate and exact solutions for a sufficiently high order of approximation, while a high order approximation may be useless due to the computational burden. For such cases in the estimation process, instead of studying the convergence of roots of approximation denominator to the non-analytical points of OLF, we propose a method in Section 5 to put the estimated domain inside the exact domain. However, no analytical proof is given for the convergence of roots of denominator of the approximation to the non-analytical points of the OLF. This paper is organized as follows. Section 2 is dedicated to preliminary materials. Constructive introduction to Padé-like function approximation (PFA) is considered in Section 3. The proposed method for $$D_A$$ estimation is presented in Section 4. Some stability analysis for estimation result is posed in Section 5. Finally, this paper is ended with conclusion in Section 6. 2. Preliminary The following nomenclature is used frequently in the sequel. A multi-index $$k$$ is an element of $$({\mathbb{N} \cup \{0\}})^n$$ or loosely, $$k=(k_1,k_2,\dots,k_n)$$. Some standard multi-index notations with $$x \in \mathbb{R}^n, \quad i,k \in ({\mathbb{N} \cup \{0\}})^n, \quad n\in \mathbb{N}$$ are as follows (see Krantz & Parks, 2002): |k| =k1+k2+⋯+kn,k! =k1!k2!…kn!,i<kifij<kj∀j=1,2,…,n,i≤kifij≤kj∀j=1,2,…,n,(ki) =(k1i1)(k2i2)…(knin),xk =x1k1x2k2…xnkn,∂|k|∂xk =∂k1∂xk1∂k2∂xk2…∂kn∂xkn. (2.1) Also, the following notation can simplify the expressions: ci =c(i), c:(N∪{0})n→R∑i=0k(.) ≡∑i1=0k1∑i2=0k2…∑in=0kn(.)∑|i|=m∈N(.) ≡∑i1=0m∑i2=0m−i1…∑in=m−i1−i2...−in−1(.)∑|i|=0m∈N(.) ≡∑i1=0m∑i2=0m−i1…∑in=0m−∑j=0n−1ij(.). (2.2) It is convenient to set a convention: ci=0,∀i∉(N∪{0})n. (2.3) The expression like $$x^k$$ where $$x\in \mathbb{R}^n,k\in ({\mathbb{N} \cup \{0\}})^n$$ is called a monomial of order $$|k|$$. It is straight to show that there are (n+|k|−1n−1)=(n+|k|−1)!(|k|)!(n−1)! (2.4)$$n$$-variable monomials of order $$|k|$$. Also, there are (n+|k|n)$$n$$-variable monomials of order less than or equal to $$|k|$$. Similarly, if $$V(x)$$ is at least a $$|k|$$-times continuously differentiable function, then it has $${n+|k|-1 \choose n-1}$$ derivatives of order $$|k|$$ such as: ∂|k|V∂xk=∂|k|V∂k1x1,∂k2x2,…,∂knxn. Throughout this paper, $$\|.\|$$ is used for 2-norm, i.e.: ∀x∈Rn,‖x‖=x12+x22+⋯+xn2 and a ball centred at the point $$x_0$$ with radius $$\epsilon \in \mathbb{R}_{\geq 0}$$ is considered as follows: Bϵ(x0)={x∈Rn|‖x−x0‖≤ϵ}. A function $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ is said to be real analytic at point $$\alpha\in\mathbb{R}^n$$ if there exist a neighbour of $$\alpha$$ such as $$U\subset\mathbb{R}^n$$ and a convergent power series such that (see Krantz & Parks, 2002): ∀x∈U;V(x)=∑|k|=0∞ck(x−α)k,ck∈R. This series representation is called Taylor series expansion about $$\alpha$$. Proposition 2.1 Let $$V$$ be a real analytic function at point $$\alpha$$. Then, its Taylor series coefficients satisfy the following (see Krantz & Parks, 2002): ck=1k! ∂|k|V∂xk|x=0. Proposition 2.2 Function $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ is real analytic on an open set $$U\subset \mathbb{R}^n$$ if and only if for each $$\alpha \in U$$ there are an open ball $$B_{\epsilon}(\alpha) \subset U$$ and constants $$C>0,R>0$$ such that (see Krantz & Parks, 2002): ∀x∈Bϵ(α),|∂|k|V∂xk|≤Ck!R|k|. For functions $$f,\phi: \mathfrak{S}\rightarrow \mathbb{R}$$, it is said that $$f(x)=O(\phi(x))$$ if there is a positive constant $$c_o$$ such that (see Bruijn, 1961): |f(x)|≤co|ϕ(x)|∀x∈S. Consequently, for function $$f(x):\mathbb{R}^n \rightarrow \mathbb{R}$$, it is said that $$f(x)=O\left(\|x\|^p\right)$$ as $$\|x\| \rightarrow 0$$ if there exist constants $$M$$ and $$c_o$$ such that: |f(x)|≤co‖x‖p∀‖x‖≤M. (2.5) In this paper, the definition of domain of attraction is considered as follows. Definition 2.1 (Khalil, 2002) Let $$x=0$$ be an asymptotic stable equilibrium point of: x˙=f(x),f:D⊆Rn→Rn, (2.6) where $$f$$ is locally Lipschitz and $$D$$ contains the origin. Let $$\phi(x;t)$$ be the solution of this system at time $$t$$ with initial condition $$x$$. The domain of attraction $$(D_A)$$ is defined as: DA={x∈D|ϕ(x;t)→0 as t→+∞}. Here, the theorem which relates the exact $$D_A$$ to the domain of analyticity of a special Lyapunov function is reviewed (see Balint, 1985; Balint et al., 2009). Theorem 2.1 (Balint et al., 2009) Consider system (2.6). Let $$f$$ be real analytic and $${\partial{f} \over \partial{x}}|_{x=0}$$ be Hurwitz. Then, domain of attraction of $$x=0$$ coincides with the domain of analyticity of the real analytical function $$V$$ defined by: ⟨∇V,f⟩=−‖x‖2V(0)=0. (2.7)$$V$$ is strictly positive on $$D_A-\{0\}$$ and $$V(x)\rightarrow\infty$$ as $$x\rightarrow\partial D_A$$ or $$\|x\|\rightarrow\infty$$. The above unique function V is called OLF. The following lemma which is based on Lyapunov theory is the foundation of several estimation methods. Lemma 2.1 (Vidayasagar, 1993) Assume there exists a smooth function $$V:\mathbb{R}^n \rightarrow \mathbb{R}$$ and a positive constant $$c$$ for system (2.6) such that: V(0)=0,V(x)>0∀x∈Rn−0 and ΩV,c:={x∈Rn|V(x)≤c} is bounded. Now, if the condition: ΩV,c⊂{x∈Rn|∇V(x).f(x)≤0} is satisfied then $$\it{\Omega}_{V,c}$$ is an invariant subset of $$D_A$$. 3. PFA Despite the existence of Padé approximation for one-variable functions, the current extensions for multivariable functions are not effectual in the line of this work’s goal. So, here, an attempt is made to propose a new extension of Padé approximation which is suited for OLF approximation. This extension like the conventional Padé approximation is built from Taylor series coefficients. So, in this section, it is assumed that the expansion of $$V$$ in the form: V(x)=∑|i|=0∞cixi (3.1) is in hand, and then this function is tried to be approximated via an appropriate rational function. Let, W(x)=∑|i|=0naaixi∑|i|=0nbbixi (3.2) be the required rational approximation. The parameter $$b_0$$ is set to $$1$$ for uniqueness of approximation. Now, in some manner, coefficients $$a_i$$s and $$b_i$$s must be determined. According to the previous argument, there are $${n+n_a} \choose {n}$$ and $${{n+n_b} \choose {n}}-1$$ unknown coefficients for numerator and denominator of (3.2), respectively. In this regard, (n+nan)+(n+nbn)−1 (3.3) independent equations are needed. $$a_i$$s and $$b_i$$s must be determined so that: V(x)−W(x)=0. (3.4) However, it should be noted that, it is usually impossible to annihilate the right-hand side of the above equation; so, $$a_i$$s and $$b_i$$s must be determined such that the right-hand side of (3.4) is possibly free of lower order terms. In this connection, consider: W(x)−V(x)=∑|i|=0naaixi∑|i|=0nbbixi−V(x)=0⇒(∑|i|=0naaixi)−(∑|i|=0nbbixi)V(x)=0 or equivalently: (∑|i|=0naaixi)−(∑|i|=0nbbixi)(∑|i|=0∞cixi)=0. (3.5) The above equation may be rearranged as: ∑|i|=0na(ai−di)xi−∑|i|=na+1∞dixi=0, (3.6) where di={∑j1=0i1∑j2=0i2…∑jn=0inbj.ci−jif i≤nb∑|j|=0nbbj.ci−jif i>nb.  With convention (2.3), we can consider $$d_i$$ for all $$i$$ as: di=∑|j|=0nbbj.ci−j. Since $$b_0=1$$, one has: di=ci+∑|j|=1nbbj.ci−j. (3.7) By setting, $$a_i-d_i=0 \quad \text{for} \quad 0\leq |i| \leq n_a$$, the first sum in Equation (3.6) will disappear. Due to the coefficients of $$a_i$$s, these sets of equations are linearly independent and consequently solvable. Actually, if $$b_i$$s are known, then equations $$a_i-d_i=0$$ determine $$a_i$$s. So, further $${{n+n_b} \choose {n}}-1$$ independent equations must be chosen from the second sum in Equation (3.6). Considering Equation (3.7), when $$|i|$$ is fixed, it is better to choose $$d_i=0$$, whose corresponding $$c_i$$ is greater in the absolute value sense. On the other hand, lower powers of $$|i|$$ have preference. So, at first, $$|i|=n_a+1$$ is considered and due to the descending order of $$c_i$$s’ absolute value, the corresponding $$d_i=0$$ is added to the equation set such that the recent equation is independent and consistent from the previous ones. Henceforth, if the number of equations is inadequate, this operation is iterated with the next of $$|i|$$. This operation continues until an enough number of equations are obtained. The method can be summarized in the following pseudo algorithm. PFA algorithm: 1. Consider (3.2) as a rational approximation of (3.1). The parameters $$a_i$$s and $$b_i$$s are unknown variables which must be determined. In the sequel, $$\mathscr{S}$$ represents the set of equations needed for finding the parameters $$a_i$$s and $$b_i$$s. 2. For $$\quad 0\leq|i|\leq n_a$$, add $$a_i-d_i=0$$ to the equation set $$\mathscr{S}.$$ 3. $$k=n_a+1$$. 4. For $$|i|=k$$, arrange $$c_i$$ from largest to the smallest absolute value. Consider $$c_j, j \in \mathbb{N}$$, where $$j=1,2,\dots,{{n-k+1} \choose {k}}$$ shows the mentioned descending order. (a) $$j=1$$. (b) Consider $$c_i$$ and $$d_i$$ corresponding to $$c_j$$. If $$d_i=0$$ is independent and consistent to the equation set $$\mathscr{S}$$; add it to $$\mathscr{S}$$. (c) If number of equations in set $$\mathscr{S}$$ is not equal to (3.3) and $$j<{{n-k+1} \choose {k}}$$, then $$j \leftarrow j+1$$ and go to 4(b). Else, go to (5). 5. If number of equations in set $$\mathscr{S}$$ is not equal to (3.3), then $$k \leftarrow k+1$$ and go to (5). 6. Solve the linear equation set $$\mathscr{S}$$ for unknown variables $$a_i$$s and $$b_i$$s. Remark 3.1 If enough equations cannot be found in steps (5) and (5) in the above algorithm, it can be concluded that the order of approximation is higher than the main function. As an exaggerated example, consider $$V(x)=\frac{1}{1-x}=\sum_{i=0}^{\infty}{x^i}$$ and $$W(x)=\frac{a_1x+a_0}{b_2x^2+b_1x+1}$$. It is impossible to find four independent equations using the aforementioned algorithm. In general, in such situations, $$V(x)$$ itself is a rational function and, by setting extra coefficients (higher order ones) to zero, $$W$$ will be computed exactly as $$V$$. When the aforementioned algorithm is used to approximate a function introduced by a PDE equation like OLF and actual function is rational, the significance of Remark 3.1 becomes clear. Remark 3.2 For one-variable functions such as $$V(x):\mathbb{R}\rightarrow\mathbb{R}$$, Padé-like approximation is reduced to conventional Padé approximation for one-variable functions. Remark 3.3 It is worth noting that if the Taylor series coefficients of a function are available, then what remains for this method is solving a set of linear equations. In other words, this method needs solving equations like (3.6) which are linear in unknown variables $$a_i$$s and $$b_i$$s. This equation solving and also all steps of Padé-like approximation algorithm can be implemented by numerical programs such as MATLAB, with no need to symbolic computations. Theorem 3.1 Let $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ be an analytic function in some region about the origin. Suppose that $$W(x):\mathbb{R}^n\rightarrow\mathbb{R}$$ is the approximation of $$V(x)$$ in the form of (3.2) obtained via the PFA algorithm. Then $$V(x)-W(x)=O(\|x\|^p)$$ as $$\|x\|\rightarrow 0$$ for some $$p\geq n_a$$. Proof It may be found at Appendix B that for all monomials like $$x^k$$, the following holds: ∀x∈Rn,|xk|≤‖x‖|k|, (3.8) where $$k$$ is a multi-index. Based on the second part of the PFA algorithm, the expression $$D(x)V(x)-N(x)$$ is free from all monomials $$x^k$$ with $$|k| \leq n_a$$, where $$D(x)$$ and $$N(x)$$ are the denominator and numerator of $$W(x)$$, respectively. Also, some monomials with higher value of $$|k|$$ are annihilated due to the step (5) of the PFA algorithm. Therefore, the remainder may be represented as: D(x)V(x)−N(x)=∑k∈arg(S¯)dkxk, (3.9) where $${\rm arg}(\bar{\mathscr{S}})=\{i \ | \{d_i=0\} \not \in \mathscr{S} \}$$. Applying triangular inequality to (3.9), it is concluded that: |D(x)V(x)−N(x)|≤∑k∈arg(S¯)|dk||xk| and then, by (3.8) one has: |D(x)V(x)−N(x)|≤∑j=j0∞(∑|k|=j|dk|)‖x‖j, (3.10) where $$j_0=\min_{k\in {\rm arg}(\bar{\mathscr{S}})}\{|k|\}$$. Propositions 2.1 and 2.2 state that there are positive constants $$R,C$$ such that: |ck|≤CR|k| and due to (3.7), the following inequality holds: |dk|=|ck+∑|i|=1nbbi.ck−i|≤|ck|+∑|i|=1nb|bi.ck−i|≤CR|k|(1+∑|i|=1nb|bi|.R|i|)=C~R|k|, (3.11) where $$\tilde{C}$$ is defined as: $$\tilde{C}=C\left({1+\sum_{|i|=1}^{n_b}{|b_i|.R^{|i|}}}\right)$$. Since $$n_b$$, $$R$$ and $$b_i$$s are bounded, therefore, $$\tilde{C}$$ is bounded and well-defined. Inequality (3.10) may be simplified using (2.4) and (3.11) as follows: |D(x)V(x)−N(x)|≤∑j=j0∞(n+j−1n−1)C~Rj‖x‖j. (3.12) The recent expression can be rewritten as: ∑j=j0∞(n+j−1n−1)C~Rj‖x‖j=∑j=j0∞(n+j−1n−1)C~2j‖2xR‖j. (3.13) The following inequalities are obvious for all $$r\in \mathbb{R}$$ and $$|r| \leq 1$$: |r|≥|r|2≥|r|3≥⋯ Consequently, it could be concluded that: ∀x∈B1(0),‖x‖1≥‖x‖2≥‖x‖3≥⋯ and also: ∀x∈BR2(0),‖2xR‖1≥‖2xR‖2≥‖2xR‖3≥⋯ (3.14) Relations (3.12)–(3.14) result in the following for all $$x \in B_{R \over 2}(0)$$: |D(x)V(x)−N(x)| ≤∑j=j0∞(n+j−1n−1)C~2j‖2xR‖j≤(∑j=j0∞(n+j−1n−1)C~2j)‖2xR‖j0. (3.15) Now, consider the following generating function for a binomial series (see Graham et al., 1989): 1(1−z)n+1=∑j=0∞(n+jn)zj,∀z∈C,|z|≤1,n∈N. (3.16) Based on (3.16), (3.15) may be reduced to: ∑j=j0∞(n+j−1n−1)C~2j≤C~∑j=0∞(n+j−1n−1)(12)j=C~(1−12)n. (3.17) The relations (3.15) and (3.17) result in: |D(x)V(x)−N(x)|≤C~(1−12)n‖2xR‖j0=C^‖x‖j0, (3.18) where $$\hat{C}={\tilde{C}2^{j_0} \over (1-{1 \over 2})^{n}R^{j_0}}$$. Consider the denominator of approximation (3.2). Due to the fact that $$b_0=1$$, it is resulted that $$D(0)=1$$. Since $$D(x)$$ is a finite term polynomial and consequently is a continuous function, it is possible to select a ball $$B_{\tilde{R}}(0)$$ around the origin such that: minx∈BR~(0)D(x)≥C′, (3.19) where $$0<C'<1$$. Finally, from (3.18) and (3.19), the following relation holds for all $$\|x\| \leq M=\min{\{{R \over 2},\tilde{R}\}}$$. |V(x)−W(x)|=|D(x)V(x)−N(x)D(x)|≤C^C′‖x‖j0. Setting $$p=j_0$$ and $$c_o={\hat{C} \over C'}$$, one completes the proof with regard to (2.5). □ 4. Estimation of domain of attraction As stated in previous sections, the computational backbone of the proposed method for estimation of $$D_A$$ is the OLF introduced in Balint et al. (2009). The OLF is presented by a PDE, whose solution is generally difficult to be found. On the other hand, it is the non-analytical points of OLF that indicate the exact $$D_A$$. The mentioned points produce a ground for using a rational approximation of the OLF. The coefficients of Taylor series expansion of OLF, $$V$$, about the origin can be computed without solving the Equation (2.7), which characterizes the Theorem 2.1. Lemma 4.1 The coefficients of Taylor series expansion of the OLF introduced in Theorem 2.1 can be determined exactly up to any order. Proof See Appendix A. □ Theorem 2.1 and Lemma 4.1 allow one to consider the OLF in the form of (3.1). On the other hand, PFA algorithm gives an approximation in the form of (3.2). So, due to Theorem 3.1, this $$W(x)$$ itself is a Lyapunov function for system (2.6). $$D_A$$ of system (2.6) may be estimated by considering this $$W(x)$$ as a Lyapunov function; i.e. a region provided by Lemma 2.1. It should be noted however, that the main aim of the proposed approach is to consider the roots of denominator of $$W(x)$$ as an approximation of non-analytical points of $$V(x)$$ and, consequently, to consider the roots of denominator of $$W(x)$$ as an estimation of boundary of $$D_A$$ for (2.6). In other words, a region containing the origin and surrounded by the roots of denominator of $$W(x)$$ is an estimate for $$D_A$$. So far, we have not discussed the order of approximation. Indeed, $$n_a$$ and $$n_b$$ remains as the degrees of freedom for the user. Since estimation of $$D_A$$ is determined by denominator of the approximated function, the degree of denominator must be rich enough. For example, for $$n_b=2$$, denominator can only show elliptic, parabolic and hyperbolic shapes. For more complicated shapes, higher order $$n_b$$ is needed. After $$n_b$$ is set, it is preferred to set $$n_a$$ as the minimum value satisfying the following inequality. (n+nan−1)≥(n+nbn), (4.1) where $$n$$ is the number of states of the system. Here, the method is reviewed in the cast of following pseudo algorithm. Padé-like estimation of domain of attraction algorithm: 1. Consider (2.6) as the system equation. 2. Using Equation (2.7) and Lemma 4.1, calculate Taylor expansion coefficients of OLF up to a sufficient order and consider OLF as (3.1). 3. Consider (3.2) as an approximation of OLF. Set $$n_b$$ and $$n_a$$. The relation (4.1) and related argument may be considered as a suggestion. 4. Using the PFA algorithm, determine unknown coefficients $$a_i$$s and $$b_i$$s. 5. Consider the region containing the origin and surrounded by roots of denominator of (3.2) as an estimation of $$D_A$$ for system (2.6). Remark 4.1 As stated in Remark 3.3, if Taylor series coefficients of OLF are available, then the rest of the algorithm can be implemented by numerical programs. Equations (A.12) and (A.15) in the proof of Lemma 4.1 (see Appendix A) show that for some systems, obtaining Taylor series coefficients of OLF can be implemented by numerical programs. Polynomial systems are placed in this class of systems. Actually, (A.12) and (A.15) are linear equations with the elements of $${\partial^i{V(x)} \over \partial{x^j}}{|_{x=0}}$$ as unknown variables. Therefore, if factors of these unknowns ($${\partial^k{f(x)} \over \partial{x^l}}{|_{x=0}}$$s) are available, which is the case for polynomial systems, then $${\partial^i{V(x)} \over \partial{x^j}}{|_{x=0}}$$ can be linearly obtained by a numerical program. Also, this privilege holds for non-polynomial systems for which factors of $${\partial^k{f(x)} \over \partial{x^l}}{|_{x=0}}$$ may be computed numerically, e.g. a system containing terms such as $$\sin$$, $$\cos$$, etc. However, in general, computing Taylor series coefficients of OLF may need some symbolic computations. First, let us apply this method to a simple scalar system in order to illustrate the method. Consider the following scalar system: x˙=−sin⁡(x)−sin⁡(2x). The origin of the system is locally exponentially stable and the set {x∈R | −23π<x<23π≃2.094} is its domain of attraction (Fig. 1 shows $$\dot{x}$$ vs. $$x$$). Now, suppose that we want to approximate OLF $$V(x)$$ by a rational approximation with $$n_a=n_b=2$$ which means: W(x)=a0+a1x+a2x21+b1x+b2x2. Fig. 1. View largeDownload slide Plot of $$\dot{x}=-\sin(x)-\sin(2 x)$$ shows the domain of attraction. Fig. 1. View largeDownload slide Plot of $$\dot{x}=-\sin(x)-\sin(2 x)$$ shows the domain of attraction. It needs the Taylor series expansion of $$V(x)$$. A simple calculation results in: V(x)=16x2+124x4+O(x6). Noting Remark 3.2, one finds that the Padé-like approximation for this function is as the conventional version of the Padé approximation. Therefore, $$W(x)$$ is obtained as: W(x)=16x21−14x2. Roots of denominator of $$W$$ gives an estimate for the $$D_A$$. The estimation reads as the set: {x∈R | −2<x<2}. Now, the efficiency of the proposed method can be shown through some well-known examples. Example 4.1 Consider van der Pol system in the form of: {x˙1=−x2x˙2=x1−x2+x12x2.  $$D_A$$ of this system is surrounded by an unstable limit cycle. This system is a common example for comparison of the existing methods (see Tan & Packard, 2008, Balint et al., 2009). In Padé-like estimation algorithm, some values are first suggested for $$n_a$$ and $$n_b$$, following of which, other possible values may be tried in the line of enhancing the estimation. The results are shown in Fig. 2. For comparison, Fig. 3 simultaneously shows the exact $$D_A$$, results of the proposed method and of some other previous ones. Order of approximation in Fig. 3 for the proposed method is $$n_a=28, \ n_b=6$$ and the the corresponding denominator of the Padé-like estimation can be found in the Appendix C. Fig. 2. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. This simulation shows the improvement of estimation when order of approximation increases. Fig. 2. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. This simulation shows the improvement of estimation when order of approximation increases. Fig. 3. View largeDownload slide The Padé-like estimation result is close to the boundary of exact $$D_A$$ of Example 4.2, which is an unstable limit cycle. Fig. 3. View largeDownload slide The Padé-like estimation result is close to the boundary of exact $$D_A$$ of Example 4.2, which is an unstable limit cycle. Example 4.2 As another system which has a bounded $$D_A$$, the following system can be considered: {x˙1=x2x˙2=−2x1−x2−x13+x1x24+x25.  Different values of $$n_a$$ and $$n_b$$ are tried. Figure 4 shows the result of different approximations based on different $$n_a$$s and $$n_b$$s. Similar to the previous one, this simulation intuitively shows the improvement of the estimation when order of approximation increases. For approximation with $$n_a=28, n_b=6$$, the expression of denominator of approximation is given in Appendix C. Fig. 4. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. Fig. 4. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. Example 4.3 When facing with unbounded $$D_A$$, most estimations methods experience a conservative behaviour. One of the advantages of the proposed method is potential of providing the precise estimation of $$D_A$$, even for unbounded $$D_A$$s. As an example, consider a dynamical system in the form of: {x˙1=x2x˙2=−x1−x2+x13.  The proposed method can be tried for different values of $$n_a$$ and $$n_b$$. The quality of improvement of approximation is shown in Fig. 5. Now, the corresponding result can be compared to those of previous ones (See Fig. 6). It is observed that, the proposed method gives the best estimation in comparison with Lyapunov-based estimation methods (see Fermin et al., 2009) and even some of numerical methods such as trajectory reversing method. Expression of denominator for $$n_a=26,n_b=6$$ can be found in the Appendix C. Fig. 5. View largeDownload slide Estimated $$D_A$$ for Example 4.3 with different values of $$n_a$$ and $$n_b$$. Fig. 5. View largeDownload slide Estimated $$D_A$$ for Example 4.3 with different values of $$n_a$$ and $$n_b$$. Fig. 6. View largeDownload slide The Padé-like estimation technique gives the best result for unbounded $$D_A$$ among different methods. For this simulation, $$n_a=26$$ and $$n_b=6$$ (Example 4.3). Fig. 6. View largeDownload slide The Padé-like estimation technique gives the best result for unbounded $$D_A$$ among different methods. For this simulation, $$n_a=26$$ and $$n_b=6$$ (Example 4.3). Example 4.4 This method can give exact OLF when the OLF itself is rational. Consider a system in the form of: {x˙1=−x1+x12x2x˙2=−x2.  The OLF for this example is obtained as: V(x)=2x22+2x12−x1x232−x1x2. (4.2) For all approximations with $$n_a\geq 4$$ and $$n_b\geq 2$$, and considering Remark 3.1, Equation (4.2) is accurately obtained as a result of utilizing the PFA. Example 4.5 The proposed method can be applied to systems having more than two states. Consider the following system: {x˙1=−x1(1−0.5x12−2x22−3x32)x˙2=−x2(1−0.5x12−2x22−3x32)x˙3=−x3(1−0.5x12−2x22−3x32).  Domain of attraction of this system is the set $$\{(x_1, x_2, x_3)\in \mathbb{R}^3 \ | \ 1-0.5 x_1^2-2x_2^2-3x_3^2>0\}$$. Results of the application of the proposed method to this example results are depicted in Fig. 7. Fig. 7. View largeDownload slide Exact and estimated $$D_A$$ for Example 4.5 with different values of $$n_a$$ and $$n_b$$. Fig. 7. View largeDownload slide Exact and estimated $$D_A$$ for Example 4.5 with different values of $$n_a$$ and $$n_b$$. 5. Stability properties Up to this section, a rough estimation of $$D_A$$ has been obtained, which is sufficient in applications where only approximated borders of $$D_A$$ is needed. One of the major benefits of the proposed method is its ability to present unbounded domains. Also, compared to previous techniques, the geometry and shape of the estimated region in the proposed approach is more similar to those of the exact $$D_A$$. On the other hand, several applications need a guarantee for the estimation region to be inside the stability domain, which is not supported by Theorem 3.1. In this section, through Lemma 2.1, such a guarantee is provided. {However, it is worth noting that for high orders of $$n_a$$ and $$n_b$$, the result of estimation is very close to the exact domain. In the followings, $$\bar{D}_{W,c}$$ means the simple connected region that contains the origin and all of its points satisfy the condition $$W(x)<c$$. The term $$D_{W,c}$$ represents the union of all possible $$\bar{D}_{W,c}$$s. Based on Lemma 2.1, if we find a constant $$c$$ such that the $$D_{W,c}$$ where $$W$$ is a Lyapunov function, falls inside the region where $$\dot{W}<0$$, then $$D_{W,c}$$ is an estimation of $$D_A$$. Obviously, the larger value of $$c$$ gives the larger estimation. Suitable value of $$c$$ may be found via solving a non-linear optimization problem. Theorem 3.1 shows that the function $$W$$ in that theorem is a Lyapunov function. %Therefore, $$W$$ may be used via Lemma 2.1 for a guaranteed estimation of $$D_A$$. Consider the following constrained optimization problem: co= minxW(x)s.t.  W˙(x)=0,x≠0. (5.1) It is obvious that $$\dot{W}$$ is negative in the region $$D_{W,c_o}$$ (excluding the origin), otherwise there exists $$\hat{x}\neq 0$$ such that $$\dot{W}(\hat{x})=0$$ and $$W(\hat{x})=\hat{c}_o<c_o$$, which is in contradiction with the fact that $$c_o$$ is the result of optimization problem. Therefore, the region $$D_{W,c_o}$$ is an estimation of $$D_A$$, which is inside the stability region. The optimization problem (5.1) is a non-linear and (usually) non-convex problem. Since for Lemma 2.1, the global solution of (5.1) is needed, the optimal value of $$c_o$$ is found through a numerical algorithm. Now, let us illustrate the optimization-based method via a simple example. The approximation $$W$$ for Example 4.1 with the order of approximation as $$n_a=n_b=2$$ is as follows: W=32x12−x1x2+x221−14x12+314x1x2−114x22. The optimization problem (5.1) has the global result $$c_o=2.4650786$$ which occurs at points $$(x_1,x_2) = (-0.685248, 0.758646)$$ and $$(x_1,x_2) =(0.758646,-0.685248)$$. See Fig. 8 for details. Fig. 8. View largeDownload slide The region $$D_{W,c}$$ which is inside the region $$\dot{W}=0$$. Fig. 8. View largeDownload slide The region $$D_{W,c}$$ which is inside the region $$\dot{W}=0$$. Figure 9 shows the result of this optimization-based method for Example 4.1. All of these estimations are necessarily inside the stability region. This figure shows that higher order of approximation gives a more precise result.. Fig. 9. View largeDownload slide Result of optimization-based method applied to Example 4.1. The legend shows the orders of numerator and denominator of approximation $$W$$. Fig. 9. View largeDownload slide Result of optimization-based method applied to Example 4.1. The legend shows the orders of numerator and denominator of approximation $$W$$. 6. Conclusion It was found that the estimation of $$D_A$$ based on the proposed, i.e. PFA technique gives more precise results mainly due to the fact that rational polynomials have the potential of revealing the possible non-analytic nature of the functions. Compared to other relevant techniques, the superiority of the proposed approach becomes clear when dealing with unbounded $$D_A$$s. Coefficients of Taylor series expansion of OLF may be computed by a recursive algorithm as stated in the proof of Lemma 4.1. It is worth noting that all computations, including those for coefficients of Taylor series expansion, and those for function approximation, are done in a linear framework. In the meantime, since implementation of the proposed technique needs no symbolic software, it is straightforwardly applicable to high-order systems. Unlike the posed examples, there is no need for system equations, which are supposed to be analytic, to be polynomial. In this paper, no analytic proof for the convergence of roots of denominator of the approximation to the non-analytical points of the main function was given; a subject which is in the line of our future research. Acknowledgements The authors would like to thank members of ACSL (Advanced Control Systems Laboratory) of University of Tehran for the fruitful discussion on the subject. References Baker G. A. ( 1975 ) Essentials of Padé Approximants . New York : Academic Press . Bakhtiari R. & Yazdanpanah M. J. ( 2005 ) Designing a linear controller for polynomial systems with the largest domain of attraction via LMIs. International Conference on Control and Automation, IEEE, Budapest . 1 , 449 – 453 . Balint S. ( 1985 ) Considerations concerning the maneuvering of some physical systems. Ann. Univ. Timisoara . Cambridge, UK : Cambridge Scientific Publishers , 23 , 8 – 16 . Balint S. Braescu L. & Kaslik E. ( 2009 ) Region of Attraction and its Application in Control Theory . Cambridge Scientific Publishers . Benzaouia A. ( 2005 ) Constrained stabilization: an enlargement technique of positively invariant sets. IMA J. Math. Control Inf. Amsterdam : North Holland Publishing Co., 22 , 109 – 118 . Google Scholar CrossRef Search ADS Bruijn N. ( 1961 ) Asymptotic Methods in Analysis . North Holland Publishing Co. Chesi G. ( 2004 ) On the estimation of the domain of attraction for uncertain polynomial systems via LMIs. Proceedings of the 43rd IEEE Conference on Decision Control (CDC) , IEEE, Nassau, Vol.1. pp. 881 – 886 . Chesi G. ( 2009 ) Estimating the domain of attraction for non-polynomial systems via LMI optimizations. Automatica , 45 , 1536 – 1541 , ISSN 0005-1098. Google Scholar CrossRef Search ADS Chesi G. Garulli A. Tesi A. & Vicino A ( 2005 ) LMI-based computation of optimal quadratic Lyapunov functions for odd polynomial systems. Int. J. Robust Nonlinear Control , 15 , 35 – 49 . Google Scholar CrossRef Search ADS Chiang H. D. Hirsch M. & Wu F. ( 1988 ) Stability regions of nonlinear autonomous dynamical systems. IEEE Trans. Autom. Control , 33 , 16 – 27 . ISSN 0018-9286. Google Scholar CrossRef Search ADS Cuyt A. ( 1983 ) The QD-algorithm and multivariate Padé-approximants. Numerische Mathematik , 42 , 259 – 269 . Google Scholar CrossRef Search ADS Evans L. C. ( 1998 ) Partial Differential Equations . Providence, RI : American Mathematical Society . Fermin W. Guerrero- Castellanos J. F. & Alexandrov V. V. ( 2009 ) A computational method for the determination of attraction regions. International Conference on Electrical Engineering, Computing Science and Automatic Control (CCE) . IEEE, Toluca, pp. 1 – 7 . ISBN 978-1-4244-4688-9. Genesio R. ( 1984 ) New techniques for constructing asymptotic stability regions for nonlinear systems. IEEE Trans. Circuits Syst. , 31 , 574 – 581 . Google Scholar CrossRef Search ADS Genesio R. Tartaglia M. & Vicino A. ( 1985 ) On the estimation of asymptotic stability regions: state of the art and new proposals. IEEE Trans. Autom. Control , 30 , 747 – 755 . ISSN 0018-9286. Google Scholar CrossRef Search ADS Graham L. R. Knuth E. D. & Patashnik O. ( 1989 ) Concrete Mathematics, A Foundation for Computer Science . Menlo Park, CA : Addison-Wesley Publishind Company . Guillaume P. ( 1997 ) Nested multivariate Padé approximants. J. Comput. Appl. Math. , 82 , 149 – 158 . ISSN 0377-0427 . Google Scholar CrossRef Search ADS Guillaume P. ( 1998 ) Convergence of the nested multivariate Padé approximants. J. Approx. Theory , 94 , 455 – 466 , ISSN 0021-9045 . Google Scholar CrossRef Search ADS Guillaume P. & Huard A. ( 2000 ) Multivariate Padé approximation. J. Comput. Appl. Math. , 121 , 197 – 219 . ISSN 0377-0427. Google Scholar CrossRef Search ADS Guillaume P. Huard A. & Robin V. ( 1998 ) Generalized multivariate Padé approximants. J. Approx. Theory , 95 , 203 – 214 . ISSN 0021-9045. Google Scholar CrossRef Search ADS Hachicho O. ( 2007 ) A novel LMI-based optimization algorithm for the guaranteed estimation of the domain of attraction using rational Lyapunov functions . J. Franklin Inst. , 344 , 535 – 552 . ISSN 00160032 . Google Scholar CrossRef Search ADS Hashemzadeh F. & Yazdanpanah M. J. ( 2006 ) Semi-global enlargement of domain of attraction for a class of affine nonlinear systems . Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications , IEEE, Munich , Vol. 1 , pp. 2257 – 2262 . Jones R. H. ( 1976 ) General rational approximants in N-variables . J. Approx. Theory , 233 , 201 – 233 . Google Scholar CrossRef Search ADS Kaslik E. Balint A. & Balint S. ( 2003 ) Gradual approximation of the domain of attraction by gradual extension of the “embryo” of the transformed optimal Lyapunov function . Nonlinear Stud. , 10 , 67 – 78 . Khalil H. K. ( 2002 ) Nonlinear Systems , 3rd edn . Englewood Cliffs, NJ: Prentice Hall . Kida S. ( 1990 ) Padé-type and Padé approximants in several variables . Appl. Numer. Math. , 6 , 371 – 391 . ISSN 0168-9274 . Google Scholar CrossRef Search ADS Krantz S. G. & Parks H. R. ( 2002 ) A Primer of Real Analytic Functions . Boston: Birkhauser . ISBN 0817642641 . Google Scholar CrossRef Search ADS Melchor-Aguilar D. & Silviu-Iulian N. ( 2007 ) Estimates of the attraction region for a class of nonlinear time-delay systems . IMA J. Math. Control Inform. , 24 , 523 – 550 . Google Scholar CrossRef Search ADS Steele J. M. ( 2004 ) The Cauchy-Schwarz Master Class . New York: Cambridge University Press . Google Scholar CrossRef Search ADS Tan W. & Packard A. ( 2008 ) Stability region analysis using polynomial and composite polynomial Lyapunov functions and sum-of-squares programming . IEEE Trans. Autom. Control , 53 , 565 – 571 . Google Scholar CrossRef Search ADS Tibken B. ( 2000 ) Estimation of the domain of attraction for polynomial systems via LMIs . Proceedings of the Conference on Decision Control. IEEE, Sydney, NSW , Vol. 4 , pp. 3860 – 3864 . Topcu U. & Packard A. ( 2007 ) Stability region analysis for uncertain nonlinear systems . Proceedings of the IEEE Conference on Decision Control , IEEE, New Orleans, LA , Vol. 1 , pp. 1693 – 1698 . Topcu U. Packard A. Seiler P. & Balas G. ( 2010 ) Robust Region of Attraction Estimation . IEEE Trans. Autom. Control , 55 , 137 – 142 . Google Scholar CrossRef Search ADS Valmórbida G. ( 2009 ) Region of attraction estimates for polynomial systems . Proceedings of the IEEE Conference on Decision Control , IEEE, Shanghai , Vol. 1 , pp. 5947 – 5952 . Vannelli A. & Vidyasagar M. ( 1985 ) Maximal Lyapunov functions and domains of attraction for autonomous nonlinear systems . Automatica , 21 , 69 – 80 . Google Scholar CrossRef Search ADS Vidayasagar M. ( 1993 ) Nonlinear Systems Analysis , 2nd edn . Englewood Cliffs, NJ: Prentice Hall . Yazdanpanah M. J. Khorasani K. & Patel R. V. ( 1998 ) Uncertainty compensation for a flexible-link manipulator using nonlinear $$H_\infty$$ control . Int. J. Control , 69 , 753 – 771 . Google Scholar CrossRef Search ADS Yazdanpanah M. J. Khorasani K. & Patel R. V. ( 1999 ) On the estimate of the domain of validity of non-linear $$H_\infty$$ control . Int. J. Control , 72 , 1097 – 1105 . Google Scholar CrossRef Search ADS Appendix A. Proof of Lemma 4.1 Before proving the lemma, statement of some technical lemmas is needed. Now, the previous lemma is extended for multivariable functions. Now, Lemma 4.1 can be proved. According to Proposition 2.1 for any $$k$$, one has: $$c_{k}=\frac{1}{k!} \left.{\partial^{|k|}{V} \over {\partial{x}^k}}\right|_{x=0}$$. Therefore, it is sufficient to show that $$\left.{\partial^{|k|}{V} \over {\partial{x}^k}}\right|_{x=0}$$ may be computed for any $$k\in (\mathbb{N}\cup\{0\})^n$$. An induction can be established on $$|k|$$ as follows. Start with $$|k|=1$$. Using Equation (2.4), there are $${n\choose n-1}=n$$ derivatives of order $$|k|=1$$. So, for $$l=1,2,\dots,n$$, two sides of Equation (2.7) are driven as: ∂∂xl(⟨∇V,f⟩)=∂∂xl(−‖x‖2). (A.8) Equation (A.8) can be expanded as follows: ∂∂xl(∑j=1n∂V(x)∂xj.fj(x))=−2∑j=1nxj. (A.9) Thereupon: (∑j=1n∂2V(x)∂xjxl.fj(x)+∂V(x)∂xj.∂∂xlfj(x))=−2∑j=1nxj. (A.10) Evaluation of Equation (A.10) at origin and noting that $$f(0)=0$$ results in: ∑j=1n( ∂V(x)∂xj|x=0. ∂fj(x)∂xl|x=0)=0. (A.11) Equation (A.11) is an equation in which the derivatives of $$V$$ are unknowns. Considering (A.11) for $$l=1,2,\dots,n$$, a system with $$n$$ linear equation in $$n$$ unknowns is obtained. They can be rewritten in the following matrix form: [∂f1(x)∂x1∂f2(x)∂x1…∂fn(x)∂x1∂f1(x)∂x2∂f2(x)∂x2…∂fn(x)∂x2⋮⋮⋱⋮∂f1(x)∂xn∂f2(x)∂xn…∂fn(x)∂xn]x=0[∂V(x)∂x1∂V(x)∂x2⋮∂V(x)∂xn]x=0=[00⋮0]. (A.12) Coefficients matrix in (A.12) is equal to $${\partial{f} \over \partial{x}}^T|_{x=0}$$. Since, based on the assumption of Theorem 2.1, the matrix $${\partial{f} \over \partial{x}}|_{x=0}$$ is Hurwitz; therefore, it is invertible. Invertibility of coefficients matrix in (A.12) means that this equation has a unique solution and this solution is:$${\partial{V(x)} \over \partial{x_l}}|_{x=0}=0$$ for $$l=1,2,\dots,n$$. This result was expectable due to positive definiteness of $$V(x)$$. Assume that for all $$|k|=\mathfrak{n-1}$$, $${\partial^{|k|} \over \partial{x^k}}V(x)|_{x=0}$$ can be computed. Set $$|k|=\mathfrak{n}$$. For all $$|k|=\mathfrak{n}$$, the equations which contain derivatives of order $$\mathfrak{n}$$ of $$V(x)$$ as unknowns are obtainable through: ∂|k|∂xk(⟨∇V,f⟩)=∂|k|∂xk(−‖x‖2). (A.13) The left-hand side of (A.13) can also be expanded as follows: ∂|k|∂xk(∑j=1n∂V(x)∂xj.fj(x)) =∑j=1n∂|k|∂xk(∂V(x)∂xj.fj(x)) =∑j=1n(∑i=0k(ki)(∂|k|−|i|∂xk−i∂V(x)∂xj).(∂|i|∂xifj(x))) =∂|k|∂xk(−‖x‖2). (A.14) The right-hand side of the last equality in (A.14) for $$\mathfrak{n}>2$$ is equal to zero. Also, for $$\mathfrak{n}=2$$, when one of the elements of $$k$$ is $$2$$ such as $$k=(0,0,\dots,2,\dots,0)$$, then, the right-hand side of the last equality in (A.14) is equal to $$-2$$; otherwise, it is like $$k=(0,\dots,1,0,\dots,1,\dots,0)$$, which is equal to zero. However, due to $$f(0)=0$$, (A.14) can be rewritten as:  ∑j=1n(∑|i|=1(ki)(∂|k|−|i|∂xk−i∂V(x)∂xj).(∂|i|∂xifj(x))) =∗∗∗ (A.15) The $$***$$ on the right-hand side of the above equation represents all derivatives of $$V(x)$$ whose order are lower than $$\mathfrak{n}$$, higher order derivatives of $$f$$ and constant terms. When considering $$|k|=\mathfrak{n}$$, all terms in $$***$$ are known and do not affect solvability of (A.15). These equations are linearly independent and consequently are solvable. Actually, the matrix of coefficients of above equations is a block diagonal matrix whose blocks are leading principle minors of $${\partial{f} \over \partial{x}}|_{x=0}$$. Therefore, it is invertible and consequently the mentioned equations are solvable. Appendix B. Proof of (3.8) First, consider the following lemma from Steele (2004). In order to prove (3.8), it is sufficient to show that: (|xk|)2≤(‖x‖|k|)2 or equivalently: (|xk|)2|k|≤(‖x‖)2. (B.2) Inequality (B.2) may be expanded as follows: (x12)k1|k|(x22)k2|k|…(xn2)kn|k|≤(x12+x22+⋯+xn2). (B.3) The left-hand side of (B.3) is in the form of (B.1) with $$x_i^2$$ as $$a_i$$ and $${k_i \over |k|}$$ as $$p_i$$ and $${k_1 \over |k|}+{k_2 \over |k|}+\dots+{k_n \over |k|}={|k| \over |k|}=1$$. Therefore, (B.1) states that: (x12)k1|k|(x22)k2|k|…(xn2)kn|k| ≤k1|k|.x12+k2|k|.x22+⋯+kn|k|.xn2 ≤x12+x22+⋯+xnn and the last inequality holds because $${k_i \over |k|}\leq 1$$. Appendix C. Some simulation expressions Expression of denominator of approximation for Example 4.1 with $$n_a=28,n_b=6$$: $$1.0+ 0.1065412730\,x_{{1}}x_{{2}}+ 0.1362438857\,{x_{{1}}}^{3}x_{{2}} - 0.2242705150\,{x_{{1}}}^{2}{x_{{2}}}^{2}+ 0.1013585171\,x_{{1}}{x_{{2}}}^{3}- 0.04151182324\,{x_{{1}}}^{5}x_{{2}}+ 0.04835446107\,{x_{{1}}}^{4}{x_{{2}}}^{2}- 0.02548402153\,{x_{{1}}}^{3}{x_{{2}}}^{3}+0.009070713464\,{x_{{1}}}^{2}{x_{{2}}}^{4}\\- 0.002330734027\,x_{{1}}{x _{{2}}}^{5}- 0.1480416345\,{x_{{1}}}^{2}- 0.1098080759\,{x_{{2}}}^{2}-0.01336682711\,{x_{{1}}}^{4}\\- 0.01896447641\,{x_{{2}}}^{4}-0.003112455890\,{x_{{1}}}^{6}+ 0.0002849304531\,{x_{{2}}}^{6}$$ Expression of denominator of approximation for Example 4.1 with $$n_a=28,n_b=6$$: $$1.0- 0.2790014167\,x_{{1}}x_{{2}}+ 0.2368419702\,{x_{{1}}}^{3}x_{{2}}- 0.7081714653\,{x_{{1}}}^{2}{x_{{2}}}^{2}+ 0.3648060899\,x_{{1}}{x_{{ 2}}}^{3}+ 0.1404756418\,{x_{{1}}}^{5}x_{{2}}+ 0.05357554804\,{x_{{1}}}^{4}{x_{{2}}}^{2}- 0.05682371974\,{x_{{1}}}^{3}{x_{{2}}}^{3}+0.5660169893\,{x_{{1}}}^{2}{x_{{2}}}^{4}+ 0.1787257727\,x_{{1}}{x_{{2}}}^{5}+ 0.07123538082\,{x_{{2}}}^{6}- 0.8135596832\,{x_{{1}}}^{2}- 0.05029478641\,{x_{{2}}}^{2}\\- 0.5301539075\,{x_{{1}}}^{4}- 0.7668262956\,{x_{{2}}}^{4}+ 0.005283318119\,{x_{{1}}}^{6}$$ Expression of denominator of approximation for Example 4.3 with $$n_a=26, n_b=6$$: $$1.0- 0.002314998116\,{x_1}^{4}+ 0.001060978122\,{x_2}^{4}- 0.002242078600\,{x_1}^{6}- \\ 0.000002985127786\,{x_2}^{6}- 0.2781498140\,{x_1}^{2}- 0.1527167633\,{x_2}^{2}+\\ 0.004763572498\,x_1{x_2}^{3}- 0.001975789673\,{x_1}^{5}x_2- 0.0009491287752\,{x_1}^{4}{x_2}^{2}-\\ 0.0003555944435\,{x_1}^{3}{x_2}^{3}- 0.0001177270691\,{x_1}^{2}{x_2}^{4}- \\ 0.00002591183460\,x_1{x_2}^{5}- 0.01080307989\,{x_1}^{3}x_2+ 0.005802793836\,{x_1}^{2}{x_2}^{2}- 0.3086406808\,x_1x_2$$ © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Mathematical Control and Information Oxford University Press

# Padé-like approximation and its application in domain of attraction estimation

IMA Journal of Mathematical Control and Information, Volume Advance Article (2) – Jan 11, 2017
27 pages

Publisher
Oxford University Press
ISSN
0265-0754
eISSN
1471-6887
D.O.I.
10.1093/imamci/dnw071
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this paper, a rational approximation of multivariable functions is proposed. Since the proposed approach, hereafter called Padé-like approximation, uses the coefficients of Taylor series expansion at a point, it can be considered as a novel extension of conventional Padé approximation for multivariable functions. This extension is free from some of previous extensions imperfections. Similar to the conventional Padé approximation for one-variable functions, the outcome of the proposed method can reveal the points where the multi-variable function is non-analytic. The proposed method, whose accuracy is analytically guaranteed, could find a wide range of applications wherever an accurate approximation is needed. In the line of estimating domain of attraction in non-linear systems, Padé-like approximation of an optimal Lyapunov function (OLF) is obtained from the corresponding partial differential equation. Using the approximate OLF, an estimation of domain of attraction is obtained which is less conservative compared to those based on the previous techniques. 1. Introduction For one-variable functions, the Padé approximation is a well-known method and is considered as the best rational function approximation (see Baker, 1975). Padé approximation could be valid even in the region where the Taylor series does not converge (see Guillaume & Huard, 2000). Due to the rationality of approximation, it can be used for approximating the non-analytic points of a function. Also, construction by Taylor series coefficients is another great property of Padé approximation. That is because, in some cases, Taylor series coefficients of an even non-explicitly known function may be obtained from the equation (possibly a partial differential equation, PDE) satisfying it. It should be noted that approximate solution of an equation is usually different from the expansion of the exact solution. The aforementioned property of the Padé approximation, however, causes no difference between the equation and the function approximations. It is important to note that for multi-variable functions, the matter is not so simple, which is emergence of the fact that the extension of Padé approximation to multi-variable functions is not straightforward. There are some extensions of Padé approximation for multi-variable functions such as those posed in Jones (1976), Cuyt (1983), Guillaume (1997, 1998), Guillaume et al. (1998), Kida (1990), Guillaume & Huard (2000). With regard to the estimation of domain of attraction, none of the mentioned methods are sufficient. Methods like the one in Jones (1976) need some symmetry property in coefficients of Taylor series and do not work for some functions, especially those vanish at the origin. On the other hand, methods like the one posed in Cuyt (1983) give a high-order polynomial, even for a low order approximation, and in the meantime, the origin is always one of the roots of approximation denominator. The recent fact makes this method useless for those purposes which need the study of non-analyticity property of functions. This matter becomes important for our application, i.e. the computation of $$D_A$$ (domain of attraction), which may be related to non-analytic points of a multi-variable function (see Balint et al., 2009). Guillaume (1997, 1998) has proposed a nested Padé approximation technique which suffers from non-homogeneity of approximation degree for each variable. In the meantime, the denominator is not a polynomial and furthermore, it could cause some computational costs. Guillaume et al. (1998) have used the least square technique for obtaining the coefficients of numerator and denominator of the rational approximation. Through this method, Taylor series expansion coefficients of the difference of a function and its approximation become marginal, but not zero. These circumstances convinced the authors to develop a novel rational approximation and use its denominators roots as an estimation of domain of attraction. By computation of $$D_A$$, we mean the solution of the problem of finding initial conditions which converge to an asymptotic stable equilibrium point. A linear model is usually an approximation of a real non-linear system about its operating point. Hence, the validity of all subsequent analyses and designs is local. As a result, after stabilization, one has to determine the corresponding domain of attraction (see Yazdanpanah et al., 1999). From control engineering point of view, the importance of finding the $$D_A$$ seems obvious in analysis and synthesis. From application point of view, computation of $$D_A$$ has found numerous applications in different fields of science and technology. Power plant systems, chemical reactors, aircrafts and robotic systems are some typical examples of $$D_A$$ applications in engineering. Biological systems, economy and ecology are some other applied fields where $$D_A$$ plays a crucial role in the ground of behavioural analysis (see Genesio et al., 1985, Yazdanpanah et al., 1998). Moreover, enlargement of the closed loop $$D_A$$ has been a goal for some controller design strategies (see Bakhtiari & Yazdanpanah, 2005, Benzaouia, 2005, Hashemzadeh & Yazdanpanah, 2006). $$D_A$$ has rich and challenging theoretical aspects as it has a practically important and beneficial insight. It is clear that $$D_A$$ is a complicated set and, except very special cases, it has no analytical and simple representation. Therefore, massive research effort has been done on the $$D_A$$ estimation during past decades. $$D_A$$ estimation is extensively considered in the literature and many different methods have been applied for it. These methods can be divided into two general categories: Lyapunov-based methods and non-Lyapunov-based methods (see Genesio et al., 1985). This classification is very common and is almost used by every paper in the field. Lyapunov-based methods try to find a Lyapunov function based on which an invariant set is proposed as the estimation of $$D_A$$. Trajectory reversing method (see Genesio, 1984) and topological-based method (see Chiang et al., 1988) can be mentioned as some non-Lyapunov-based methods. Genesio et al. (1985) have reviewed the relevant works published till 1985. Papers related to the subject have reviewed briefly these methods in their introduction sections. Zubov theorem and the estimation methods based on that are early approaches in this field (see Genesio et al., 1985, Melchor-Aguilar & Niculescu, 2007 and Balint et al., 2009 for a class of time-delay systems). Vannelli & Vidyasagar (1985) introduced maximal Lyapunov function and proposed an estimation method which is categorized under Lyapunov-based methods. This method is one of the few ones which uses a rational Lyapunov function. Hachicho (2007) has extended this method through an optimization technique with LMI (linear matrix inequality) constraints. An important branch of Lyapunov-based method which is concentrated in the past decade is optimization-based methods. These methods try to formulate the finding of a Lyapunov function which maximizes $$D_A$$ estimation—towards its exact domain—as an optimization problem. This optimization problem, in general, is a non-convex problem (see Tibken, 2000). These methods are applied to odd polynomial systems (see Chesi & Garulli, 2005), polynomial systems (see Tibken, 2000, Chesi, 2004, Valmórbida, 2009) and some non-polynomial systems (see Chesi, 2009). Square programming is also a new and powerful optimization-based method that is used for $$D_A$$ estimation (see Topcu & Packard, 2007, Tan & Packard, 2008, Topcu et al., 2010). Although $$D_A$$ estimation is extensively studied, almost all of the obtained results are conservative. For example, the results regarding unbounded $$D_A$$s are not satisfactory. Consequently, $$D_A$$ estimation is still an open problem for research. Balint et al. have presented one of the interesting theorems in this area which presents the conditions for exact determination of $$D_A$$ (see Balint, 1985, Kaslik et al., 2003, Balint et al., 2009). Simply, Balint et al. (2009) states that exact $$D_A$$ for a system is coincides with the domain of analyticity of a special Lyapunov function. This Lyapunov function which is introduced by a PDE is called the optimal Lyapunov function (OLF). It is worth noting that such an OLF is different from the similarly called function in the Zubov theore. Like many other PDEs, obtaining the closed-form solution of the OLF—but in a very special cases—generally is not possible. Furthermore, for this type of PDEs, solving PDE by characteristic method (Evans, 1998)—a known analytical way—is equivalent to solving the dynamical equations of the system. Hence, finding the approximate solution for this PDE becomes significant. Rational polynomials have the potential of revealing the possible non-analytic nature of the functions. Consequently, it seems that using a rational function as an approximation of the OLF can approximate non-analytical points of OLF based on which an estimation of $$D_A$$ may be obtained. This fact is the linkage key of Padé-like approximation and estimation of $$D_A$$. In the line of aforementioned justification, here, first, a novel rational approximation is developed for multi-variable functions. This approximation, like conventional Padé approximation, uses Taylor series coefficients and preserves the great property of conventional Padé approximation discussed above. This method is quite different from previous ones and resolves some of their deficiencies. In this paper, this method is used for the purpose of estimation of $$D_A$$, but like other approximators, has a wide potential applications. The convergence of this method is shown mathematically in Theorem 3.1. In the next step, this approximation method is applied to approximate an OLF. The roots of denominator of this rational approximation are near to the non-analytical points of OLF. Thus, the denominator of this approximation is used to present an estimation of $$D_A$$. Applying this method to some known sample systems shows its power, especially for unbounded $$D_A$$s. Some applications of $$D_A$$ estimation need a guarantee for the estimated region to be inside the exact domain. One possible way to present this guarantee is to prove that the roots of approximation denominator converge to the non-analytical points of the OLF. This paper does not pursue the mentioned line of proof because of the following reason. Usually, convergence theorems of approximation methods discuss the differences between approximate and exact solutions for a sufficiently high order of approximation, while a high order approximation may be useless due to the computational burden. For such cases in the estimation process, instead of studying the convergence of roots of approximation denominator to the non-analytical points of OLF, we propose a method in Section 5 to put the estimated domain inside the exact domain. However, no analytical proof is given for the convergence of roots of denominator of the approximation to the non-analytical points of the OLF. This paper is organized as follows. Section 2 is dedicated to preliminary materials. Constructive introduction to Padé-like function approximation (PFA) is considered in Section 3. The proposed method for $$D_A$$ estimation is presented in Section 4. Some stability analysis for estimation result is posed in Section 5. Finally, this paper is ended with conclusion in Section 6. 2. Preliminary The following nomenclature is used frequently in the sequel. A multi-index $$k$$ is an element of $$({\mathbb{N} \cup \{0\}})^n$$ or loosely, $$k=(k_1,k_2,\dots,k_n)$$. Some standard multi-index notations with $$x \in \mathbb{R}^n, \quad i,k \in ({\mathbb{N} \cup \{0\}})^n, \quad n\in \mathbb{N}$$ are as follows (see Krantz & Parks, 2002): |k| =k1+k2+⋯+kn,k! =k1!k2!…kn!,i<kifij<kj∀j=1,2,…,n,i≤kifij≤kj∀j=1,2,…,n,(ki) =(k1i1)(k2i2)…(knin),xk =x1k1x2k2…xnkn,∂|k|∂xk =∂k1∂xk1∂k2∂xk2…∂kn∂xkn. (2.1) Also, the following notation can simplify the expressions: ci =c(i), c:(N∪{0})n→R∑i=0k(.) ≡∑i1=0k1∑i2=0k2…∑in=0kn(.)∑|i|=m∈N(.) ≡∑i1=0m∑i2=0m−i1…∑in=m−i1−i2...−in−1(.)∑|i|=0m∈N(.) ≡∑i1=0m∑i2=0m−i1…∑in=0m−∑j=0n−1ij(.). (2.2) It is convenient to set a convention: ci=0,∀i∉(N∪{0})n. (2.3) The expression like $$x^k$$ where $$x\in \mathbb{R}^n,k\in ({\mathbb{N} \cup \{0\}})^n$$ is called a monomial of order $$|k|$$. It is straight to show that there are (n+|k|−1n−1)=(n+|k|−1)!(|k|)!(n−1)! (2.4)$$n$$-variable monomials of order $$|k|$$. Also, there are (n+|k|n)$$n$$-variable monomials of order less than or equal to $$|k|$$. Similarly, if $$V(x)$$ is at least a $$|k|$$-times continuously differentiable function, then it has $${n+|k|-1 \choose n-1}$$ derivatives of order $$|k|$$ such as: ∂|k|V∂xk=∂|k|V∂k1x1,∂k2x2,…,∂knxn. Throughout this paper, $$\|.\|$$ is used for 2-norm, i.e.: ∀x∈Rn,‖x‖=x12+x22+⋯+xn2 and a ball centred at the point $$x_0$$ with radius $$\epsilon \in \mathbb{R}_{\geq 0}$$ is considered as follows: Bϵ(x0)={x∈Rn|‖x−x0‖≤ϵ}. A function $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ is said to be real analytic at point $$\alpha\in\mathbb{R}^n$$ if there exist a neighbour of $$\alpha$$ such as $$U\subset\mathbb{R}^n$$ and a convergent power series such that (see Krantz & Parks, 2002): ∀x∈U;V(x)=∑|k|=0∞ck(x−α)k,ck∈R. This series representation is called Taylor series expansion about $$\alpha$$. Proposition 2.1 Let $$V$$ be a real analytic function at point $$\alpha$$. Then, its Taylor series coefficients satisfy the following (see Krantz & Parks, 2002): ck=1k! ∂|k|V∂xk|x=0. Proposition 2.2 Function $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ is real analytic on an open set $$U\subset \mathbb{R}^n$$ if and only if for each $$\alpha \in U$$ there are an open ball $$B_{\epsilon}(\alpha) \subset U$$ and constants $$C>0,R>0$$ such that (see Krantz & Parks, 2002): ∀x∈Bϵ(α),|∂|k|V∂xk|≤Ck!R|k|. For functions $$f,\phi: \mathfrak{S}\rightarrow \mathbb{R}$$, it is said that $$f(x)=O(\phi(x))$$ if there is a positive constant $$c_o$$ such that (see Bruijn, 1961): |f(x)|≤co|ϕ(x)|∀x∈S. Consequently, for function $$f(x):\mathbb{R}^n \rightarrow \mathbb{R}$$, it is said that $$f(x)=O\left(\|x\|^p\right)$$ as $$\|x\| \rightarrow 0$$ if there exist constants $$M$$ and $$c_o$$ such that: |f(x)|≤co‖x‖p∀‖x‖≤M. (2.5) In this paper, the definition of domain of attraction is considered as follows. Definition 2.1 (Khalil, 2002) Let $$x=0$$ be an asymptotic stable equilibrium point of: x˙=f(x),f:D⊆Rn→Rn, (2.6) where $$f$$ is locally Lipschitz and $$D$$ contains the origin. Let $$\phi(x;t)$$ be the solution of this system at time $$t$$ with initial condition $$x$$. The domain of attraction $$(D_A)$$ is defined as: DA={x∈D|ϕ(x;t)→0 as t→+∞}. Here, the theorem which relates the exact $$D_A$$ to the domain of analyticity of a special Lyapunov function is reviewed (see Balint, 1985; Balint et al., 2009). Theorem 2.1 (Balint et al., 2009) Consider system (2.6). Let $$f$$ be real analytic and $${\partial{f} \over \partial{x}}|_{x=0}$$ be Hurwitz. Then, domain of attraction of $$x=0$$ coincides with the domain of analyticity of the real analytical function $$V$$ defined by: ⟨∇V,f⟩=−‖x‖2V(0)=0. (2.7)$$V$$ is strictly positive on $$D_A-\{0\}$$ and $$V(x)\rightarrow\infty$$ as $$x\rightarrow\partial D_A$$ or $$\|x\|\rightarrow\infty$$. The above unique function V is called OLF. The following lemma which is based on Lyapunov theory is the foundation of several estimation methods. Lemma 2.1 (Vidayasagar, 1993) Assume there exists a smooth function $$V:\mathbb{R}^n \rightarrow \mathbb{R}$$ and a positive constant $$c$$ for system (2.6) such that: V(0)=0,V(x)>0∀x∈Rn−0 and ΩV,c:={x∈Rn|V(x)≤c} is bounded. Now, if the condition: ΩV,c⊂{x∈Rn|∇V(x).f(x)≤0} is satisfied then $$\it{\Omega}_{V,c}$$ is an invariant subset of $$D_A$$. 3. PFA Despite the existence of Padé approximation for one-variable functions, the current extensions for multivariable functions are not effectual in the line of this work’s goal. So, here, an attempt is made to propose a new extension of Padé approximation which is suited for OLF approximation. This extension like the conventional Padé approximation is built from Taylor series coefficients. So, in this section, it is assumed that the expansion of $$V$$ in the form: V(x)=∑|i|=0∞cixi (3.1) is in hand, and then this function is tried to be approximated via an appropriate rational function. Let, W(x)=∑|i|=0naaixi∑|i|=0nbbixi (3.2) be the required rational approximation. The parameter $$b_0$$ is set to $$1$$ for uniqueness of approximation. Now, in some manner, coefficients $$a_i$$s and $$b_i$$s must be determined. According to the previous argument, there are $${n+n_a} \choose {n}$$ and $${{n+n_b} \choose {n}}-1$$ unknown coefficients for numerator and denominator of (3.2), respectively. In this regard, (n+nan)+(n+nbn)−1 (3.3) independent equations are needed. $$a_i$$s and $$b_i$$s must be determined so that: V(x)−W(x)=0. (3.4) However, it should be noted that, it is usually impossible to annihilate the right-hand side of the above equation; so, $$a_i$$s and $$b_i$$s must be determined such that the right-hand side of (3.4) is possibly free of lower order terms. In this connection, consider: W(x)−V(x)=∑|i|=0naaixi∑|i|=0nbbixi−V(x)=0⇒(∑|i|=0naaixi)−(∑|i|=0nbbixi)V(x)=0 or equivalently: (∑|i|=0naaixi)−(∑|i|=0nbbixi)(∑|i|=0∞cixi)=0. (3.5) The above equation may be rearranged as: ∑|i|=0na(ai−di)xi−∑|i|=na+1∞dixi=0, (3.6) where di={∑j1=0i1∑j2=0i2…∑jn=0inbj.ci−jif i≤nb∑|j|=0nbbj.ci−jif i>nb.  With convention (2.3), we can consider $$d_i$$ for all $$i$$ as: di=∑|j|=0nbbj.ci−j. Since $$b_0=1$$, one has: di=ci+∑|j|=1nbbj.ci−j. (3.7) By setting, $$a_i-d_i=0 \quad \text{for} \quad 0\leq |i| \leq n_a$$, the first sum in Equation (3.6) will disappear. Due to the coefficients of $$a_i$$s, these sets of equations are linearly independent and consequently solvable. Actually, if $$b_i$$s are known, then equations $$a_i-d_i=0$$ determine $$a_i$$s. So, further $${{n+n_b} \choose {n}}-1$$ independent equations must be chosen from the second sum in Equation (3.6). Considering Equation (3.7), when $$|i|$$ is fixed, it is better to choose $$d_i=0$$, whose corresponding $$c_i$$ is greater in the absolute value sense. On the other hand, lower powers of $$|i|$$ have preference. So, at first, $$|i|=n_a+1$$ is considered and due to the descending order of $$c_i$$s’ absolute value, the corresponding $$d_i=0$$ is added to the equation set such that the recent equation is independent and consistent from the previous ones. Henceforth, if the number of equations is inadequate, this operation is iterated with the next of $$|i|$$. This operation continues until an enough number of equations are obtained. The method can be summarized in the following pseudo algorithm. PFA algorithm: 1. Consider (3.2) as a rational approximation of (3.1). The parameters $$a_i$$s and $$b_i$$s are unknown variables which must be determined. In the sequel, $$\mathscr{S}$$ represents the set of equations needed for finding the parameters $$a_i$$s and $$b_i$$s. 2. For $$\quad 0\leq|i|\leq n_a$$, add $$a_i-d_i=0$$ to the equation set $$\mathscr{S}.$$ 3. $$k=n_a+1$$. 4. For $$|i|=k$$, arrange $$c_i$$ from largest to the smallest absolute value. Consider $$c_j, j \in \mathbb{N}$$, where $$j=1,2,\dots,{{n-k+1} \choose {k}}$$ shows the mentioned descending order. (a) $$j=1$$. (b) Consider $$c_i$$ and $$d_i$$ corresponding to $$c_j$$. If $$d_i=0$$ is independent and consistent to the equation set $$\mathscr{S}$$; add it to $$\mathscr{S}$$. (c) If number of equations in set $$\mathscr{S}$$ is not equal to (3.3) and $$j<{{n-k+1} \choose {k}}$$, then $$j \leftarrow j+1$$ and go to 4(b). Else, go to (5). 5. If number of equations in set $$\mathscr{S}$$ is not equal to (3.3), then $$k \leftarrow k+1$$ and go to (5). 6. Solve the linear equation set $$\mathscr{S}$$ for unknown variables $$a_i$$s and $$b_i$$s. Remark 3.1 If enough equations cannot be found in steps (5) and (5) in the above algorithm, it can be concluded that the order of approximation is higher than the main function. As an exaggerated example, consider $$V(x)=\frac{1}{1-x}=\sum_{i=0}^{\infty}{x^i}$$ and $$W(x)=\frac{a_1x+a_0}{b_2x^2+b_1x+1}$$. It is impossible to find four independent equations using the aforementioned algorithm. In general, in such situations, $$V(x)$$ itself is a rational function and, by setting extra coefficients (higher order ones) to zero, $$W$$ will be computed exactly as $$V$$. When the aforementioned algorithm is used to approximate a function introduced by a PDE equation like OLF and actual function is rational, the significance of Remark 3.1 becomes clear. Remark 3.2 For one-variable functions such as $$V(x):\mathbb{R}\rightarrow\mathbb{R}$$, Padé-like approximation is reduced to conventional Padé approximation for one-variable functions. Remark 3.3 It is worth noting that if the Taylor series coefficients of a function are available, then what remains for this method is solving a set of linear equations. In other words, this method needs solving equations like (3.6) which are linear in unknown variables $$a_i$$s and $$b_i$$s. This equation solving and also all steps of Padé-like approximation algorithm can be implemented by numerical programs such as MATLAB, with no need to symbolic computations. Theorem 3.1 Let $$V:\mathbb{R}^n\rightarrow\mathbb{R}$$ be an analytic function in some region about the origin. Suppose that $$W(x):\mathbb{R}^n\rightarrow\mathbb{R}$$ is the approximation of $$V(x)$$ in the form of (3.2) obtained via the PFA algorithm. Then $$V(x)-W(x)=O(\|x\|^p)$$ as $$\|x\|\rightarrow 0$$ for some $$p\geq n_a$$. Proof It may be found at Appendix B that for all monomials like $$x^k$$, the following holds: ∀x∈Rn,|xk|≤‖x‖|k|, (3.8) where $$k$$ is a multi-index. Based on the second part of the PFA algorithm, the expression $$D(x)V(x)-N(x)$$ is free from all monomials $$x^k$$ with $$|k| \leq n_a$$, where $$D(x)$$ and $$N(x)$$ are the denominator and numerator of $$W(x)$$, respectively. Also, some monomials with higher value of $$|k|$$ are annihilated due to the step (5) of the PFA algorithm. Therefore, the remainder may be represented as: D(x)V(x)−N(x)=∑k∈arg(S¯)dkxk, (3.9) where $${\rm arg}(\bar{\mathscr{S}})=\{i \ | \{d_i=0\} \not \in \mathscr{S} \}$$. Applying triangular inequality to (3.9), it is concluded that: |D(x)V(x)−N(x)|≤∑k∈arg(S¯)|dk||xk| and then, by (3.8) one has: |D(x)V(x)−N(x)|≤∑j=j0∞(∑|k|=j|dk|)‖x‖j, (3.10) where $$j_0=\min_{k\in {\rm arg}(\bar{\mathscr{S}})}\{|k|\}$$. Propositions 2.1 and 2.2 state that there are positive constants $$R,C$$ such that: |ck|≤CR|k| and due to (3.7), the following inequality holds: |dk|=|ck+∑|i|=1nbbi.ck−i|≤|ck|+∑|i|=1nb|bi.ck−i|≤CR|k|(1+∑|i|=1nb|bi|.R|i|)=C~R|k|, (3.11) where $$\tilde{C}$$ is defined as: $$\tilde{C}=C\left({1+\sum_{|i|=1}^{n_b}{|b_i|.R^{|i|}}}\right)$$. Since $$n_b$$, $$R$$ and $$b_i$$s are bounded, therefore, $$\tilde{C}$$ is bounded and well-defined. Inequality (3.10) may be simplified using (2.4) and (3.11) as follows: |D(x)V(x)−N(x)|≤∑j=j0∞(n+j−1n−1)C~Rj‖x‖j. (3.12) The recent expression can be rewritten as: ∑j=j0∞(n+j−1n−1)C~Rj‖x‖j=∑j=j0∞(n+j−1n−1)C~2j‖2xR‖j. (3.13) The following inequalities are obvious for all $$r\in \mathbb{R}$$ and $$|r| \leq 1$$: |r|≥|r|2≥|r|3≥⋯ Consequently, it could be concluded that: ∀x∈B1(0),‖x‖1≥‖x‖2≥‖x‖3≥⋯ and also: ∀x∈BR2(0),‖2xR‖1≥‖2xR‖2≥‖2xR‖3≥⋯ (3.14) Relations (3.12)–(3.14) result in the following for all $$x \in B_{R \over 2}(0)$$: |D(x)V(x)−N(x)| ≤∑j=j0∞(n+j−1n−1)C~2j‖2xR‖j≤(∑j=j0∞(n+j−1n−1)C~2j)‖2xR‖j0. (3.15) Now, consider the following generating function for a binomial series (see Graham et al., 1989): 1(1−z)n+1=∑j=0∞(n+jn)zj,∀z∈C,|z|≤1,n∈N. (3.16) Based on (3.16), (3.15) may be reduced to: ∑j=j0∞(n+j−1n−1)C~2j≤C~∑j=0∞(n+j−1n−1)(12)j=C~(1−12)n. (3.17) The relations (3.15) and (3.17) result in: |D(x)V(x)−N(x)|≤C~(1−12)n‖2xR‖j0=C^‖x‖j0, (3.18) where $$\hat{C}={\tilde{C}2^{j_0} \over (1-{1 \over 2})^{n}R^{j_0}}$$. Consider the denominator of approximation (3.2). Due to the fact that $$b_0=1$$, it is resulted that $$D(0)=1$$. Since $$D(x)$$ is a finite term polynomial and consequently is a continuous function, it is possible to select a ball $$B_{\tilde{R}}(0)$$ around the origin such that: minx∈BR~(0)D(x)≥C′, (3.19) where $$0<C'<1$$. Finally, from (3.18) and (3.19), the following relation holds for all $$\|x\| \leq M=\min{\{{R \over 2},\tilde{R}\}}$$. |V(x)−W(x)|=|D(x)V(x)−N(x)D(x)|≤C^C′‖x‖j0. Setting $$p=j_0$$ and $$c_o={\hat{C} \over C'}$$, one completes the proof with regard to (2.5). □ 4. Estimation of domain of attraction As stated in previous sections, the computational backbone of the proposed method for estimation of $$D_A$$ is the OLF introduced in Balint et al. (2009). The OLF is presented by a PDE, whose solution is generally difficult to be found. On the other hand, it is the non-analytical points of OLF that indicate the exact $$D_A$$. The mentioned points produce a ground for using a rational approximation of the OLF. The coefficients of Taylor series expansion of OLF, $$V$$, about the origin can be computed without solving the Equation (2.7), which characterizes the Theorem 2.1. Lemma 4.1 The coefficients of Taylor series expansion of the OLF introduced in Theorem 2.1 can be determined exactly up to any order. Proof See Appendix A. □ Theorem 2.1 and Lemma 4.1 allow one to consider the OLF in the form of (3.1). On the other hand, PFA algorithm gives an approximation in the form of (3.2). So, due to Theorem 3.1, this $$W(x)$$ itself is a Lyapunov function for system (2.6). $$D_A$$ of system (2.6) may be estimated by considering this $$W(x)$$ as a Lyapunov function; i.e. a region provided by Lemma 2.1. It should be noted however, that the main aim of the proposed approach is to consider the roots of denominator of $$W(x)$$ as an approximation of non-analytical points of $$V(x)$$ and, consequently, to consider the roots of denominator of $$W(x)$$ as an estimation of boundary of $$D_A$$ for (2.6). In other words, a region containing the origin and surrounded by the roots of denominator of $$W(x)$$ is an estimate for $$D_A$$. So far, we have not discussed the order of approximation. Indeed, $$n_a$$ and $$n_b$$ remains as the degrees of freedom for the user. Since estimation of $$D_A$$ is determined by denominator of the approximated function, the degree of denominator must be rich enough. For example, for $$n_b=2$$, denominator can only show elliptic, parabolic and hyperbolic shapes. For more complicated shapes, higher order $$n_b$$ is needed. After $$n_b$$ is set, it is preferred to set $$n_a$$ as the minimum value satisfying the following inequality. (n+nan−1)≥(n+nbn), (4.1) where $$n$$ is the number of states of the system. Here, the method is reviewed in the cast of following pseudo algorithm. Padé-like estimation of domain of attraction algorithm: 1. Consider (2.6) as the system equation. 2. Using Equation (2.7) and Lemma 4.1, calculate Taylor expansion coefficients of OLF up to a sufficient order and consider OLF as (3.1). 3. Consider (3.2) as an approximation of OLF. Set $$n_b$$ and $$n_a$$. The relation (4.1) and related argument may be considered as a suggestion. 4. Using the PFA algorithm, determine unknown coefficients $$a_i$$s and $$b_i$$s. 5. Consider the region containing the origin and surrounded by roots of denominator of (3.2) as an estimation of $$D_A$$ for system (2.6). Remark 4.1 As stated in Remark 3.3, if Taylor series coefficients of OLF are available, then the rest of the algorithm can be implemented by numerical programs. Equations (A.12) and (A.15) in the proof of Lemma 4.1 (see Appendix A) show that for some systems, obtaining Taylor series coefficients of OLF can be implemented by numerical programs. Polynomial systems are placed in this class of systems. Actually, (A.12) and (A.15) are linear equations with the elements of $${\partial^i{V(x)} \over \partial{x^j}}{|_{x=0}}$$ as unknown variables. Therefore, if factors of these unknowns ($${\partial^k{f(x)} \over \partial{x^l}}{|_{x=0}}$$s) are available, which is the case for polynomial systems, then $${\partial^i{V(x)} \over \partial{x^j}}{|_{x=0}}$$ can be linearly obtained by a numerical program. Also, this privilege holds for non-polynomial systems for which factors of $${\partial^k{f(x)} \over \partial{x^l}}{|_{x=0}}$$ may be computed numerically, e.g. a system containing terms such as $$\sin$$, $$\cos$$, etc. However, in general, computing Taylor series coefficients of OLF may need some symbolic computations. First, let us apply this method to a simple scalar system in order to illustrate the method. Consider the following scalar system: x˙=−sin⁡(x)−sin⁡(2x). The origin of the system is locally exponentially stable and the set {x∈R | −23π<x<23π≃2.094} is its domain of attraction (Fig. 1 shows $$\dot{x}$$ vs. $$x$$). Now, suppose that we want to approximate OLF $$V(x)$$ by a rational approximation with $$n_a=n_b=2$$ which means: W(x)=a0+a1x+a2x21+b1x+b2x2. Fig. 1. View largeDownload slide Plot of $$\dot{x}=-\sin(x)-\sin(2 x)$$ shows the domain of attraction. Fig. 1. View largeDownload slide Plot of $$\dot{x}=-\sin(x)-\sin(2 x)$$ shows the domain of attraction. It needs the Taylor series expansion of $$V(x)$$. A simple calculation results in: V(x)=16x2+124x4+O(x6). Noting Remark 3.2, one finds that the Padé-like approximation for this function is as the conventional version of the Padé approximation. Therefore, $$W(x)$$ is obtained as: W(x)=16x21−14x2. Roots of denominator of $$W$$ gives an estimate for the $$D_A$$. The estimation reads as the set: {x∈R | −2<x<2}. Now, the efficiency of the proposed method can be shown through some well-known examples. Example 4.1 Consider van der Pol system in the form of: {x˙1=−x2x˙2=x1−x2+x12x2.  $$D_A$$ of this system is surrounded by an unstable limit cycle. This system is a common example for comparison of the existing methods (see Tan & Packard, 2008, Balint et al., 2009). In Padé-like estimation algorithm, some values are first suggested for $$n_a$$ and $$n_b$$, following of which, other possible values may be tried in the line of enhancing the estimation. The results are shown in Fig. 2. For comparison, Fig. 3 simultaneously shows the exact $$D_A$$, results of the proposed method and of some other previous ones. Order of approximation in Fig. 3 for the proposed method is $$n_a=28, \ n_b=6$$ and the the corresponding denominator of the Padé-like estimation can be found in the Appendix C. Fig. 2. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. This simulation shows the improvement of estimation when order of approximation increases. Fig. 2. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. This simulation shows the improvement of estimation when order of approximation increases. Fig. 3. View largeDownload slide The Padé-like estimation result is close to the boundary of exact $$D_A$$ of Example 4.2, which is an unstable limit cycle. Fig. 3. View largeDownload slide The Padé-like estimation result is close to the boundary of exact $$D_A$$ of Example 4.2, which is an unstable limit cycle. Example 4.2 As another system which has a bounded $$D_A$$, the following system can be considered: {x˙1=x2x˙2=−2x1−x2−x13+x1x24+x25.  Different values of $$n_a$$ and $$n_b$$ are tried. Figure 4 shows the result of different approximations based on different $$n_a$$s and $$n_b$$s. Similar to the previous one, this simulation intuitively shows the improvement of the estimation when order of approximation increases. For approximation with $$n_a=28, n_b=6$$, the expression of denominator of approximation is given in Appendix C. Fig. 4. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. Fig. 4. View largeDownload slide Estimated $$D_A$$ for Example 4.1 with different values for $$n_a$$ and $$n_b$$. Example 4.3 When facing with unbounded $$D_A$$, most estimations methods experience a conservative behaviour. One of the advantages of the proposed method is potential of providing the precise estimation of $$D_A$$, even for unbounded $$D_A$$s. As an example, consider a dynamical system in the form of: {x˙1=x2x˙2=−x1−x2+x13.  The proposed method can be tried for different values of $$n_a$$ and $$n_b$$. The quality of improvement of approximation is shown in Fig. 5. Now, the corresponding result can be compared to those of previous ones (See Fig. 6). It is observed that, the proposed method gives the best estimation in comparison with Lyapunov-based estimation methods (see Fermin et al., 2009) and even some of numerical methods such as trajectory reversing method. Expression of denominator for $$n_a=26,n_b=6$$ can be found in the Appendix C. Fig. 5. View largeDownload slide Estimated $$D_A$$ for Example 4.3 with different values of $$n_a$$ and $$n_b$$. Fig. 5. View largeDownload slide Estimated $$D_A$$ for Example 4.3 with different values of $$n_a$$ and $$n_b$$. Fig. 6. View largeDownload slide The Padé-like estimation technique gives the best result for unbounded $$D_A$$ among different methods. For this simulation, $$n_a=26$$ and $$n_b=6$$ (Example 4.3). Fig. 6. View largeDownload slide The Padé-like estimation technique gives the best result for unbounded $$D_A$$ among different methods. For this simulation, $$n_a=26$$ and $$n_b=6$$ (Example 4.3). Example 4.4 This method can give exact OLF when the OLF itself is rational. Consider a system in the form of: {x˙1=−x1+x12x2x˙2=−x2.  The OLF for this example is obtained as: V(x)=2x22+2x12−x1x232−x1x2. (4.2) For all approximations with $$n_a\geq 4$$ and $$n_b\geq 2$$, and considering Remark 3.1, Equation (4.2) is accurately obtained as a result of utilizing the PFA. Example 4.5 The proposed method can be applied to systems having more than two states. Consider the following system: {x˙1=−x1(1−0.5x12−2x22−3x32)x˙2=−x2(1−0.5x12−2x22−3x32)x˙3=−x3(1−0.5x12−2x22−3x32).  Domain of attraction of this system is the set $$\{(x_1, x_2, x_3)\in \mathbb{R}^3 \ | \ 1-0.5 x_1^2-2x_2^2-3x_3^2>0\}$$. Results of the application of the proposed method to this example results are depicted in Fig. 7. Fig. 7. View largeDownload slide Exact and estimated $$D_A$$ for Example 4.5 with different values of $$n_a$$ and $$n_b$$. Fig. 7. View largeDownload slide Exact and estimated $$D_A$$ for Example 4.5 with different values of $$n_a$$ and $$n_b$$. 5. Stability properties Up to this section, a rough estimation of $$D_A$$ has been obtained, which is sufficient in applications where only approximated borders of $$D_A$$ is needed. One of the major benefits of the proposed method is its ability to present unbounded domains. Also, compared to previous techniques, the geometry and shape of the estimated region in the proposed approach is more similar to those of the exact $$D_A$$. On the other hand, several applications need a guarantee for the estimation region to be inside the stability domain, which is not supported by Theorem 3.1. In this section, through Lemma 2.1, such a guarantee is provided. {However, it is worth noting that for high orders of $$n_a$$ and $$n_b$$, the result of estimation is very close to the exact domain. In the followings, $$\bar{D}_{W,c}$$ means the simple connected region that contains the origin and all of its points satisfy the condition $$W(x)<c$$. The term $$D_{W,c}$$ represents the union of all possible $$\bar{D}_{W,c}$$s. Based on Lemma 2.1, if we find a constant $$c$$ such that the $$D_{W,c}$$ where $$W$$ is a Lyapunov function, falls inside the region where $$\dot{W}<0$$, then $$D_{W,c}$$ is an estimation of $$D_A$$. Obviously, the larger value of $$c$$ gives the larger estimation. Suitable value of $$c$$ may be found via solving a non-linear optimization problem. Theorem 3.1 shows that the function $$W$$ in that theorem is a Lyapunov function. %Therefore, $$W$$ may be used via Lemma 2.1 for a guaranteed estimation of $$D_A$$. Consider the following constrained optimization problem: co= minxW(x)s.t.  W˙(x)=0,x≠0. (5.1) It is obvious that $$\dot{W}$$ is negative in the region $$D_{W,c_o}$$ (excluding the origin), otherwise there exists $$\hat{x}\neq 0$$ such that $$\dot{W}(\hat{x})=0$$ and $$W(\hat{x})=\hat{c}_o<c_o$$, which is in contradiction with the fact that $$c_o$$ is the result of optimization problem. Therefore, the region $$D_{W,c_o}$$ is an estimation of $$D_A$$, which is inside the stability region. The optimization problem (5.1) is a non-linear and (usually) non-convex problem. Since for Lemma 2.1, the global solution of (5.1) is needed, the optimal value of $$c_o$$ is found through a numerical algorithm. Now, let us illustrate the optimization-based method via a simple example. The approximation $$W$$ for Example 4.1 with the order of approximation as $$n_a=n_b=2$$ is as follows: W=32x12−x1x2+x221−14x12+314x1x2−114x22. The optimization problem (5.1) has the global result $$c_o=2.4650786$$ which occurs at points $$(x_1,x_2) = (-0.685248, 0.758646)$$ and $$(x_1,x_2) =(0.758646,-0.685248)$$. See Fig. 8 for details. Fig. 8. View largeDownload slide The region $$D_{W,c}$$ which is inside the region $$\dot{W}=0$$. Fig. 8. View largeDownload slide The region $$D_{W,c}$$ which is inside the region $$\dot{W}=0$$. Figure 9 shows the result of this optimization-based method for Example 4.1. All of these estimations are necessarily inside the stability region. This figure shows that higher order of approximation gives a more precise result.. Fig. 9. View largeDownload slide Result of optimization-based method applied to Example 4.1. The legend shows the orders of numerator and denominator of approximation $$W$$. Fig. 9. View largeDownload slide Result of optimization-based method applied to Example 4.1. The legend shows the orders of numerator and denominator of approximation $$W$$. 6. Conclusion It was found that the estimation of $$D_A$$ based on the proposed, i.e. PFA technique gives more precise results mainly due to the fact that rational polynomials have the potential of revealing the possible non-analytic nature of the functions. Compared to other relevant techniques, the superiority of the proposed approach becomes clear when dealing with unbounded $$D_A$$s. Coefficients of Taylor series expansion of OLF may be computed by a recursive algorithm as stated in the proof of Lemma 4.1. It is worth noting that all computations, including those for coefficients of Taylor series expansion, and those for function approximation, are done in a linear framework. In the meantime, since implementation of the proposed technique needs no symbolic software, it is straightforwardly applicable to high-order systems. Unlike the posed examples, there is no need for system equations, which are supposed to be analytic, to be polynomial. In this paper, no analytic proof for the convergence of roots of denominator of the approximation to the non-analytical points of the main function was given; a subject which is in the line of our future research. Acknowledgements The authors would like to thank members of ACSL (Advanced Control Systems Laboratory) of University of Tehran for the fruitful discussion on the subject. References Baker G. A. ( 1975 ) Essentials of Padé Approximants . New York : Academic Press . Bakhtiari R. & Yazdanpanah M. J. ( 2005 ) Designing a linear controller for polynomial systems with the largest domain of attraction via LMIs. International Conference on Control and Automation, IEEE, Budapest . 1 , 449 – 453 . Balint S. ( 1985 ) Considerations concerning the maneuvering of some physical systems. Ann. Univ. Timisoara . Cambridge, UK : Cambridge Scientific Publishers , 23 , 8 – 16 . Balint S. Braescu L. & Kaslik E. ( 2009 ) Region of Attraction and its Application in Control Theory . Cambridge Scientific Publishers . Benzaouia A. ( 2005 ) Constrained stabilization: an enlargement technique of positively invariant sets. IMA J. Math. Control Inf. Amsterdam : North Holland Publishing Co., 22 , 109 – 118 . Google Scholar CrossRef Search ADS Bruijn N. ( 1961 ) Asymptotic Methods in Analysis . North Holland Publishing Co. Chesi G. ( 2004 ) On the estimation of the domain of attraction for uncertain polynomial systems via LMIs. Proceedings of the 43rd IEEE Conference on Decision Control (CDC) , IEEE, Nassau, Vol.1. pp. 881 – 886 . Chesi G. ( 2009 ) Estimating the domain of attraction for non-polynomial systems via LMI optimizations. Automatica , 45 , 1536 – 1541 , ISSN 0005-1098. Google Scholar CrossRef Search ADS Chesi G. Garulli A. Tesi A. & Vicino A ( 2005 ) LMI-based computation of optimal quadratic Lyapunov functions for odd polynomial systems. Int. J. Robust Nonlinear Control , 15 , 35 – 49 . Google Scholar CrossRef Search ADS Chiang H. D. Hirsch M. & Wu F. ( 1988 ) Stability regions of nonlinear autonomous dynamical systems. IEEE Trans. Autom. Control , 33 , 16 – 27 . ISSN 0018-9286. Google Scholar CrossRef Search ADS Cuyt A. ( 1983 ) The QD-algorithm and multivariate Padé-approximants. Numerische Mathematik , 42 , 259 – 269 . Google Scholar CrossRef Search ADS Evans L. C. ( 1998 ) Partial Differential Equations . Providence, RI : American Mathematical Society . Fermin W. Guerrero- Castellanos J. F. & Alexandrov V. V. ( 2009 ) A computational method for the determination of attraction regions. International Conference on Electrical Engineering, Computing Science and Automatic Control (CCE) . IEEE, Toluca, pp. 1 – 7 . ISBN 978-1-4244-4688-9. Genesio R. ( 1984 ) New techniques for constructing asymptotic stability regions for nonlinear systems. IEEE Trans. Circuits Syst. , 31 , 574 – 581 . Google Scholar CrossRef Search ADS Genesio R. Tartaglia M. & Vicino A. ( 1985 ) On the estimation of asymptotic stability regions: state of the art and new proposals. IEEE Trans. Autom. Control , 30 , 747 – 755 . ISSN 0018-9286. Google Scholar CrossRef Search ADS Graham L. R. Knuth E. D. & Patashnik O. ( 1989 ) Concrete Mathematics, A Foundation for Computer Science . Menlo Park, CA : Addison-Wesley Publishind Company . Guillaume P. ( 1997 ) Nested multivariate Padé approximants. J. Comput. Appl. Math. , 82 , 149 – 158 . ISSN 0377-0427 . Google Scholar CrossRef Search ADS Guillaume P. ( 1998 ) Convergence of the nested multivariate Padé approximants. J. Approx. Theory , 94 , 455 – 466 , ISSN 0021-9045 . Google Scholar CrossRef Search ADS Guillaume P. & Huard A. ( 2000 ) Multivariate Padé approximation. J. Comput. Appl. Math. , 121 , 197 – 219 . ISSN 0377-0427. Google Scholar CrossRef Search ADS Guillaume P. Huard A. & Robin V. ( 1998 ) Generalized multivariate Padé approximants. J. Approx. Theory , 95 , 203 – 214 . ISSN 0021-9045. Google Scholar CrossRef Search ADS Hachicho O. ( 2007 ) A novel LMI-based optimization algorithm for the guaranteed estimation of the domain of attraction using rational Lyapunov functions . J. Franklin Inst. , 344 , 535 – 552 . ISSN 00160032 . Google Scholar CrossRef Search ADS Hashemzadeh F. & Yazdanpanah M. J. ( 2006 ) Semi-global enlargement of domain of attraction for a class of affine nonlinear systems . Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications , IEEE, Munich , Vol. 1 , pp. 2257 – 2262 . Jones R. H. ( 1976 ) General rational approximants in N-variables . J. Approx. Theory , 233 , 201 – 233 . Google Scholar CrossRef Search ADS Kaslik E. Balint A. & Balint S. ( 2003 ) Gradual approximation of the domain of attraction by gradual extension of the “embryo” of the transformed optimal Lyapunov function . Nonlinear Stud. , 10 , 67 – 78 . Khalil H. K. ( 2002 ) Nonlinear Systems , 3rd edn . Englewood Cliffs, NJ: Prentice Hall . Kida S. ( 1990 ) Padé-type and Padé approximants in several variables . Appl. Numer. Math. , 6 , 371 – 391 . ISSN 0168-9274 . Google Scholar CrossRef Search ADS Krantz S. G. & Parks H. R. ( 2002 ) A Primer of Real Analytic Functions . Boston: Birkhauser . ISBN 0817642641 . Google Scholar CrossRef Search ADS Melchor-Aguilar D. & Silviu-Iulian N. ( 2007 ) Estimates of the attraction region for a class of nonlinear time-delay systems . IMA J. Math. Control Inform. , 24 , 523 – 550 . Google Scholar CrossRef Search ADS Steele J. M. ( 2004 ) The Cauchy-Schwarz Master Class . New York: Cambridge University Press . Google Scholar CrossRef Search ADS Tan W. & Packard A. ( 2008 ) Stability region analysis using polynomial and composite polynomial Lyapunov functions and sum-of-squares programming . IEEE Trans. Autom. Control , 53 , 565 – 571 . Google Scholar CrossRef Search ADS Tibken B. ( 2000 ) Estimation of the domain of attraction for polynomial systems via LMIs . Proceedings of the Conference on Decision Control. IEEE, Sydney, NSW , Vol. 4 , pp. 3860 – 3864 . Topcu U. & Packard A. ( 2007 ) Stability region analysis for uncertain nonlinear systems . Proceedings of the IEEE Conference on Decision Control , IEEE, New Orleans, LA , Vol. 1 , pp. 1693 – 1698 . Topcu U. Packard A. Seiler P. & Balas G. ( 2010 ) Robust Region of Attraction Estimation . IEEE Trans. Autom. Control , 55 , 137 – 142 . Google Scholar CrossRef Search ADS Valmórbida G. ( 2009 ) Region of attraction estimates for polynomial systems . Proceedings of the IEEE Conference on Decision Control , IEEE, Shanghai , Vol. 1 , pp. 5947 – 5952 . Vannelli A. & Vidyasagar M. ( 1985 ) Maximal Lyapunov functions and domains of attraction for autonomous nonlinear systems . Automatica , 21 , 69 – 80 . Google Scholar CrossRef Search ADS Vidayasagar M. ( 1993 ) Nonlinear Systems Analysis , 2nd edn . Englewood Cliffs, NJ: Prentice Hall . Yazdanpanah M. J. Khorasani K. & Patel R. V. ( 1998 ) Uncertainty compensation for a flexible-link manipulator using nonlinear $$H_\infty$$ control . Int. J. Control , 69 , 753 – 771 . Google Scholar CrossRef Search ADS Yazdanpanah M. J. Khorasani K. & Patel R. V. ( 1999 ) On the estimate of the domain of validity of non-linear $$H_\infty$$ control . Int. J. Control , 72 , 1097 – 1105 . Google Scholar CrossRef Search ADS Appendix A. Proof of Lemma 4.1 Before proving the lemma, statement of some technical lemmas is needed. Now, the previous lemma is extended for multivariable functions. Now, Lemma 4.1 can be proved. According to Proposition 2.1 for any $$k$$, one has: $$c_{k}=\frac{1}{k!} \left.{\partial^{|k|}{V} \over {\partial{x}^k}}\right|_{x=0}$$. Therefore, it is sufficient to show that $$\left.{\partial^{|k|}{V} \over {\partial{x}^k}}\right|_{x=0}$$ may be computed for any $$k\in (\mathbb{N}\cup\{0\})^n$$. An induction can be established on $$|k|$$ as follows. Start with $$|k|=1$$. Using Equation (2.4), there are $${n\choose n-1}=n$$ derivatives of order $$|k|=1$$. So, for $$l=1,2,\dots,n$$, two sides of Equation (2.7) are driven as: ∂∂xl(⟨∇V,f⟩)=∂∂xl(−‖x‖2). (A.8) Equation (A.8) can be expanded as follows: ∂∂xl(∑j=1n∂V(x)∂xj.fj(x))=−2∑j=1nxj. (A.9) Thereupon: (∑j=1n∂2V(x)∂xjxl.fj(x)+∂V(x)∂xj.∂∂xlfj(x))=−2∑j=1nxj. (A.10) Evaluation of Equation (A.10) at origin and noting that $$f(0)=0$$ results in: ∑j=1n( ∂V(x)∂xj|x=0. ∂fj(x)∂xl|x=0)=0. (A.11) Equation (A.11) is an equation in which the derivatives of $$V$$ are unknowns. Considering (A.11) for $$l=1,2,\dots,n$$, a system with $$n$$ linear equation in $$n$$ unknowns is obtained. They can be rewritten in the following matrix form: [∂f1(x)∂x1∂f2(x)∂x1…∂fn(x)∂x1∂f1(x)∂x2∂f2(x)∂x2…∂fn(x)∂x2⋮⋮⋱⋮∂f1(x)∂xn∂f2(x)∂xn…∂fn(x)∂xn]x=0[∂V(x)∂x1∂V(x)∂x2⋮∂V(x)∂xn]x=0=[00⋮0]. (A.12) Coefficients matrix in (A.12) is equal to $${\partial{f} \over \partial{x}}^T|_{x=0}$$. Since, based on the assumption of Theorem 2.1, the matrix $${\partial{f} \over \partial{x}}|_{x=0}$$ is Hurwitz; therefore, it is invertible. Invertibility of coefficients matrix in (A.12) means that this equation has a unique solution and this solution is:$${\partial{V(x)} \over \partial{x_l}}|_{x=0}=0$$ for $$l=1,2,\dots,n$$. This result was expectable due to positive definiteness of $$V(x)$$. Assume that for all $$|k|=\mathfrak{n-1}$$, $${\partial^{|k|} \over \partial{x^k}}V(x)|_{x=0}$$ can be computed. Set $$|k|=\mathfrak{n}$$. For all $$|k|=\mathfrak{n}$$, the equations which contain derivatives of order $$\mathfrak{n}$$ of $$V(x)$$ as unknowns are obtainable through: ∂|k|∂xk(⟨∇V,f⟩)=∂|k|∂xk(−‖x‖2). (A.13) The left-hand side of (A.13) can also be expanded as follows: ∂|k|∂xk(∑j=1n∂V(x)∂xj.fj(x)) =∑j=1n∂|k|∂xk(∂V(x)∂xj.fj(x)) =∑j=1n(∑i=0k(ki)(∂|k|−|i|∂xk−i∂V(x)∂xj).(∂|i|∂xifj(x))) =∂|k|∂xk(−‖x‖2). (A.14) The right-hand side of the last equality in (A.14) for $$\mathfrak{n}>2$$ is equal to zero. Also, for $$\mathfrak{n}=2$$, when one of the elements of $$k$$ is $$2$$ such as $$k=(0,0,\dots,2,\dots,0)$$, then, the right-hand side of the last equality in (A.14) is equal to $$-2$$; otherwise, it is like $$k=(0,\dots,1,0,\dots,1,\dots,0)$$, which is equal to zero. However, due to $$f(0)=0$$, (A.14) can be rewritten as:  ∑j=1n(∑|i|=1(ki)(∂|k|−|i|∂xk−i∂V(x)∂xj).(∂|i|∂xifj(x))) =∗∗∗ (A.15) The $$***$$ on the right-hand side of the above equation represents all derivatives of $$V(x)$$ whose order are lower than $$\mathfrak{n}$$, higher order derivatives of $$f$$ and constant terms. When considering $$|k|=\mathfrak{n}$$, all terms in $$***$$ are known and do not affect solvability of (A.15). These equations are linearly independent and consequently are solvable. Actually, the matrix of coefficients of above equations is a block diagonal matrix whose blocks are leading principle minors of $${\partial{f} \over \partial{x}}|_{x=0}$$. Therefore, it is invertible and consequently the mentioned equations are solvable. Appendix B. Proof of (3.8) First, consider the following lemma from Steele (2004). In order to prove (3.8), it is sufficient to show that: (|xk|)2≤(‖x‖|k|)2 or equivalently: (|xk|)2|k|≤(‖x‖)2. (B.2) Inequality (B.2) may be expanded as follows: (x12)k1|k|(x22)k2|k|…(xn2)kn|k|≤(x12+x22+⋯+xn2). (B.3) The left-hand side of (B.3) is in the form of (B.1) with $$x_i^2$$ as $$a_i$$ and $${k_i \over |k|}$$ as $$p_i$$ and $${k_1 \over |k|}+{k_2 \over |k|}+\dots+{k_n \over |k|}={|k| \over |k|}=1$$. Therefore, (B.1) states that: (x12)k1|k|(x22)k2|k|…(xn2)kn|k| ≤k1|k|.x12+k2|k|.x22+⋯+kn|k|.xn2 ≤x12+x22+⋯+xnn and the last inequality holds because $${k_i \over |k|}\leq 1$$. Appendix C. Some simulation expressions Expression of denominator of approximation for Example 4.1 with $$n_a=28,n_b=6$$: $$1.0+ 0.1065412730\,x_{{1}}x_{{2}}+ 0.1362438857\,{x_{{1}}}^{3}x_{{2}} - 0.2242705150\,{x_{{1}}}^{2}{x_{{2}}}^{2}+ 0.1013585171\,x_{{1}}{x_{{2}}}^{3}- 0.04151182324\,{x_{{1}}}^{5}x_{{2}}+ 0.04835446107\,{x_{{1}}}^{4}{x_{{2}}}^{2}- 0.02548402153\,{x_{{1}}}^{3}{x_{{2}}}^{3}+0.009070713464\,{x_{{1}}}^{2}{x_{{2}}}^{4}\\- 0.002330734027\,x_{{1}}{x _{{2}}}^{5}- 0.1480416345\,{x_{{1}}}^{2}- 0.1098080759\,{x_{{2}}}^{2}-0.01336682711\,{x_{{1}}}^{4}\\- 0.01896447641\,{x_{{2}}}^{4}-0.003112455890\,{x_{{1}}}^{6}+ 0.0002849304531\,{x_{{2}}}^{6}$$ Expression of denominator of approximation for Example 4.1 with $$n_a=28,n_b=6$$: $$1.0- 0.2790014167\,x_{{1}}x_{{2}}+ 0.2368419702\,{x_{{1}}}^{3}x_{{2}}- 0.7081714653\,{x_{{1}}}^{2}{x_{{2}}}^{2}+ 0.3648060899\,x_{{1}}{x_{{ 2}}}^{3}+ 0.1404756418\,{x_{{1}}}^{5}x_{{2}}+ 0.05357554804\,{x_{{1}}}^{4}{x_{{2}}}^{2}- 0.05682371974\,{x_{{1}}}^{3}{x_{{2}}}^{3}+0.5660169893\,{x_{{1}}}^{2}{x_{{2}}}^{4}+ 0.1787257727\,x_{{1}}{x_{{2}}}^{5}+ 0.07123538082\,{x_{{2}}}^{6}- 0.8135596832\,{x_{{1}}}^{2}- 0.05029478641\,{x_{{2}}}^{2}\\- 0.5301539075\,{x_{{1}}}^{4}- 0.7668262956\,{x_{{2}}}^{4}+ 0.005283318119\,{x_{{1}}}^{6}$$ Expression of denominator of approximation for Example 4.3 with $$n_a=26, n_b=6$$: $$1.0- 0.002314998116\,{x_1}^{4}+ 0.001060978122\,{x_2}^{4}- 0.002242078600\,{x_1}^{6}- \\ 0.000002985127786\,{x_2}^{6}- 0.2781498140\,{x_1}^{2}- 0.1527167633\,{x_2}^{2}+\\ 0.004763572498\,x_1{x_2}^{3}- 0.001975789673\,{x_1}^{5}x_2- 0.0009491287752\,{x_1}^{4}{x_2}^{2}-\\ 0.0003555944435\,{x_1}^{3}{x_2}^{3}- 0.0001177270691\,{x_1}^{2}{x_2}^{4}- \\ 0.00002591183460\,x_1{x_2}^{5}- 0.01080307989\,{x_1}^{3}x_2+ 0.005802793836\,{x_1}^{2}{x_2}^{2}- 0.3086406808\,x_1x_2$$ © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

### Journal

IMA Journal of Mathematical Control and InformationOxford University Press

Published: Jan 11, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create folders to

Export folders, citations