A fast sampling method for estimating the domain of attraction

Most stabilizing controllers designed for nonlinear systems are valid only within a speciﬁc region of the state space, called the domain of attraction (DoA). Computation of the DoA is usually costly and time-consuming. This paper proposes a computationally effective sampling approach to estimate the DoAs of nonlinear systems in real time. This method is validated to approximate the DoAs of stable equi-libria in several nonlinear systems. In addition, it is implementedforthepassivity-basedlearningcontroller designed for a second-order dynamical system. Simulation and experimental results show that, in all cases studied, the proposed sampling technique quickly estimates the DoAs, corroborating its suitability for real-time applications.


Introduction
The domain of attraction (DoA) of a stable equilibrium in a nonlinear system is a region of the state space from which each trajectory starts and eventually converges to the equilibrium itself. In the literature, the DoA is also known as the region of attraction or basin of attraction [1,33]. The DoA of an equilibrium and its computation is of main importance in control applications. However, in most cases, computation of the DoA is quite costly. This paper aims to approximate the DoAs of nonlinear systems in real time by introducing a sampling approach.
Several techniques have been proposed in the literature to compute an inner approximation for the DoA [8], which can broadly be classified into Lyapunovbased and non-Lyapunov methods [11]. Lyapunovbased approaches include, for instance, sum of squares (SOS) programming [4], methods that apply both simulation and SOS programming [32], procedures that use theory of moments [13]. In this approach, first, a candidate Lyapunov function is chosen to show asymptotic stability of the system within a small neighborhood of the equilibrium. Next, the largest sublevel set of this Lyapunov function, in which its time derivative is negative definite, is computed as an estimate for the DoA [25]. Non-Lyapunov methods include, for instance, trajectory reversing [11,24], determining reachable sets of the system [2], and occupation measures [14,18].  Although Lyapunov-based techniques have been successfully implemented for estimating the DoAs of various nonlinear systems [8], there are still two main issues with using these approaches. The first is that most of the existing methods are limited to polynomial systems [12,31]. As such, in the case of nonpolynomial systems, first, the equations of motion are approximated by using the Taylor's expansion and then the DoA is computed based on the approximated polynomial equations. The second is that the available methods are usually computationally costly and timeconsuming which make them unsuitable for real-time applications [7]. This paper proposes a fast sampling approach for Lyapunov-based techniques to estimate the DoAs of various nonlinear systems. This method is computationally effective and is beneficial for real-time applications. In this procedure, once a candidate Lyapunov function is chosen, a sampling algorithm searches for the largest sublevel set of the Lyapunov function such that its time derivative is negative definite throughout the obtained sublevel set. The proposed sampling approach is applied to approximate the DoAs of several nonlinear systems, which have been already investigated in the literature, to validate its capability in comparison with the existing methods. In addition, we go beyond these examples and implement it to compute the DoAs of the passivity-based learning controller [28] designed for an inverted pendulum. This paper is organized as follows. Section 2 reviews the process of estimating the DoAs of nonlinear systems using Lyapunov-based techniques. Section 3 describes the sampling approach and provides a comparison between the estimated DoAs computed by the sampling method and by the existing optimizationbased methods. Section 4 illustrates the DoAs approxi-mated for the controller learned for an inverted pendulum. Finally, Sect. 5 concludes the paper after a short discussion on the properties of the sampling algorithm.

Estimating the domain of attraction using Lyapunov-based methods
Consider the dynamical systeṁ where x ∈ X ⊆ R n is the state vector, u ∈ U ⊆ R m is the control input, and f : X × U → R n is the system dynamics. For a specific state-feedback controller Φ(x) the closed-loop system is described bẏ If x * is a stable equilibrium of the closed-loop system and x(t, x 0 ) denotes the solution of (2) at time t with respect to the initial condition, the DoA of controller Φ is defined by the set An analytical method to approximate the DoA is defined via Lyapunov stability theory as follows [6,16]. If the equilibrium is nonzero, without loss of generality, we can replace the variable x by z = x −x * , wherē x * is the nonzero equilibrium. As such, we can study the stability of the associated zero equilibrium [1]. The conditions of Theorem 1 ensure that the approximated set M is certainly contained in the DoA. The choice of a candidate Lyapunov function is not a trivial task and the DoA approximation relies on the shape of the Lyapunov function's level sets. A procedure to find an appropriate Lyapunov function has been proposed in [10], where gradient search algorithms are implemented to compute a candidate Lyapunov function. Moreover, using composite polynomial Lyapunov functions [29] and rational Lyapunov functions instead of quadratic ones might lead to better approximations, since these have a richer representation power (see, e.g., [9,34]). Quadratic Lyapunov functions restrict the estimates to ellipsoids which are quite conservative [30]. A rational Lyapunov function is written in the form where R i (x) and Q i (x) are homogeneous polynomials of degree i, which are constructed by solving an optimization problem [34]. The sublevel set V(c) of the Lyapunov function V (x) is defined by According to Theorem 1, any sublevel set of a candidate Lyapunov function that satisfies the locally asymptotic stability of the equilibrium can be an estimate for the DoA if the time derivative of the Lyapunov function is negative everywhere within the sublevel set. Since the largest sublevel set provides a more accurate estimate, the problem of approximating the DoA is converted to the problem of finding the largest sublevel set of a given Lyapunov function [15]. To attain the largest estimate for the DoA, one needs to find the maximum value c ∈ R for V(c) such that the computed set satisfies the conditions of Theorem 1.
Theorem 2 [8] The invariant set V(c * ), which is a sublevel set of the Lyapunov function V (x), is the largest estimate of the DoA for the origin of system (2) This can be approached as an optimization problem that has been solved by using SOS programming, methods that apply both simulation and SOS programming, and methods that use theory of moments. However, these techniques are typically restricted to systems and Lyapunov functions described by polynomial equations. In this paper, we present an alternative approach using the sampling approach.

