Construction of Lyapunov functions for the estimation of basins of attraction

TECHNICAL PAPERS

^[e-mail] TU Darmstadt, System Reliability and Machine Acoustics, Dynamics and Vibrations Group, Germany

^[e-mail] TU Berlin, Department of Applied Mechanics, Chair of Mechatronics and Machine Dynamics, Germany

^III [e-mail]

^IV [e-mail]

[e-mail] TU Darmstadt, System Reliability and Machine Acoustics, 5Dynamics and Vibrations Group, Germany

^[e-mail] PUC-Rio, Department of Mechanical Engineering, Brazil

ABSTRACT

Technical systems are often modeled through systems of differential equations in which the parameters and initial conditions are subject to uncertainties. Usually, special solutions of the differential equations like equilibrium positions and periodic orbits are of importance and frequently the corresponding equations are only set up with the intent to describe the behavior in the vicinity of a limit cycle or an equilibrium position. For the validity of the analysis it must therefore be assumed that the initial conditions lie indeed in the basins of attraction of the corresponding attractors. In order to estimate basins of attraction, Lyapunov functions can be used. However, there are no systematic approaches available for the construction of Lyapunov functions with the goal to achieve a good approximation of the basin of attraction. The present paper suggests a method for defining appropriate Lyapunov functions using insight from center manifold theory. With this approach, not only variations in the initial conditions, but also in the parameters can be studied. The results are used to calculate the likelihood for the system to reach a certain attractor assuming different random distributions for the initial conditions.

Keywords: Keywords: Lyapunov functions, basins of attraction, center manifold theory

Introduction

In a number of technical applications the behavior of structures is strongly nonlinear and the stationary motions depend on the initial conditions and on a number of parameters. Well known examples of such systems are the snap-through of arches and shells, squealing states in brake systems, self-excited vibrations in paper machinery and many more. In order to accurately model and predict the behavior of a structure, it is important to determine possible stationary solutions and their basins of attraction. A common approach to estimate basins of attraction is the use of Lyapunov functions. For autonomous systems, it is well known how to construct Lyapunov functions for the linearized system in order to prove stability of the solution of the nonlinear problem. However, these Lyapunov functions often only yield crude and technically insufficient estimations of the basins of attraction. Therefore, the goal of this paper is to develop Lyapunov functions from which the basins of attraction of solutions can be estimated more accurately. Construction methods for Lyapunov functions for stability investigations have been developed, for example, by Aizermann and Schultz-Gibson (see Unbehauen, 2000) or by Vannelli and Vidyasagar (1985) or Giesl (2007, 2007, 2012). However, it turns out that the procedures are sometimes inconvenient in practice.

In this paper, we use insight gained from center manifold theory to construct Lyapunov functions, making use of the fact that the stability behavior of a system is often determined on a low dimensional manifold. Since the behavior of a nonlinear structure strongly depends on its initial conditions and on its parameters, which are not accurately known, uncertainties have to be taken into account. Therefore it is also studied here how the basins of attraction change due to small changes in the parameters. If initial conditions are randomly distributed the results can be used to calculate the likelihood for reaching an attractor.

Estimation of Basins of Attraction Through Lyapunov Functions

The task of calculating basins of attraction arises for time continuous and time discrete systems. It is well known that basins of attraction can be studied through Lyapunov functions. In this paper we concentrate on the time continuous case. The underlying theorem can be formulated as:

Theorem: Let V (x) be a Lyapunov function for the time continuous system

with f(0) = 0. The domain

where c is a positive real constant, which belongs to the basin of attraction G of x = 0 (see La Salle, 1967).

The proof is based on the fact that the conditions on S certify that any solution vector starting in S strictly monotonically approaches the origin, as can also be easily visualized by geometric intuition.

It is important to note that the theorem gives only a one-sided estimate and the basin of attraction can be much larger than S. The theorem does not say anything on how to choose the function V. A common, straightforward choice is to construct a quadratic Lyapunov function for the linearized system and then to apply it to the nonlinear system and check for the basin of attraction. This procedure is also the basis of Lyapunov's indirect method in which the stability of the origin can be studied through the linearized system. However, the standard quadratic Lyapunov functions often yield poor approximations for the domain of attraction.

As an example we use the well known Van der Pol equation¹:

