Exact Region of Attraction Determination Based on Manifold Method

The exact region of attraction plays an important role in autonomous nonlinear system, while the results based on the conventional method, such as Lyapunov function approach, are always conservative. However, results via the manifold method, which is the main approach studied, are exact. This method optimizes the distribution of points on the circle through modifying the end point of the former trajectory and inserting/deleting point on the circle on the basis of trajectory arc length method to improve the accuracy and efficiency. First, the basic theory of manifold method is introduced. Secondly, a methodology for determining stable manifold are proposed, which is the core of the manifold method in stability boundary determining. Finally, on this basis, three examples about academic model, power system and aviation system are taken to illustrate the advantages of the method. The results show that the method can improve the accuracy and significantly reduce the calculation time, and can be widely used in engineering systems.


I. INTRODUCTION
The region of attraction (ROA), also named stability region, which plays an important role in electrical power system [1]- [3], aviation [4], [5], economics, ecology, etc., is a domain around an asymptotically stable equilibrium point, in which trajectories starting will converge to the equilibrium point. Of cause, there are lots of methods for ROA determination, such as Lyapunov function approach [6], reachable set [3], [7], Monte Carlo [4] and so on. But all of them have valuable advantages as well as inevitable shortcomings. For instance, Lyapunov function approach as a well-known and general method which can obtain the explicit function for the boundary of the ROA, and it can easily distinguish whether a point is inside the stability boundary or not. Whereas, the Lyapunov function is usually too complex to determine by computer and the conservatism of the ROA also can't be acceptable. Theoretically, Monte Carlo approach as a numerical method utilizes a number of points in state space to estimate the stability region. However, the efficiency of this method is quite low.
The associate editor coordinating the review of this manuscript and approving it for publication was Meng Huang .
The exact ROA determination is of critical importance in improving the system's performance, however, the results based on the methods above are always conservative. Another non-Lyapunov and numerical method for finding the ROA based on its topological properties to some extent solves the problem. In 1988 and 1989, Hsiao-Dong Chiang [8], etc. proposed that the stability boundary of dynamic system consists of the union of the stable manifolds of all equilibrium points on the stability boundary and the method presented can obtain the exact ROA. And it can be also seen in our recent study [5] about an aircraft to determine the dynamic envelope based on the method. Nowadays, with the development of the control theory, many advanced control theories are wildly used in nonlinear systems, For instance, fuzzy adaptive finite-time faulttolerant control for strict-feedback nonlinear systems [9], event-triggered robust fuzzy adaptive finite-time control of nonlinear systems with prescribed performance [10], fuzzy output tracking control and filtering for nonlinear discretetime descriptor systems [11],and some theories even been applied to a robotic airship against sensor faults [12] and attitude tracking of hypersonic vehicle [13]. ROA as the stability boundary of the nonlinear system, once the system exceeds the ROA, the state of the system will diverge, which is important to guide the application of nonlinear control theories. However, consulting similar references about ROA determination, we found that the method proposed was not widely used from then on because of the difficulty in computing the stable manifolds for the unstable equilibrium points on the stability boundary. With the development of the computer technology, the determination for two or even more dimensional manifold is not ever a difficult problem, which gives the method an opportunity in exact ROA determination and that is the reason why the method can be proposed once again.
At present, there are several acceptable methods to compute ROA for a dynamic system, of which three are the most commonly used. The first one is reachability analysis [7], which can be formulated as an optimal control problems, and the ROA can be characterized using variants of the Hamilton Jacobi Bellman (HJB) or Isaacs (HJI) partial differential equations, the accuracy of this method relies on grids of state space, as the grids of the states grow, the computational load will increase exponentially, therefore, the calculation time of this method is very long. The other two methods are trajectory arc length method (it's also called trajectory reversing method) and geodesic circle method which are both based on manifold theory to estimate the ROA [14]. Geodesic circle method determines the grid position by solving the boundary value problem, the distance between adjacent circle is determined by the point of the fastest manifold on the former circle, which ensures the high calculating precision, however, it has to adjust the distance and recalculate all the points on the circle when individual points don't satisfy the angle constraints. Therefore, it needs to be calculated repeatedly, which adds extra computation [15]. Trajectory arc length method first determine the initial manifold near the unstable equilibrium point, and make the manifold grow uniformly in all directions, and adjust the points on the circle after a certain integral time, eventually form the entire manifold. This method has the merits of simple calculation and high efficiency, while the accuracy of the calculation decreases rapidly with the increase of the length of manifold trajectory, and this makes the manifold far from the equilibrium point less accurate [16].
This paper optimizes the distribution of points on the circle through modifying the end point of the former trajectory and inserting/deleting point on the circle on the basis of trajectory arc length method to ensure the accuracy and efficiency. Meanwhile, this paper takes a parallel process to speed up the computation process. The simulation results show the feasibility and accuracy of this method, and the examples shows that this method can be widely used in the engineering nonlinear system.
The organization of the paper is as follows. Section II presents introduction of fundamental concepts of manifold method for determining ROA. And the algorithm for estimating the stability regions of nonlinear dynamic system is also presented in the section. The numerical method for manifold computation, also regarded as calculation for the part of the stability boundary, is proposed in the section III. Several engineering examples are proposed in Section IV. Finally, conclusions are drawn in Section V.

II. BASIC THEORY FOR MANIFOLD
Consider a nonlinear autonomous dynamical system described byẋ (1) where f : R n → R n is continuous and sufficiently smooth, and satisfies Lipschitz condition in x.
Supposex is the hyperbolic equilibrium point of f , satisfying that f x = 0 and the eigenvalues of Jacobian matrix Df has no zero real part.
There are two classes of equilibrium points (EP), stable equilibrium point (SEP) and unstable equilibrium point (UEP). The anterior one is an equilibrium point, at which every eigenvalues of the Jacobian matrix only possess negative real part. And the latter one is an equilibrium point, at which eigenvalues of the Jacobian matrix possess positive real part. If the number of eigenvalues with positive real part is equal to n, the dimension of the system, the UEP is called a source; or, the UEP is called a saddle. In addition, type of the UEP is equal to the number of the positive real part, for example, if a UEP only possess a positive real part at his eigenvalues; the unstable equilibrium point is called type-1 UEP.
For a saddle x u , the stable manifold W s (x u ) and unstable manifold W u (x u ) exist which can be defined as the following, where t (x) is called the flow induced by the vector field f , also the trajectory. Moreover, for a SEP x s , there is a stability region, denoted by A (x s ), that is the set of all point which the trajectory started from converges to x s . And the stability region can be expressed by, And the boundary of the stability region is denoted by ∂A (x s ). Furthermore, it can be sure that an equilibrium point is a UEP on the boundary under the condition that W u And based on the method mentioned above for determining the stability region, the boundary of stability region ∂A (x s ) is the union of the stable manifolds of equilibrium points on the stability boundary; ∂A (x s ) can be expressed as following where x i , i = 1, 2, · · · are the equilibrium points on the stability boundary. And the proof of this theorem can be seen in Ref. 8. The details of how to determine the stability boundary are as follows. VOLUME 8, 2020 Step 1: Find all the equilibrium points via solving the equations f (x) = 0.
Step 2: Identify the stable equilibrium point x s whose eigenvalues of Jacobian matrix only have negative real part.
Step 3: Identify the unstable equilibrium points x ui , i = 1, 2, · · · on the stability boundary via judging whose unstable manifolds contain trajectories approaching the stable equilibrium point x s .
Step 4: Determine the stability boundary posed by all the stable manifolds of every unstable equilibrium points obtained from step 3.
From the steps above, we know that step.1 to step.3 are carried out easily, which will not be discussed in detail, while there is something difficulty in step.4 due to dynamic characteristic of the system, that is the important reason why the method of ROA determination is not attached any importance by most researchers up to now. And the part called manifold computation will be discussed in Section III.

III. NUMERICAL METHOD OF MANIFLOD COMPUTION
As known from the manifold theorem, the global manifold can be evolved to sweep out from a local manifold in it. And flow belonging to the unstable manifold on the boundary of the region of attraction, which is also an invariant manifold, will not depart from it forever. Thus, the global manifold can be evolved from the initial set around the unstable equilibrium point on the boundary. The details of the algorithm utilized in the paper are as following.
Step.1 calculate the eigenvalues ε 1 , ε 2 , ε 3 and its eigenvectorsv 1 , v 2 , v 3 of the Jacobian matrix Jx, (suppose the real part of ε 1 > 0, the real part of ε 2 , ε 3 ≤ 0) Step.2 Determine the initial normal vector η of the manifold. If ε 2 , ε 3 are negative real number, η = v 1 × v 2 ; and if ε 2 and ε 3 are conjugate complex roots, Step.3 Determine the initial circle around the UEP. Considering that local stable manifold around a UEP belongs to the stable subspace posed by the stable eigenvectors of the dynamical system at the UEP, the initial circle as a local stable manifold can be a circle around the UEP with a very small radius in the space posed by eigenvectors, of which eigenvalue possess negative real part.
Step.4 Calculate and modify the end point of the former trajectory to get a smooth circular trajectory of the next generation circle. The next generation circle is posed by the end points integrated in reverse time from the points in the former generation.
(a) Calculate the endpoint correct time t add where l is the limiting length of the trajectory, l real is the real length of the trajectory when the calculation is stopped.
(b) calculate the correction point of the terminal point x new In the paper, the growing length of the trajectories can be various rather than constant to accelerate the computation and strengthen stability of calculation results.
Step.5 Calculate the distance between two adjacent points and recompose points in the circle. Due to the different dynamic characteristic in each dimension, distance of adjacent points in next generation may be too small or large rather than equivalent, even though points in the former generation circle are uniform distribution. Thus, it is necessary to recompose the point in the circle. The method is inserting points in two adjacent points, if the distance is too large, and deleting a point in two adjacent points, if the distance is too small.
Define a parameter , which represents the ideal distance between adjacent points, if 0.5 < d < 1.5 , keep this point; if d > 1.5 , insert round(d/ ) − 1 points; if d < 0.5 , delete this point.
Step.6 Build the global manifold by repeating the step.4 and step.5 until the total increment length is equal to the premeditated value.
In order to verify the accuracy and precision of the algorithm for manifold computation, we take following examples. The first one is Lorenz system [9], [10] studied by most researchers for two-dimensional stable or unstable manifold computation. And the other is a four-dimensional Hamiltonian system studied by Hinke M Osinga [18] for computing the two-dimensional invariant manifolds in four-dimensional dynamical systems. For the Lorenz system, the state equations are as follow.
x = a (y − x) where, a=10, b=8/3, and c=28. In the case, origin is one of the type-one UEP and has a two-dimensional stable manifold. Based on the method above, we visualize the manifold of Lorenz system in Figure 1(a), of which the initial circle of radius R = 2 and the total trajectory air length L = 120. At the same time, to illustrate the merits of this method, it's compared with the geodesic circle method, which has been verified to have high accuracy. The stable manifold of Lorenz system based on geodesic circle method is shown in Figure 1(b).
Where, N represents the number of points on the circle, T IVP represents the time required to calculate the initial value of each point.
As we see in Figure 1 and Table 1, the accuracy of stable manifold calculated by the manifold method is nearly consistent with the geodesic circle method, however, the calculation of manifold method is significantly less than the geodesic circle method. For the four-dimensional Hamiltonian system arising from control theory, the state equations are as follow.
where,    In the case, origin is one of the type-two UEP and has a two-dimensional stable manifold. Based on the method above, we visualize the manifolds of four-dimensional Hamiltonian system in Figure 2 and Figure 3, of which the initial circle of radius R = 1 and the total trajectory air length L = 30. Considering the similar of the projection in various 3-D state space, which are (x 1 , x 2 , p 1 ), (x 1 , x 2 , p 2 ), (x 1 , p 1 , p 2 ), and(x 2 , p 1 , p 2 ), the two kind of projections are illustrated in figure 2 for (x 1 , x 2 , p 1 ) 3-D state space and figure 3 for (x 1 , p 1 , p 2 ) 3-D state space.
From the Ref. 19, we can sure that the method for manifold computation is correct and has high precision, because of the same manifolds of the Lorenz system and Hamiltonian system based on other approaches and this method. And this is the important precondition for an exact ROA determination via manifold method. where, In this example, there are four equilibrium points showed in the Table 2. And the point A is an SEP, while the other points are UEP, whose type can be seen in Table 1. In addition, it can be verified that the point B is on the stability boundary of the point A. In order to verify the accuracy of the ROA based on the manifold method, the ROA of the same dynamic models mentioned above is also determined via Monte Carlo method, which can yield an exact ROA around the SEP at a certain range. And for illustrating clearly, the ROA based on the manifold method is drawn as line and that of the Monte Carlo is plotted as colorful surface. And the results are drawn in Figure 4. And in the figure, results determined based on manifold method and Monte Carlo have the same trend, which can be regarded as the proof for feasible and accurate of the manifold method.

B. ROA FOR POWER SYSTEM
This is an application example about single machine infinite bus system (SMIB) in power system, of which dynamic model can be represented as follow. where, where, ω is the relative speed of the generator; ω b is the reference speed of the generator; δ is the power angle of the generator; X d is d-axis reactance of the generator; X q is q-axis reactance of the generator; X d is d-axis transient reactance of the generator; X t is reactance of the transformer; X l is the reactance of the transmission line; T do is the time constant of the excitation winding; T j is the inertia constant; E q is the q-axis transient voltage; P m is the mechanical input power; P e is the active power delivered by the generator; Id is the d-axis current; I q is the q-axis current; D is the damping coefficient; U is the infinite bus voltage. In addition, ω, δ and E q are the variable state of this dynamic system and the other parameters above are constant, of which the values are as follow.
X q = 0.72 T j = 8.75 There are two equilibrium points showed in the Table 3. And the point A is an SEP, while the point B is a typeone UEP. Figure 5 is the ROA determined via the manifold method.
In addition, the manifold method is suitable for the system with some uncertain parameters. And in order to provide a   good proof, the same SMIB system, of which the infinite bus voltage is under ±10 uncertainty, is taken as an example. Of course, the method is also suitable for the system of other uncertain parameters. The stability boundaries with uncertain parameters are both drawn in the Figure 6. As we can see in the figure, the rule of ROA change is coincident with the effect of the infinite bus voltage change, and the ROA enlarges when the infinite bus voltage increases and the ROA reduces when the infinite bus voltage decreases.

C. ROA FOR NASA'S GENERIC TRANSPORT MODEL
This is an example about a condition, about which there are more than two equilibrium points on the stability boundary and one of these points are not the type-1 UEP. The example is about the NASA's Generic Transport Model (GTM) in aviation, of which longitudinal dynamic model can be expressed as follow. (13) where V t is the flight velocity, α is the angle of attack, q is the pitch rate, θ is the pitch angle, m is the mass,c is the mean aerodynamic chord, ρ is the atmospheric density, g is the acceleration of gravity, S ref is the reference wing surface area, J y is the moment inertias along aircraft Y axis,q = 1 2 ρV 2 t represents the dynamic pressure. The main structure parameters can be obtained by consulting full-scale model of NASA's GTM-T2 [20] and the detailed polynomial model parameters of the coefficients C x , C z , C m can be found in paper which is obtained from flight test [21].
In addition, a pitch hold controller is designed for the longitudinal dynamic model to bring the aircraft state to the reference pitch attitude and improve the flight quality; of cause, maneuvering envelope (stability boundary) of an aircraft with PID controller also can be determined by this method. The control law simulated in this paper is expressed below: where α = α 0 −α is state error; k α , k q and k θ are the control coefficients for the state feedback controller; θ ref is the pitch angle command. k α = −2, k q = 1 and k θ = −1 are utilized in order to make the aircraft have better flight quality during the landing phase. Consider a landing phase, the speed V t = 85 m/s, and the altitude H = 400 m: In this example, we regard V t as a constant. Calculate the equilibrium point of the aircraft during the landing phase, and on this basis, estimate the ROA for the three-order system of α, θ , q.
Except the SEP, there are four unstable equilibrium points showed in the Table 4. And the type-1 UEP C, E and type-2 UEP B, D are the UEPs on the stability boundary of the SEP.     In the figure, the red surface is the stability manifold of the type-1 UEP E; the blue surface is the stability manifold of the type-1 UEP C. Shown in the figure, the boundary of the ROA is the union of the stability manifold of the type-1 UEP C and E. To illustrate the accuracy of this method, the ROA of the aircraft is also determined via Monte Carlo method, as we can find in Figure 8, the ROA obtained by this method is completely conformity with the Monte Carlo method. Table 5 shows the calculation time between different methods, and the method proposed in this paper can significantly reduce the calculation time, and improve the calculation efficiency.

V. CONCLUSION
In the paper, the main content is the ROA determination based on the manifold method, which refers that the stability boundary is the union of the stability manifold of the UEP on the boundary. The advantages of manifold method are as follows, which obtained via studying three examples.
(1) Comparing with Monte Carlo method, the ROA computation based on the manifold method is faster, and comparing other methods is more accurate.
(2) The manifold method can be widely used in the engineering nonlinear system, such as aviation, power system and so on.
(3) The manifold method can also be applied to estimate the ROA of the nonlinear system with uncertain parameters.
However, there exists some problem in this method, for instance, when the dimension of the system is higher, the accuracy of this method will decrease, at the same time, it is difficult to visualize the ROA.
In the future work, efforts will be focused on the further improvement of the algorithm and the application of algorithm in engineering practice. In the part of improved algorithm, we will use adaptive factors to further optimize the distance between the adjacent circles and the position of the points in the circle to improve the accuracy. In engineering practice, we will take the new method to estimate the safe flight envelope of aircraft quickly and accurately, which can be used to guide the pilot's manipulates under abnormal conditions to avoid loss of control. And it can also be used to estimate the maneuvering envelope of hypersonic flight envelope, which is the basis of the application of advanced modern control theory.