Sampling method for estimating the domain of attraction
The sampling approach presented in this paper has the same goal as the Lyapunov-based optimization approaches have: Find the largest sublevel set of a candidate Lyapunov function to approximate the DoA. We explicitly evaluate the conditions stated in Theorem 1 for a given Lyapunov function with respect to a randomly chosen state x i . The level sets associated with the sample x i with positive derivative of the Lyapunov function are discarded. We propose two sampling methods, memoryless and with a memory, designed to achieve tighter estimates.

Memoryless sampling
This method searches for the upper bound of the parameter c * in (6). First, a state x i is randomly chosen within X or its user-defined subset and the conditions of Theorem 1 are checked for V (x i ) andV (x i ). If these conditions are not satisfied, the upper bound of c * , denotedĉ * , is decreased to the valueĉ * = V (x i ) and the sublevel set V(ĉ * ) is computed as an overestimation for the DoA. At the beginning of the algorithm,ĉ * is initialized atĉ * = ∞. As the sampling proceeds for a large number of samples (n s ) throughout the state space, the value ofĉ * converges to c * from above and the obtained largest sublevel set V(ĉ * ) will be very close to V(c * ).
Since this procedure just focuses on the upper bound of c * , the achieved estimates are not tight enough and the condition ofV (x) < 0 may not be satisfied for some regions of the attained sublevel set as the computed valueĉ * is actually larger than the real value c * . This algorithm may exceptionally not exclude very small regions whereV (x) ≥ 0 from the DoA approximated. However, the empirical evidence arising from extensive simulations suggests that, in practice, the proposed algorithm converges to the exact level set for a sufficiently large number of samples. Based on the practical results, we found that this technique is very fast and its result is very close to the reported estimates in the literature for various classes of systems. Moreover, it does not require computer memory to save the results computed since once a new value is computed forĉ * , its current value is replaced by the new value. Algorithm 1 summarizes this method for estimating the DoA of a given stable equilibrium.

Algorithm 1 Memoryless sampling method for estimating the DoA
As an example, consider a pendulum described by the following nonlinear dynamic equations where x 1 is the angle of the pendulum measured from the vertical axis and x 2 is the angular velocity. The state vector is defined by x = [x 1 x 2 ] T . We use the sampling method with a uniform distribution to approximate the DoA of the stable equilibrium x = (0, 0). To compute a candidate Lyapunov function, first the dynamic Eq. (7) are linearized around the equilibrium and then the candidate Lyapunov function is computed in the form V (x) = x T Px, where P is the solution of the Lyapunov equation A T P + P A + Q = 0 with the identity matrix Q. In this example, the candidate Lyapunov function is obtained as Figure 2 illustrates the evolution ofĉ * of the sampling approach with n s = 500 samples. The real value c * for the candidate Lyapunov function (8), calculated by solving the optimization problem (6), is c * = 9.287 and the value computed by our method isĉ * = 9.702.