For the parameter the corresponding phase diagram is shown in Fig. 1. Using V(q₁,q₂) = with q₁ = q and q₂ = q, the corresponding estimate for the basin of attraction is shown as the shaded area in the figure. Clearly it is a poor approximation of the exact domain of attraction, which corresponds to the area within the unstable limit cycle. The reason for this is of course that a quadratic Lyapunov function can only yield an ellipsoid in the phase space, whereas the basin of attraction can have a very different shape.

In order to obtain a more accurate estimate of the basin of attraction we need to construct a more suitable Lyapunov function. In the above theorem we observe that in the case of a one-dimensional system the basin of attraction coincides exactly with the estimate S. Therefore, the goal is to reduce the number of the states of the problem as much as possible. It is well known that in bifurcation problems of sufficiently smooth systems the long term dynamics of the system are dominated by a low dimensional center manifold (Troger and Steindl, 1991). There are different possibilities to calculate the governing equations for the center manifold. One possibility is to use a polynomial ansatz for the expressions of the center manifold. After substitution into the governing equations the coefficients can be calculated by comparison. In this paper we use a different approach by making use of normal form reduction as suggested by Hochlenert (2012). As usual, any time autonomous dynamical system, as the Van der Pol equation for example, can be written as

where

is a diagonal matrix of the eigenvalues of the system linearized about the trivial solution, which corresponds to the stability boundary. The matrices F_i are coefficient matrices so that terms of the form F_ixⁱ denote polynomials of order i in x. In this representation the symbols xⁱ represent vectors of all monomials of the order i

and the F_i contain the corresponding coefficients. We now try to simplify equations (4) by a near identity transformation

where again denotes polynomials of order i in y. The coefficients of g(y) are chosen such that the resulting equations in normal form

are as simple as possible, meaning that they contain as few terms as possible. In order to calculate the coefficients of g(y), we insert (7) and (8) in (4) and obtain

which is a partial differential equation. Using the expansions for g(y) and h(y)one can sort for the orders of yⁱ and thus obtain

The first of the equations is trivially fulfilled, since linear terms are not changed by the near identity transformation. A comparison of the coefficients for the equations of second order yields

where j is the index for the row and the two indices K and I cover all quadratic terms. Details can be found in Hochlenert (2012). In order to simplify (8) as much as possible we try to choose G_j,kI o that H_j,kI vanishes. This is possible unless the term in the curly bracket vanishes, which can be written as

and is the so called resonance condition. The resulting coefficients are

If the resonance condition is met we cannot eliminate the term and choose

Once H₂ and G₂ are determined one can continue the same process to higher orders. For expressions of order N one obtains

where the tilde is used because have been obtained from the manipulation of lower order terms and are not the original coefficients from Eq. (4). By the same reasoning as above we derive the general resonance condition

If the resonance condition is not met we obtain

and otherwise

At the bifurcation point where the system (3) loses stability, a pair of complex conjugate eigenvalues has a zero real part and all other eigenvalues have a negative real part. From the resonance condition (20), as found in Hochlenert (2012), one can see that the critical equations corresponding to rows with critical eigenvalues, i.e. those with zero real parts, decouple from the rest of the equations of the normal form to arbitrary order. Therefore these equations represent the behavior of the system on the center manifold. This dimension reduction of the system due to the normal form transformation implies that the following procedure for the construction of a Liapunov function can be used for systems with an arbitrary number of degrees of freedom undergoing a codimension one bifurcation with a purely imaginary pair of eigenvalues. The Van der Pol equation (3) studied here should be considered as a prototype system possibly embedded in a system with many degrees of freedom. The corresponding normal form (decoupled from the stable subsystem) in either case reads

where only terms up to order four in are stated. Introducing the polar coordinates , the normal form can be written as

Considering the normal form in polar coordinates, we can define a Lyapunov function

which reproduces exactly the basin of attraction of the trivial solution up to the order to which the normal form has been calculated. In agreement with the well known inverse function theorem the transformation from the physical coordinates q_i to the coordinates y_iof the normal form can be inverted in a neighborhood of the equilibrium position q_i = 0.

For this the key is the inversion of . According to the inverse function theorem we can expand

The coefficients K_i can again be calculated by comparison of coefficients. Substitution of (31) in (7) yields

and hence

The physical coordinates are then recovered by inverting the diagonalization of the linear system matrix.

