Hangos, Katalin

Rozgonyi, Szabolcs

Szederkényi, Gábor

SZTAKI Publication Repository

Improved Estimation Method of Region ofStability for Nonlinear Autonomous SystemsSzabolcs Rozgonyi, Katalin M. Hangos, Gábor SzederkényiResearch Report SCL-002/2006Improved Estimation Method of Region of Stabilityfor Nonlinear Autonomous SystemsSzabolcs Rozgonyi1,2 Katalin Hangos3E-mail: rozgonyi.szabolcs@chello.hu E-mail: hangos@scl.sztaki.huGábor Szederkényi3E-mail: szeder@sztaki.hu1Department of Computer Science, Pannon University, H8201 Veszprém, HungaryH-8201 Veszprém, P.O. Box 158, Hungary2Continental Teves Hungary Ltd., H-8201 Veszprém, HungaryH-8200 Veszprém, Házgyári u. 6-8., Hungary3Process Control Research Group, Systems and Control Research Laboratory,Computer and Automation Institute HASH-1518 Budapest, P.O. Box 63, HungaryAbstractThe stability and the stability region of a (controlled) industrial system is an important propertyto be determined. It is highly desirable to have a system that is not able to become unstable thatis why the stability theory is being continuously in both the theoretical researchers’s and industrialpartitioner’s focus. The aim of this paper is to examine the stability region of an existing industrialplant of great importance using stability region estimation methods can be found in literature (see[1],[2],[3],[4]).Research Report SCL-002/2006 1Contents1 Introduction 22 Basic notions 23 Maximal Lyapunov functions 24 Computation of the Domain of Attraction (DOA) 45 Properties of the recursive algorithm 66 Case studies 76.1 Van der Pole-equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76.2 Lotka-Volterra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86.3 Fermentation system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Conclusion and Future work 8List of Figures1 Pairs of Qi and Ri which should be multiplied in Equation 5.8 . . . . . . . . . . . . . . 102 Domain of attraction of the Van der Pole system . . . . . . . . . . . . . . . . . . . . . 113 Direct check of stability of the Van der Pole-equation . . . . . . . . . . . . . . . . . . . 114 Domain of attraction of the Lotka-Volterra system . . . . . . . . . . . . . . . . . . . . 125 Direct check of stability of the Lotka-Volterra system . . . . . . . . . . . . . . . . . . . 126 Stability region of the fermentation system . . . . . . . . . . . . . . . . . . . . . . . . . 13Research Report SCL-002/2006 21 IntroductionIn this article an improved algorithm will be shown to estimate the region of attraction of nonlinearautonomous systems. The concept is based on [3]. The algorithm will also be analyzed and a numberof examples will be shown.2 Basic notionsIn the following we will consider the nonlinear autonomous system supposing that the origin is itsasymptotically stable equilibrium point:x˙(t) = f(x(t)) (2.1)By the region or domain of attraction S(M) of the set M (which need not to be an attractor) wemean the union of all trajectories with the property that their limit sets are non-empty and containedby M itself. Based on this composition we can define the domain of attraction of the origin as a sethaving only one element.Definition 2.1. The domain of attraction of the origin is the setS(0) = {x0 : x(t, x0)→ 0 as t→∞} (2.2)where x(t, x0) denotes the solution of the system in Equation 2.1 corresponding to the initial conditionx(0) = x0.In this paper it is shown by following the ideas of [3] that there exists a sequence of special kindof Lyapunov functions Vm that can be used to estimate the set S(0) through estimating a Lyapunovfunction of special kind. An iterative method will be given to find these appropriate Vms. The givenalgorithm is able to find unbounded domains of attraction, too. Usually, the first few number ofiterations can show if the domain is bounded or not.Throughout the paper it is assumed that the function f is smooth enough that Equation 2.1 hasunique solution corresponding to each initial condition x(0) = x0. If the solution is Lyapunov-stablethen it depends continuously on x0. Being stable at t0 concludes stability at any other t0′, too.According to a known theorem in [4] if the function f does not depend on t or periodic in t then thestability of the origin implies uniform stability and asymptotical stability implies uniform asymptoticalstability. For the definitions of different types of attraction and stability see for example [5].3 Maximal Lyapunov functionsThe proofs of the statements in this section can be found in [3].Theorem 3.1. Suppose we can find a set A ⊆ Rn containing the origin in its interior, a continuousfunction V : A→ R+ and a positive definite function φ such that1. V (0) = 0, V (x) > 0∀x ∈ A \ {0}2. The functionV˙ (x0) = limt→0+V (x(t, x0))− V (x0)t(3.1)is well defined at all x ∈ A and satisfiesV˙ (x) = −φ(x),∀x ∈ A. (3.2)Research Report SCL-002/2006 33. V (x)→∞ as x→ ∂A and/or ‖x‖ → ∞.Then A = S.Corollary 3.2. Suppose we can find a set A ⊆ Rn containing the origin in its interior, a continuouslydifferentiable function V : A→ R+ and a positive definite function φ such that1. V (0) = 0, V (x) > 0∀x ∈ A\{0}2. ∇V (x)′f(x) = −φ(x)∀x ∈ A3. V (x)→∞ as x→ ∂A and/or ‖x‖ → ∞.Then A = S.Definition 3.3. We say that a function φ : R+ → R+ is of class K if φ is continuous, φ(0) = 0 andφ is monotonically increasing.The next results show that the conditions on V imposed in Theorem 3.1 are reasonable.Theorem 3.4. Suppose f is continuously differentiable in some neighborhood of 0. Then there existsa continuous function Vm : S → R+ and γ ∈ K such that1. Vm(0) = 0, Vm(x) > 0∀x ∈ S\{0}2. ˙Vm(x) = −γ(‖x‖)∀x ∈ S3. Vm(x)→∞ as x→ ∂S. Moreover, if f is Lipschitz continuous on S then Vm can be selected tobe continuously differentiable on S and4. Vm(x)→∞ as ‖x‖ → ∞.Theorem 3.5. Suppose f is Lipschitz-continuous on S. Then in order for an open set A containingthe origin to be the domain of attraction for Equation 2.1, it is necessary and sufficient that there exista continuous function V : A→ R+ and a positive definite function φ such that conditions of Theorem3.1 hold true.Suppose V is a continuous function on some ball Bδ such that V (0) = 0 and V˙ is negative definite.Then one could prove that V is positive definite. This fact shows that if we can find a function V anda positive definite function φ such that V (0) = 0 and∇V (x)′f(x) = −φ(x) (3.3)then V is guaranteed to be positive definite.Definition 3.6. A function Vm : Rn → R+∪{∞} is called maximal Lyapunov function for the systemdescribed in Equation 2.1 if1. Vm(0) = 0, Vm(x) > 0, x ∈ S\{0}2. Vm(x) <∞⇔ x ∈ S3. Vm(x)→∞ as x→ ∂S and/or ‖x‖ → ∞4. ˙Vm is well-defined and negative definite over S.Remark 3.7. Theorem 3.4 shows that a maximal Lyapunov function exist under not too strict condi-tions on f .Research Report SCL-002/2006 44 Computation of the Domain of Attraction (DOA)We look for a function V and a positive definite function φ satisfying V (0) = 0 andV˙ (x) = −φ(x) (4.1)over some neighborhood of 0. Then the boundary of the domain of attraction is defined by the limitV (x)→∞.A systematic procedure will be discussed here to solve Equation 4.1 supposing that the Taylorseries expansion exists for f around the origin. Express f asf(x) =∞∑i=1Fi(x) (4.2)where the functions Fi, i ≥ 1 are homogeneous functions of degree i. For i = 1 we haveF1(x) = Ax,A ∈ Rn×n. (4.3)For the sake of brevity let Fi(x) = 0, i ≤ 0.Our candidate Lyapunov function should excess any limit as x gets closer to the boundary of setS or as ‖x‖ → ∞. For this reason we put D(n) to the denominator, i.e.V (x) =N(x)D(x)(4.4)where N(x) and D(x) are polynomials in x. Thus, V (x) → ∞ as x → ∂S and this suggests thatx ∈ ∂S when D(x) = 0. According to the results of Section 3 the boundary of S is defined by solvingD(x) = 0 for x. We obtain a recursive technique to find this boundary by definingV (x) =∞∑i=2Ri(x)1 +∞∑i=1Qi(x)(4.5)where Ri and Qi are homogeneous functions of degree i. The most straightforward idea for φ is x′Qxwhere Q > 0. Substituting this expression into Equation 4.1 we obtainV˙ (x) = ∇V (x)′f(x) = −φ(x) = −x′Qx (4.6)Based on Equation 4.5 and Equation 4.6 we get((1 +∞∑i=1Qi)∞∑i=2∇Ri′ −(∞∑i=1∇Qi′)∞∑i=2Ri)∞∑i=1Fi = −x′Qx(1 +∞∑i=1Qi)2. (4.7)From this finally we obtain∞∑i=2∞∑k=1∇Ri′Fk+∞∑i=1∞∑j=2∞∑k=1Qi∇Rj′Fk−∞∑i=1∞∑j=2∞∑k=1Qi′RjFk = −x′Qx1 + 2 ∞∑i=1Qi +∞∑i=1∞∑j=1QiQj .(4.8)Equating the coefficients of the same degrees of the two sides of Equation 4.8 we get for degree 2 that∇R2′F1 = −x′Qx (4.9)Research Report SCL-002/2006 5and the general solution when degree k is greater than or equal to 3 isk∑i=2∇Ri′Fk+1−i +k−2∑i=1k−1∑j=2(Qi∇Rj′ −∇Qi′Rj)Fk+1−i−j = −x′Qx(2Qk−2 +k−3∑i=1QiQk−2−i). (4.10)Thus in each step of the algorithm we get the following linear under determined set of equationsas the equivalent form of the previous two equations:Any = bn (4.11)where An are matrices of appropriate dimension. Consider the nonlinear system of equationsx˙ = f(x) =∞∑i=iFi(x). (4.12)First, select homogeneous functions Rn and Qn−2,n ≥ 3, such that the coefficients of Rn and Qn solvethe constrained minimization problem yielded by Equation 4.11min en(y)s.t.An(y) = bn(4.13)where en(y) is the square of 2-norm of the coefficients of degree greater than or equal to n + 1 inthe expression of V˙n. Furthermore, according to the theorem of La-Salle about invariant sets one canchoose the largest positive value C∗ such that the level setVn =n∑i=2Ri1 +n−2∑i=1Qi= C∗ (4.14)is contained in the region given byΩ ={x : V˙n(x) ≤ 0}. (4.15)Then the setSA = {x : Vn(x) < C∗} (4.16)is contained in the region of attraction S. If this is the case then the iteration should be stopped assoon as the desired accuracy has been reached.If en(y∗) = 0 for some y∗ then the iteration can be stopped andV˙n = −x′Qx (4.17)where Q > 0. In this case the domain of attraction is defined by the equationD(x) = 0. (4.18)This means that the domain of attraction S is given by the formulaS ={x :n−2∑i=1Qi > −1}. (4.19)Would any of the afore-mentioned cases happen for the major part of the systems the iteration canbe stopped less then 10 cycles.Research Report SCL-002/2006 65 Properties of the recursive algorithmLet us show that the sequence of functions Vn are Lyapunov-functions. The result in this section isthe corrected form of the one found in the original publication [3].Theorem 5.1. Consider the system of nonlinear equations in Equation 4.12. Assume that the lin-earized systemx˙ = F1(x) = Ax (5.1)is asymptotically stable. Let the homogeneous functions Ri and Qi satisfy the following recursiveequations for second order terms∇R2′F1 = −x′Qx (5.2)and for terms of degree greater than 2∇R2′Fk−1 +k∑j=3(∇Rj′ +j−2∑i=1(Qi∇Rj−i′ −∇Qi′Rj−i))Fk−j+1 =− x′Qx(2Qk−2 +k−3∑i=1QiQk−2−i), k ≥ 3 (5.3)where Q is a fixed positive definite matrix. ThenVn =R2(x) +R3(x) + . . . +Rn(x)1 +Q1(x) + . . . +Qn−2(x)(5.4)is a Lyapunov-function for all n ≥ 2.Proof. The proof is slightly different from the one found in [3]. LetU =(1 +n−2∑i=1Qi(x))2. (5.5)ThenV˙n = ∇Vn′f =1U(1 + n−2∑i=1Qi)n∑j=2∇Rj′ −(n−2∑i=1∇Qi′) n∑j=2Rj f=1U n∑i=2∇Ri +n−2∑i=1n∑j=2(Qi∇Rj′ −∇Qi′Rj) ∞∑k=1Fk =1U(V˙nmain + V˙nresid)(5.6)whereV˙nmain = ∇R2′F1 +n∑k=3∇R2′Fk−1 + k∑j=3(∇Rj′ +j−2∑i=1(Qi∇Rj−i′ −∇Qi′Rj−i))Fk−j+1 (5.7)andV˙nresid =∞∑i=1 n−1∑j=n−i+1∇Rj+1 +n∑r=2n−2∑q=(n+1)−(i−1)−r(Qq∇Rr′ −∇Qq′Rr)Fi −∇R1′ ∞∑i=n+1Fi. (5.8)Research Report SCL-002/2006 7For simpler notation let Qi, Ri and Fi be zeros if i ≤ 0. To understand the structure of V˙nresid youcan consider that the index combinations in the double-sum have characteristic structure as you cansee it on Figure 1 where n = 5. Note that pairs with index i > n should not be taken as summands.Moreover an important property of V˙nresid is that it contains terms only of degree greater than or equalto n+1. Substituting the left hand side of Equation 5.2 and Equation 5.3 into Equation 5.7 we obtainV˙n =1U(−x′Qx− x′Qx(n2∑k=3Qk−2 +n∑k=3k−3∑i=1QiQk−2−i))+1UV˙nresid=1U(−x′Qx(1 + 2n∑k=3Qk−2 +n∑k=3k−3∑i=1QiQk−2−i))+1UV˙nresid (5.9)Because the numerator of the first term is subexpression of U with degree up to and including n − 2Equation 5.9 can be written in the form ofV˙n =1U(−x′Qx(U − Ug) + V˙nresid)= −x′Qx+x′QxUg + V˙nresidU(5.10)where Ug is the sum of terms in U of degree greater than or equal to n − 1. This shows that V˙n isnegative definite around the origin.6 Case studiesIn this section three examples will be shown. In the first and second examples the minimizationproblem in Equation 4.13 cannot be solved such that e(n) becomes small enough that we could applyEquation 4.19. Instead, Equation 4.16 is used for the estimation for the domain of attraction. In thelast example the minimization problem can be solved and one can use Equation 4.19 to get a properregion.In the pictures showing the level sets (Figure 2 and Figure 4) the red colour curve shows V = C∗while the blue one shows those points where V˙ = 0.6.1 Van der Pole-equationThe Van der Pole system we took as example is described with the equation-systemz˙1 = −z2z˙2 = z1 − z2 + z12z2The origin is an asymptotically stable equilibrium point of the system thus the method based onmaximal Lyapunov functions can be used in this case. Applying the iteration steps 9 times we getC∗ = 6.6 and the the minimum value of e(n) was found to be 0.0168945. The stability region is theinnermost area bounded by the red colored curve in Figure 2. We stopped the iteration at step 9because further steps could not increase the minimum significantly.To verify the given region we used a direct method by scanning the points over the region [−3,−3]×[−3,−3] and examining if the system remains stable or not. The resulted set can be seen in Figure 3.We can ascertain that the regions found by the two different methods are pretty the same.Research Report SCL-002/2006 86.2 Lotka-VolterraFor this example we took a Lotka-Volterra system described by the equationsz˙1 = −z1((z1 − 1) (z1 − 3) +12z2)z˙2 = z2(−2.1 + z1)Both the origin and point (2.1, 1.98) both are equilibria of the system. By shifting the secondequilibrium to the origin we get the following centralized systemx˙1 = −0.42x1 − 2.3x12 − x13 − 1.05x2 −12x1x2x˙2 = 1.98x1 + x1x2Just like the previous case after 6 steps of the algorithm we get C∗ = 1.65 and the minimum ofe(n) was 2915.040375. See the region in Figure 4 bounded by the inner red curve.By scanning the points over the region [−1, 1.5] × [−1.5, 1.5] we find that the region estimated bythe Vanelli-method is a subset of the one we found by direct search, see Figure 5.6.3 Fermentation systemIn this example we show a simple model of a fermentation system [6] described by the equationsz˙1 = −0.802228z1 +z1z20.03 + z2 + 0.5z22z˙2 = 0.802228 (10− z2)−z1z20.03 + z2 + 00.5z22which has one asymptotically stable equilibrium point: (4.89067, 0.218662). Shifting the system bythis point we get the transformed equation systemx˙1 = −0.802228 (4.49067 + x1) +2.13881 + 0.437324x1 + 9.79134x2 + 2x1x20.545137 + 2.43732x2 + x22x˙2 = 7.84686 − 0.802228x2 +−4.27762 − 0.874648x1 − 19.5627x2 − 4x1x20.545137 + 2.43732x2 + x22After 3 steps of the algorithm the minimalization problem can be solved that e(n) becomes zerothat is why we can apply Equation 4.19 and we get the region of stability as it is seen in Figure 6. Wecan see that the given region is very small and thus meets our expectations as it is described in [6].7 Conclusion and Future workIn this report an improved algorithm based on constructing maximal Lyapunov functions [3]is proposedto estimate the region of attraction of nonlinear autonomous systems. The advantage of this algorithmis that one does not have to know the solution of the system starting form different initial values; only aminimization problem (a linear programming problem) needs to be solved in each step of the recursiveapproximation procedure. Moreover, the applicable system class is wider than that of the majority ofavailable algorithms can handle (they are mainly restricted to polynomial systems). And lastly, thisalgorithm is faster than the well-known Zubov’s method as only a few steps are needed instead of a fewtimes ten. In addition, it is accurate enough that one could use it instead of the cumbersome (howevermore exact) scanning approach.Research Report SCL-002/2006 9A few problems arose during the analysis and implementation of the described algorithm, they canbe divided into two main area. One of them is the deeper investigation of the method to find broaderclass of systems which can be analyzed, for example, the class of periodic non-autonomous systems orthe so-called partial stable systems (for definitions see [5]) are promising. An other important problemis to seek for restrictive conditions under which the number of iterations could be decreased or at leastestimated.References[1] Chesi, G., Domain of attraction: estimates for non-polynomial systems via LMIs. in Proc. of16th IFAC World Congress on Automatic Control, 2005.[2] Camilli, F., & Grüne, L. (2000). A regularization of Zubov’s equation forrobust domains of attraction. Universität Bremen, Bremen. 00-06, Camillihttp://citeseer.ist.psu.edu/camilli00regularization.html[3] Vanelli, A., & Vidyasagar, M. (1985). Maximal Lyapunov functions and do-mains of attraction for autonomous nonlinear systems. Automatica, 21, 69–80.http://citeseer.ist.psu.edu/context/507042/0[4] Yoshizawa, T. (1966). Stability Theory by Lyapunov’s Second Method. The Mathematical Societyof Japan, 9, 223. http://citeseer.ist.psu.edu/context/389473/0[5] Rouche, N., Habets, P., & Laloy, M. (1977). Stability theory by Liapunov’s direct method.. NewYork: Springer-Verlag. http://citeseer.ist.psu.edu/context/320841/0[6] Szederkényi, G., Kristensen, N. R., Hangos, K. M., & Jorgensen, S. B. (2002). Nonlinear analysisand control of a continuous fermentation process. Computers and Chemical Engineering, 26, 659–670.Research Report SCL-002/2006 10R5 Q5R4 Q4R3 Q3R2 Q2R1 Q1i=3R5 Q5R4 Q4R3 Q3R2 Q2R1 Q1i=4R5 Q5R4 Q4R3 Q3R2 Q2R1 Q1i=1R5 Q5R4 Q4R3 Q3R2 Q2R1 Q1i=2Figure 1: Pairs of Qi and Ri which should be multiplied in Equation 5.8Research Report SCL-002/2006 11-3 -2 -1 1 2 3-3-2-1123Figure 2: Domain of attraction of the Van der Pole system-3. -1.8 -0.6 0.6 1.8 3.-3.-1.8-0.60.61.83.Figure 3: Direct check of stability of the Van der Pole-equationResearch Report SCL-002/2006 12-1 -0.5 0.5 1-1-0.50.51Figure 4: Domain of attraction of the Lotka-Volterra system-1 -0.5 0. 0.5 1.-1.5-0.9-0.30.30.91.5-1 -0.5 0. 0.5 1.-1.5-0.9-0.30.30.91.5Figure 5: Direct check of stability of the Lotka-Volterra systemResearch Report SCL-002/2006 13-0.06-0.03 0.03 0.06-0.09-0.06-0.030.030.060.09Figure 6: Stability region of the fermentation system

Improved estimation method of region of stability for nonlinear autonomous systems

Abstract

Similar works

Full text

Available Versions

SZTAKI Publication Repository