Sampling with memory
This method updates both the lower and the upper bounds of c * denoted c * andc * , respectively. Together, As the sampling proceeds, after a large number of samples, the value of c * increases, but not necessarily monotonically. Eventually it converges to c * and the largest sublevel set V(c * ) is obtained. Moreover, the value ofc * monotonically decreases and converges to c * from above. When the conditions of Theorem 1 are satisfied for state x i , the value of V (x i ) is stored in an array as a possible estimate for c * . This is required to guarantee that the approximated DoAs computed by the lower bound of c * always verify the conditions of Theorem 1. This leads to tighter estimates. The array, denoted E, contains 0 initially. The length of this array, without counting its initial element, is in the worst case n s − 1.
is a potential estimate for the DoA. In the caseV (x i ) ≥ 0 and V (x i ) <c * , if c * ≥c * then the algorithm looks for a new lower bound c * among the values stored in the array E. The maximum value of c * is chosen from E such that c * < c * . Selecting a previously stored lower bound satisfies the conditionV < 0 for the obtained sublevel set V(c * ).
Although the sampling algorithm with memory is a conservative method, it may exceptionally overesti-mate the DoA, for instance, when the region described byV (x) < 0 is not simply connected. In such a case, the algorithm may not exclude small holes inside the region in whichV (x) ≥ 0. A formal guarantee for convergence of this algorithm does not exist yet, but the empirical result attained from extensive simulations and experiments illustrates that the sampling technique converges to the exact level set for a sufficiently large number of samples, in practice. Algorithm 2 describes the sampling method with memory for estimating the DoA.
if c * ≥c * then 12: c * = arg max{c ∈ E : c <c * } 13: end if 14: end if 15: end for 16: return c * We apply this approach with a uniform distribution sampling to approximate the DoA for the equilibrium of the pendulum example. Figure 3 illustrates the values of the lower and upper bounds of c * throughout the sampling process with 500 samples where c * = 9.174. Figure 4 depicts the approximated DoA of the equilibrium. The black ellipsoid represents the DoA estimate with c * = 9.271, the dashed blue line, which determines the boundary of the light blue area, represents the region in whichV (x) < 0, and the arrows represent the system trajectories. If the trajectories start inside the DoA estimate, they certainly converge to the origin. The randomly chosen sampling states, which are 500 samples in this example, are represented by red points throughout the state space.

Repeatability of the sampling method
To check the repeatability of the proposed sampling approach, we run various instances of the process of estimating the DoA for the equilibrium of the pendulum example. Figure 5 illustrates the mean value of c * andc * (i.e., (c * +c * )/2) and its standard deviation by a black line and green bars, minimum of c * and maximum ofc * by blue dashed lines at each sample in a simulation where the sampling method runs 1000 iterations each

Directed sampling
In the pendulum example, we used a uniform distribution for sampling the state space or its subset. However, if the structure of the level sets of the Lyapunov function are known, other distributions can be used to avoid sampling in areas of the state space which are already known not belong to the DoA. It is desirable to sample inside the largest level set found so far, specially in its boundary.
In general, sampling with an arbitrary distribution is a challenging problem. Two main approaches exist in the literature: rejection sampling and inverse transform sampling [3]. These techniques focus on sampling the relevant locations of the state space at the cost of computational complexity. In situations where evaluating a particular sample is costly (due to a complicated Lyapunov function or system dynamics), the extra cost incurred by sampling from a complex distribution may be negligible.
To test the trade-off between the speed of convergence and the computational cost, we have applied three different sampling approaches to the pendulum example (7). The uniform sampling on a fixed box (a subset of the state space) is compared with uniform sampling mapped through polar coordinates to lie inside the largest found valid level set, and with exponential sampling mapped through polar coordinates to lie around the boundary of the largest found valid level set. Figure 6 illustrates the sampling points selected by the three types of distributions. The obtained data corroborate the hypothesis that different sampling leads to different convergence rates. Figure 7 presents the convergence statistics for 1000 iterations with 500 samples each. The exponential polar sampling converges the fastest and has the lowest variation between c * and c * while converging. This can be explained by observing Fig. 6c that most of the samples are focused around the boundary of the level set. For this particular example, the cost of evaluating the Lyapunov function and its time derivative is low, but the computation time increases with the complexity of the sampling algorithm. Table 1 shows the average computation time of each sampling method with 500 samples, implemented in the Mathematica software on an Intel core i7 2.7 GHz microprocessor.