Using this procedure, the Lyapunov function (29) and the stability boundary can be transformed back to physical coordinates. If we now use in the above theorem, we obtain the basin of attraction shown in Fig. 2 by the shaded area. Of course, it is a conservative estimate, meaning the basin of attraction cannot be overestimated. The approximation of the stability boundary obtained from inverting the stability boundary in the normal form equations almost coincides with the exact one. So far the method proposed is valid for differential equations of which only the initial conditions are of interest. If we are to study variations and uncertainties of the parameters, this can be easily done by considering them as states. An equation of the type

where q is the state vector and p are the parameters can be written as

where is the augmented state vector. This state augmentation is sometimes referred to as "suspension trick''. The parameter now occurs in the equations of the center manifold. For a variation of the parameter from 0.05 to 0.6, which can for example be caused by an uncertainty, the basins of attraction are shown in Fig. 3. At this point we note that the stability of the augmented system is defined only with respect to the physical coordinates and not with respect to the augmented state variables.

The procedure described above can be performed for time discrete systems in an analogous manner (Spelsberg-Korspeter et al., 2011). In the following we describe how the proposed method can be used to calculate the likelihood of systems to reach an attractor under random initial conditions.

Systems with Stochastic Initial Conditions

In this section we use the method proposed for the estimation of basins of attraction in order to calculate the likelihood of the solution to reach a certain attractor assuming different random distributions for the initial conditions. For simplicity, we use on the one hand a normal distribution, which is often used as a first approximation to describe real-valued random variables and on the other hand a uniform distribution, where the initial conditions are equally probable within a specified region. Nevertheless, in the following calculation, any other distribution could also be used.

The normal distribution of the initial conditions of the Van der Pol equation is described by the probability density function (pdf) p_n of the form

The parameters and are the standard deviations of the initial condition of the state variables q₁₀ and q₂₀. Here, the initial conditions are normally distributed around the origin, so that the mean of both initial conditions is zero. Figure 4 shows the normal distribution pnin the q10-q20 space. Of course, in real systems the initial conditions are bounded and cannot reach infinite values. A distribution, which fulfills this requirement is the uniform one. It can be described by the probability density function pu, which we defined as

The region, where the probability density of the initial conditions is constant and different from zero, is assumed to be circular and is characterized by the radius R_u. Figure 5 shows the probability density function p_u in the q₁₀-q₂₀ space.

The random distributions of the initial conditions described above can be used to calculate the likelihood for reaching an attractor. For any distribution described by the probability density function p the probability P of a solution to start in the domain of attraction G of the trivial solution is defined as

where the integration is performed over the basin of attraction G. If instead of the exact domain of attraction G the integration in (40) is performed over the estimated domain S G, then an estimate for this probability is obtained. In Figs. 6 and 7 the exact and the estimated basin of attraction G and S are plotted together with the contour lines of the probability density function. Since the calculated estimate for the basin of attraction is conservative, it will always be smaller than the exact basin of attraction. This also means that the probability, calculated by using the estimate of the basin of attraction, is conservative and therefore a lower bound for the exact probability ().

Since the estimate of the basin of attraction depends on the parameter in equation (3), which is not accurately known, any small changes in this parameter will cause changes in the probability. Therefore, we calculate the exact probability for each possible parameter using the exact basin of attraction as integration boundary as well as the estimated probability using the estimated basin of attraction.

In Fig. 8 the probability P is plotted over the parameter in the range from 0.05 to 0.6 for a normal distribution of the initial conditions with a standard deviation of , described by (38). Due to a better approximation of the basin of attraction for smaller values of , a better approximation of the probability is observed. The error between the exact and the estimated probability is shown in Fig. 9. The maximum error for the investigated range of the parameter is about 4.1%.

For a uniform distribution of the initial conditions described by (39) with R_u = 1.95, the exact and estimated probability is plotted in Fig. 10. attraction for uniform distribution with The estimate for the probability is quite poor since regions exist, which are not included in the estimated basin of attraction, see Fig. 7. The probability error reaches its maximum of 13.3% for = 0.6, shown in Fig. 11. If we use a different uniform distribution with R_u = 1.9, the estimated probability approximates the exact probability much better for smaller values of the parameter , see Fig. 12. For larger values of the parameter the probability error still reaches 11.4%, which is shown in Fig. 13. The good approximation for smaller values of is due to the fact that in this parameter region all initial conditions lie inside both the approximate and the exact basins of attraction, in the exact as well as in the estimated one. When the parameter reaches a certain limit the probability error starts to increase slightly. Here, it seems to be necessary to increase the order of the approximation of the basin of attraction to achieve better results for the probability.

For certain distributions of the initial conditions it is possible that for any parameter variation all possible initial conditions lie inside both basins of attraction. This is possible for example, if we use a uniform distribution with R_u = 1.3. Here, the estimated probability corresponds to the exact probability, see Fig. 14. Due to this fact the probability error is always equal to zero.

Conclusion

In this paper we propose a method to construct Lyapunov functions for the estimation of basins of attraction. Whenever uncertainties have to be taken into account, a purely linear stability analysis is not sufficient but has to be checked for its range of applicability. Since parameters can be interpreted as states of an augmented system, these questions can be answered by the calculation of basins of attraction. This is often done through the study of Lyapunov functions, for which, however, systematic construction approaches are rare. In this paper we used center manifold and normal form theory in order to construct improved Lyapunov functions to calculate estimates of the basins of attraction. Using normal form theory, every system can be transformed to a reduced one in the vicinity of a hopf bifurcation. Additionally, system parameters can be taken into account by adding trivial differential equations, or for "dynamical parameters" non-trivial differential equations. Thus, the analysis of uncertainties in parameters is reduced to an initial value problem and an estimation of basins of attraction. The construction method should also be valuable for control problems.

Using this construction method we also performed a stochastic analysis to analyze the effect of randomly distributed initial conditions on the reachability of attractors. The quality of the results depended on the calculated estimate for the basin of attraction, which can be improved by increasing the approximation order in the construction method. For further investigations it is possible to include not only uncertainties in the initial conditions, but also in system parameters.

Paper received 11 July 2012.

Paper accepted 18 August 2012.

• Giesl, P., 2012, "Construction of a finite-time Lyapunov function by Meshless Collocation'', Discrete and Continuous Dynamical Systems Series B, 17 (7), pp. 2387-2412.
• Giesl, P., 2007, "On the determination of the basin of attraction of discrete dynamical systems'', Journal of Difference Equations and Applications, 13 (6), pp. 523-546.
• Giesl, P., 2007, "Construction of Global Lyapunov Functions Using Radial Basis Functions'', Lecture Notes in Mathematics 1904, Springer, Heidelberg.
• Hochlenert, D., 2010, "Dimension reduction of nonlinear dynamical systems: Comparison of the center manifold theory and the method of multiple scales'', 81st annual meeting of the GAMM, Karlsruhe, Germany, 22nd to 26th March 2010.
• Hochlenert, D., 2012, "Normalformen und Einzugsbereiche nichtlinearer dynamischer Systeme'' (Normalforms and basins of attraction of nonlinear dynamical systems), Habilitationsschrift, TU Berlin.
• La Salle, J.L., Lefschetz, S., 1967, "Die Stabilitätstheorie von Ljapunow'', Bibliographisches Institut Mannheim.
• Spelsberg-Korspeter, G. Hochlenert, D. and Hagedorn, P., 2011, "Non-linear investigation of an asymmetric disk brake model'', Journal of Mechanical Engineering Science, Proceedings of the Institution of Mechanical Engineers Part C, DOI: 10.1177/0954406211408531.
• Troger, H. and Steindl, H., 1991, "Nonlinear Stability and Bifurcation Theory'', Springer Verlag Wien.
• Unbehauen, H., 2000, "Regelungstechnik (Control Theory)'', Vieweg Verlag.
• Vannelli, A., Vidyasagar, M., 1985, "Maximal Lyapunov Functions and Domains of Attraction for Autonomous Nonlinear Systems'', Automatica J. IFAC 21, 1, pp. 69-80.

Construction of Lyapunov functions for the estimation of basins of attraction

G. Spelsberg-Korspeter^I; D. Hochlenert^II; E. Heffel^III; A. Wagner^IV; P. Hagedorn^V; R. Sampaio^VI

1

The Van der Pol equation with

is a well known paradigm for self-excitation. Changing the sign of

corresponds to inverting time and the stability behavior of the attractors. The attractors themselves are not affected.

Publication Dates

• Publication in this collection
  24 July 2013
• Date of issue
  2012

History

• Received
  11 July 2012
• Accepted
  18 Aug 2012