Sampling method versus optimization-based methods
Both the sampling and optimization-based methods require a candidate Lyapunov function for estimating the DoA. Table 2 represents six dynamical systems with quadratic Lyapunov functions selected from the literature. The dynamic equations of the first three examples are polynomial and the equations of the last three are non-polynomial. Examples E3 and E6 are third-order systems and the others are second-order systems. For each system, the maximum possible value of c * computed by the sampling approach with 1000 samples is compared with the result of optimization-based methods, reported in the literature. The estimates attained by the sampling technique are very close to the estimates derived by optimization-based methods. In some cases, such as example E2, the result of the sampling procedure is even more accurate. The last column of Table 2 presents the simulation time for approximating the DoA of each system using the sampling approach, implemented in the MATLAB R2014a software on an Intel core i7 2.7 GHz microprocessor. Similarly, Table 3 illustrates three dynamical systems with rational Lyapunov functions selected from  the literature. Example E7 is a second-order polynomial system, E8 is a second-order non-polynomial system, and E9 is a third-order polynomial system. Table 4 presents their corresponding rational Lyapunov functions based on (4). The maximum possible value of c * obtained by the sampling approach with 1000 samples is compared with the result of optimization-based methods, reported in the literature. The result of this comparison validates the proposed sampling technique particularly for non-polynomial systems. The simulation time for approximating the DoA of each system using the sampling procedure is given in the last column of Table 3. Figure 8 depicts the approximated DoAs obtained by the sampling method for the origins of examples E1-E9. Based on the obtained results, it is concluded that the proposed sampling approach is suitable for estimating the DoAs of both polynomial and non-polynomial systems. It is computationally effective and computes the DoA estimate considerably fast. Although the sampling 2.655 2.887 8.6 1.320 1.318 16.6 method may offer less accurate estimates for the DoA at times, it is very useful for real-time applications. It is also beneficial for the control schemes applying the controllers' DoAs such as online sequential composition approaches [21][22][23].

Experimental results
Consider the inverted pendulum, shown in Fig. 9, which is modeled by the nonlinear differential equation where q is the angle of the pendulum measured from the upright position, J is the pendulum inertia, m is the mass, l is the pendulum length, and b is the viscous mechanical friction. Moreover, K is the motor constant, R is the electrical motor resistance, and u is the control input in Volts which is saturated at ±3 V. The state vector of the system is defined by x = [q p] T with p = Jq the angular momentum. Table 5 presents the values of the physical parameters of the pendulum. These values have been found partly by measuring and partly estimated using nonlinear system identification.
The algebraic interconnection and damping assignment actor-critic (A-IDA-AC) algorithm, proposed in Table 4 Rational Lyapunov functions for the systems of Table 3 Example Numerator terms of the Lyapunov function (4) Denominator terms of the Lyapunov function (4) E7 [19,26] [20], is implemented to obtain swing-up and stabilization of the pendulum at the desired upper equilibrium x d = (q d , p) = (0, 0). The goal of this algorithm is to find a proper control input after a number of learning trials. Monitoring the DoA of the learned controller at every trial provides a stopping criterion to terminate learning once the DoA is large enough to fulfill the control objective. This leads to learning in a short amount of time. The parameterized control policy is given bŷ where ϑ ∈ R n is a parameter vector, Ψ ∈ R n is a userdefined basis function vector, and γ is a unit conversion factor with the value of one. The parameter vector ϑ is updated using the actor-critic reinforcement learning (RL) method by following the procedure described in [20]. Consequently, the saturated control input of the A-IDA-AC algorithm is computed at each time step by where Δu k is a zero-mean Gaussian noise, as an exploration term.
The desired system Hamiltonian is chosen in the quadratic form We exploit the desired system Hamiltonian as a candidate Lyapunov function to approximate the DoAs of the learned controllers at each learning trial [17]. Figure 10 illustrates the approximated DoAs of the learned controllers computed by the sampling approach at seven specific trials, where the trial numbers are also given. While learning is in progress, the DoAs of the learned controllers typically enlarge centered at the up equilibrium, but not necessarily monotonically. In this example, after 35 trials, the DoA of the controller is large enough to cover the initial state; hence, the learning process can be terminated sooner instead of running for all the scheduled trials. As such, the sampling method speeds up the process of learning controllers.  x 0 x 0 Fig. 10 Approximated DoAs of the learned controllers at seven specific trials for the inverted pendulum, where the trial numbers are also indicated

Conclusions
This paper has proposed a fast sampling approach for estimating the DoAs of nonlinear systems in real time. The approximated DoAs computed by this technique have been compared with the estimates derived by optimization-based methods. It is concluded that the sampling approach is fast and computationally effective in comparison with optimization-based methods and it can be used for real-time applications. Although a formal guarantee for convergence does not exist yet, the empirical evidence arising from extensive simulations suggests that in practice this approach always converges to the exact level set for a sufficiently large number of samples. Moreover, the rate of convergence depends on the distribution function selected for sampling as well as the exploring regions of the state space. Using a more sophisticated distributed function can speed up convergence of the sampling procedure since it can avoid sampling in areas of the state space which are already known not belong to the DoA. As such, there is a trade-off between the speed of convergence and the computational cost imposed by the complexity of the sampling distribution function. In addition, the sampling approach has been applied to approximate the DoAs of a passivity-based learning controller, designed for an inverted pendulum system, at every learning trial. This online approximation can be used as a stopping criterion for the learning process. This allows learning to be terminated as soon as the controller's DoA is sufficiently large to satisfy the control objective. Thus, the proposed sampling method enables learning in a short amount of time.