Chaotic Particle Swarm Optimisation for Enlarging the Domain of Attraction of Polynomial Nonlinear Systems

: A novel technique for estimating the asymptotic stability region of nonlinear autonomous polynomial systems is established. The key idea consists of examining the optimal Lyapunov function (LF) level set that is fully included in a region satisfying the negative deﬁniteness of its time derivative. The minor bound of the biggest achievable region, denoted as Largest Estimation Domain of Attraction (LEDA), can be calculated through a Generalised Eigenvalue Problem (GEVP) as a quasi-convex Linear Inequality Matrix (LMI) optimising approach. An iterative procedure is developed to attain the optimal volume or attraction region. Furthermore, a Chaotic Particular Swarm Optimisation (CPSO) e ﬃ cient technique is suggested to compute the LF coe ﬃ cients. The implementation of the established scheme was performed using the Matlab software environment. The synthesised methodology is evaluated throughout several benchmark examples and assessed with other results of peer technique in the literature.


Introduction
It is challenging to ascertain the domain of attraction (DA) for a point in equilibrium in a non-linear dynamical system. DA represents those initial conditions where the system state converges to equilibrium asymptotically [1,2]. Incorporating DA in the design specifications is crucial while working on control synthesis to ensure system stability [3,4]. Specifically, it would be helpful to gauge the optimal Lyapunov function that maximises the DA [5][6][7]. Consequently, it is vital to determine the shape of this region every time one has to scrutinise the system's stability [7]. For this reason, we use the fundamental theory regarding Lyapunov stability (see [1,8,9]). As a matter of fact, for a specific Lyapunov, the biggest estimated area of asymptotic stability may be characterised as the largest set of the function, which is in the area where the derivate is negative [10,11].
It is widely understood that it is non-trivial to study the DA. Considering the nonlinear system, several methods have been suggested for calculating inner estimates, see e.g., [12,13]. For instance, some classical approaches like the Zubov method [14], La Salle method [15] and the trajectory reversing method are studied in [16]. Typically, such techniques use Lyapunov functions and yield inner estimates for the domain of attraction by providing a set of sublevels for a Lyapunov function. These estimates are confirmed mathematically by validating those sets as invariant.
For nonlinear polynomial systems, several techniques have been suggested to handle this step. A majority of these may be used only with this category of system, e.g., Sum of Square (SOS) relaxation-based techniques, which cause a Linear Matrix Inequality (LMI) [17], Bilinear Matrix Inequality (BMI) techniques [18], SOS programming [19], and methods employing simulations and those based on Chebychev points [20]. The interest in enlarging the stability and robustness domain for nonlinear optimal control systems with parametric and structural uncertainties is attractive and well justified for practical process implementation [21,22]. This paper further assesses the potential benefits of using the stability theory for approximating DA using a simpler shape for estimation. A Lyapunov function is used to describe the shape, and it is typically a quadratic function. For a specific Lyapunov function, the calculation of an optimal DA estimate (which is the biggest estimate for the selected shape) is equivalent to solving a non-convex distance problem.
Considering this context, a problem of significant importance is the choice of quadratic LF (QLF). The selected LF for estimating the DA significantly affects the volume of the optimal estimate. It is obvious that the optimal quadratic LF (OQLF) is that which designates the function that provides for maximum volume. Discouragingly, the calculation of the OQLF is akin to solving a non-convex optimisation [23]. The Largest Estimate of DA (LEDA) must be calculated for a given QLF after which the OQLF may be determined.
Our study uses the enlargement synthesis technique substantiated in [24] and extends it to use an evolutionary approach in order incorporate an estimator of a candidate Lyapunov function [6]. The problem of optimising the DA is converted to a generalised eigenvalue problem (GEVP). It is assumed that the Lf describing the form of the largest estimation of the attraction region is a polynomial. Using identified relaxations that have been established around the polynomial sum of square, it is revealed that a lower bound of the largest attainable attraction region can be computed through a generalised eigenvalue problem. This later is considered to be a quasi-convex LMI optimisation where the solution can be computed efficiently.
In this paper, an original strategy aiming to enlarge domains of attraction of stable equilibriums based on maximal Lyapunov functions is developed. The key idea consists in searching the coefficients of the best level set of a Lyapunov function. This latter is supposed to be fully enclosed in the region of negative definiteness of its time derivative. A heuristic optimisation methodology is established, which involves a tangency constraint between the level sets and requirements on the sign of the Lyapunov function. Such constraints help in evading a large number of potential local solutions of the nonlinear optimisation routine. In a first step, a CPSO [25] algorithm is implemented to estimate the maximal Lyapunov function coefficients. Subsequently, an iterative technique based on GEVP is set up for the resultant LF, which provides the maximal estimate for the DA.
The largest varieties of stabilising control methods considered for nonlinear systems are effective only within a definite domain of the state variables space, titled the attraction domain. The calculation of this domain is generally a delicate task and is time-consuming. This paper offers an efficient computational CPSO approach that is validated to estimate the region of stable equilibrium points for the class of nonlinear polynomial systems.
The main objective of this work consists of designing a swift optimising Lyapunov-based technique that allows the identification of the Lyapunov function coefficients. Subsequently, it exploits the resultant function in a heuristic optimising algorithm to maximise the region of attractions for the class of nonlinear polynomial systems. The synthesised technique is computationally valid and is efficient for real-time applications. In this technique, after a Lyapunov function candidate is computed, a CPSO algorithm estimates the maximal sublevel set of the Lyapunov function such that the time derivative is negative definite through the defined sublevel set. The proposed technique is implemented to estimate the attraction region of numerous benchmark nonlinear systems, which have been previously examined in the related literature, to confirm its aptitude in an analytical study in comparison with pre-existing approaches.
The motivation of this work is to develop an evolutionary synthesised technique that will be appropriate for estimating enlarged DAs of polynomial nonlinear systems. Unlike most existing works, the designed algorithm is running with non-fixed LF. Indeed, the CPSO is currently optimising the LF coefficients and continuously running a quasi-convex LMI optimisation problem. Moreover, the GEVP is considered to be the objective function for the CPSO, which is an original technique of estimation where a meta-heuristic approach is combined with a deterministic technique.
As a matter of fact, it is necessary to synthesise a computationally operative algorithm that converges to a maximised DA in a significantly reduced time. Regardless of the benefit of contributing in a better enlarged accurate DA at less time, CPSO can be very useful for real-time control problems. It is likewise valuable for the control strategies using online sequential composition formalism as described in the literature. For assessing these algorithms, certain benchmark functions are deployed. The simulation results highlight that applying the quadratic Lyapunov function rather than using random sequences enhances the performance of evolutionary algorithms.
After this introduction, the paper is organised as follows. In Section 2, we give a general description of the class of systems deliberated and recall some fundamental notions. In Section 3, we set up the LMI-based iterative algorithms to devise the Lyapunov function which maximises the DA estimate. Section 4 discusses the CPSO-based parameter selection approach. Section 5 presents the simulation outcomes of a renowned numerical example. A comparative study with other methods and a results analysis and discussion technique will be established to evaluate the efficiency of the synthesised strategy in Section 6. A brief conclusion is offered in Section 7.
Standard notation has been used in this paper, where R represents the set of real numbers, A T denotes the transpose of the real matrix A; A > 0 (A ≥ 0) positive definite (semi-definite) matrix; I n is the identity matrix n × n; det (A) represents the determinant of A; A ⊗ B refers to the Kronecker product of matrices A and B; s.t. is used to denote subject to.

LEDA and Basic Notations
Given a nonlinear system represented by the following state equation: wherein x ∈ R n is the state vector U(x) is the control input, F(·) and G(·) are smooth functions describing the system dynamics. Model (1) can be described by the autonomous form written as: .
where the system state vector and the components of f (x) are polynomials in x such that f (0) = 0. The origin should be the equilibrium point under investigation, which is asymptotically locally stable. This provides a compact description for the systems under study. Let ψ t, x 0 ∈ R n be the solution for (2) in x(t). The region of asymptotic stability of the equilibrium point is the set.
Let V(x) ∈ R represent a positive definite radially unbounded LF, which proves that the system origin is asymptotically locally stable. Then, the set is an estimation of the region of stability if ϑ(c) ⊂ Φ, where Φ = x : is the LF V(x) time derivative among the trajectories of (2).
Theorem 1 [26]. A closed set ϑ ∈ R n , which includes the original as equilibrium, may estimate the DA for the system (2) origin provided: 1. For system (2); ϑ is an invariant set.
2. V(x) which is a positive definite function can be established such that · V(x) is negative definite within ϑ. ϑ(γ) represents the Largest Estimated Domain of Attraction and In fact, for sufficiently small neighbourhoods around the origin .
V(x) is negative since V(x) establishes the asymptotic stability of the equilibrium point. Hence, in the domain such that . V(x) is negative definite, the biggest set ϑ(γ) included is attained for the smallest γ s.t. ϑ(γ) reaches a state x at which . V(x) vanishes. To obtain the maximal enlarged region of attraction a specific optimisation problem is formulated in this work. This latter specifically involves a tangency condition between the constraints on the sign the Lyapunov function and the level sets. Such constraints help in avoiding a big number of potential false solutions of the nonlinear optimisation model. The key idea consists of searching the optimal level set of a Lyapunov function which is fully included in the region of negative definiteness of its time derivative. Performing the optimisation will make it possible to determine the tangent point between V(x) and its derivative . V(x). As can been seen in the sequel of the analytical description of this paper, this particular point is providing the maximal region of the attraction region which satisfy the prescribed constraints given the stability Lyapunov theory.

Representattion of Polynomials with Complete Square Matricial Form
The comprehensive set of polynomial possible representations of expressed in a quadratic form is provided by the complete square matricial representation (CSMR) (refer to [24], where homogenous forms have the CSMR similarly defined; the CSMR is also referred to as the Gram matrix). Let the degree of a polynomial p(x) be less or equal to 2m. This is achieved when avoiding linear terms and constants. Let x {m} ∈ R ς(n,m) be the vector containing all monomials of degree less or equal to m in x (i.e., without redundant elements in the mth Kronecker power of vector x): where the dimension ς(n, m) of vector x {m} is given by: Computing the p(x) complete square matrix form with respect to the vector x {m} leads to: wherein: P(α) = P + L(α) (10) and the constant element P ∈ R ς(n,m)×ς(n,m) is any appropriate symmetric matrix such that p(x) = x {m} T P(α)x {m} , x {m} T denotes the transpose of x {m} and α ∈ R τ(n,m) is a vector of free parameters and L(α) is a linear parameterisation of the set: with: For further details on the coefficients ς(n, m) and τ(n, m) the reader can consult Ref. [24]. The matrices P(α) and P are respectively called complete square matrix form and square matrix form of the polynomial p(x). It comes out that p(x) is a sum of polynomial squares if and only if there exists α such that P(α) ≥ 0 or P(α) > 0. Therefore, a sum of polynomial squares p(x) is said to be definite if p(x) vanishes only for x = 0 n .

Calculating Lower Bounds Dependent Parameter
The initial step towards determining the domain of attraction is the introduction of a parameterisation for the square Lyapunov form V(x; G) representing the system origin Equation (2). The conditions ensuring that V(x; G) is a LF for the system origin can be stated as: V(x; G) is radially unbounded and positive definite; .
V(x; G) is negative definite locally. For a homogenous Lyapunov function, the square matricial representation (SMR) is expressed as where x {η v } ∈ R d is a vectorial function comprising all monomials wherein the degree η v in x and G = G T ∈ R d×d is a symmetric matrix comprising the parametric coefficient of the Lyapunov function.
where θ ij ; i, j = 1, . . . , d represents the weighting coefficients for the Lyapunov function. Let the ellipsoidal set associated with the positive definite QLF be defined as while the negative time derivative region is defined as Then, under that condition that ϑ(P, γ(pG)) ⊆ Φ(G), it can be said that ϑ(G, γ(G)) is an ellipsoidal estimate of the DA.
To handle the estimation problem, it should first be converted to an optimisation problem, where the LEDA should be minimised and is represented by γ(G) where This note considers that the problem is the calculation of the largest achievable LEDA in the considered matrix G.
Definition 1 Considering a feasible positive definite QLF V(x; G) for system (2), the corresponding LEDA of the origin is stated by The optimal positive definite QLF is described as the QLF that maximises the DA volume. and It may be observed that the calculation of optimal positive definite QLF is akin to solving a double non-convex optimisation problem. In fact, Ω(G), which represents the volume function, reveals not only the local maxima, but also the global, which is represented by Ω(G * ). Additionally, even the evaluation of Ω(G) necessitates the computation of γ(G), which refers to the solution for the non-convex distance problem (17).
To be able to address (19), it is evident that a specified scalar γ is a lower bound of γ * , i.e., γ < γ * is easily provable if This condition can also be outlined in the following lemma.

Lemma 1:
Let γ ∈ R be a given scalar and define the polynomial where e(x) is any positive definite polynomial, then γ ≤ γ * To check the condition in Lemma 1, the following polynomials are introduced The polynomial degrees of V(x; G) and · V(x; G) = d f (x) are 2η v and η d respectively. If we choose e(x) degree to be 2η e such that Then the polynomial r(x, γ, e(x), V(x; G)) = · V(x; G) + (γ − V(x; G))e(x) has a degree of 2η r , which is also the degree of To get an appropriate form of the optimisation problem, we use the complete square matricial representation (CSMR) and the square matricial representation (SMR) of the polynomials [24]. This implies that V(x; G) , e(x) and r(x, γ, e(x), V(x; G)) may be written using the SMR with regards to the vectors x {η v } , x {η e } , and x {η r } , that do not have the constant term, that is where G, E and R(α, γ, E, G) denote symmetric matrices having suitable dimensions. The matrix R(α, γ, E, G) is specified as where D f (α, G) denotes the CSMR matrix for d f (x), while E and Q(E, G) represent any SMR matrices of, e(x) and V(x; G)e(x), respectively , γ multiplies the parameters of E, (28) leads to a non-convex problem. To move past this difficulty, (29) is reformulated in a GEVP. To proceed, the following preliminary result is required. [24]. Suppose that there exists a polynomial

Lemma 2
where µ is a scalar and the SMR matrix is specified below: Theorem 2 [24].
Function having G > 0 and let µ denote any scalar. The lower bound γ * in (29) is specified as where t is the solution of the GEVP.
Using this GEVP technique, we can establish a domain of attraction. Nevertheless, this technique is usually limited to systems and arbitrary Lyapunov functions that can be expressed by polynomial equations. This paper presents an alternative by employing the evolutionary approach to estimate the optimal Lyapunov function.

Main Results
This section presents a strategy for choosing good parameters for the optimal positive definite QLF. The fundamental idea is to determine the matrix G, which maximises the volume of an ellipsoid having a fixed shape and is included in Φ(G) in (16), which denotes the negative time derivative area. It is reasoned that to maximise the size of the LEDA, an appropriate strategy is to broaden the set Φ(G), which itself contains the LEDA.

Computation of Optimal Positive Definite QLF: Evolutionary Approach
of the ellipsoid specified by ϑ(G, γ(G)) is proportional to (γ(G)) n det(G) , and therefore the optimisation problem (17) searches for the biggest ellipsoid inside the region Φ(G) for some feasible QLF matrix G.
The coefficients for parameters G of the Lyapunov function are evaluated using evolutionary techniques. The user-defined candidate Lyapunov function, as well as the corresponding domain and the condition specified by Theorem 1, are validated for V(x; G) being positive and · V(x; G) being negative. In the case that these conditions are not met, GEVP (34) optimisation is not feasible; hence, we estimate repeatedly the coefficients for the Lyapunov function until the GEVP optimisation has a solution. The basis for this criterion is that to obtain a large optimal estimate of the parameters G using Evolutionary Algorithms, the volume (γ(G)) n det(G) should be maximal. Considering the evolutionary algorithms, we specifically used swarm optimisation and chaotic particle swarm optimisation.
The algorithm begins with V(x; G) being initialised as the V(x; G) = V 0 (x; G 0 ) = γ 0 initial domain. Our algorithm relies on determining the matrix G that maximises the volume of a fixed-shape ellipsoid whose boundary falls inside the negative time derivative region (16).

Motivation on PSO Methods
Swarm intelligence is a crucial subject of evolutionary computing research. It investigates the mutual natural behaviour of distributed, self-organised processes. The stimulation often originates from natural, biological organisms [27]. An interesting example of swarm intelligence is the PSO. This latter is a meta-heuristic global optimisation technique that is now one of the most commonly used optimisation methods.
In this paragraph, a short comprehensive exploration of PSO is presented. As an efficient optimising technique, PSO has seen several rapid advances, such as fuzzy PSO, chaotic PSO, hybridisation with genetic algorithm, artificial bee colony, discrete optimisation, extensions to multi-objective optimisation, etc. [27].
One of the most promising fields in which PSO has been implemented is in automation and controlled systems. The PSO algorithm has been adopted by each research domain that follows a definite scheme and has approved swarm intelligence methods for attaining optimal performance. Several proposed control schemes have been based on linear matrix inequalities with parameters that are tuned by PSO [28].
Numerous research subjects have been fully addressed in the literature [29], including supply chain management, floor planning design, weapon target assignment issues, strategic industrial coordination and symbolic regression analysis.
Many other works investigate analytical analysis problems in control theory. It is essentially a matter of solving nonlinear problems that are difficult to address with conventional analytical techniques. The objective of the remaining part of this paper is to pave the way in order to implement a hybrid estimating technique merging between CPSO and an LMI technique.

Definition
PSO optimisation is a population-based technique. The method's computational efficiency and simplicity are the reason it has seen widespread use in several fields, especially for parameter optimisation [30]. In PSO, the individual agents are referred to as "particles". "Position" refers to a vector of the candidate solutions and is associated with every particle. Additionally, a "velocity" operator regulates the approach to a globally optimal solution [31]. The literature shows that PSO is one of the most efficient methods for solving non-smooth global optimisation problems. The key advantages of such a method can be summarised as follows: It is a derivative-free technique, unlike comparable meta-heuristic optimisation methods. It is simple in its coding implementation and fundamental concept compared to other heuristic optimisation methods.
It is less sensitive to the features of the objective function when compared to the conventional analytical techniques and other heuristic approaches.
Compared with other heuristic optimisation techniques, PSO has a limited number of parameters: the acceleration coefficients and the inertia weight factor. Equally importantly, the sensitivity of parameters to the optimal solutions is considered to be less delicate compared to other heuristic algorithms.
PSO appears to be less dependent on a set of initial points when evaluated with other evolutionary methods, including that the algorithm convergence is robust.
PSO methods can engender high-quality solutions within a shorter calculation time and with more stable convergence features than other stochastic techniques.
Based on the above characteristics, we opted to merge CPSO technique with an LMI-based method to obtain more accurate DA in less time. This can be very beneficial for real-time control problems. It is similarly constructive for control techniques using online sequential composition concepts. A comprehensive presentation of the PSO algorithm, variants of the algorithm, algorithm convergence analysis, and parameter tuning is specified in the references [32][33][34][35], and references cited therein. "Position" refers to a vector of the candidate solutions and is associated with every particle. Additionally, a "velocity" operator regulates the approach to the globally optimal solution.
In the case pertaining to this study, it is only necessary to ensure that a connection exists between the "Stability world" and the "PSO world" using the position of a particle [36]. The decision variables that represent the parameters to be evaluated are specified as the position vector of the ith particle, as presented below: where p is the swarm size (number of particles).
In every iteration k, the velocity for every particle in the swarm is updated, and the positions are calculated using the following equations: k denotes the current candidate solution encoded using the position of the ith particle during iteration k and G (i) k+1 denotes the updated particle position. G (i) k ∈ [L k , U k ], 1 ≤ k ≤ N where L k and U k denote the lower and upper bounds for the n th dimension, respectively. v pk denotes the best position, also called the "personal best", which the ith particle attained in the past. G (i) gk denotes the best position between the neighbours of every particle and is also referred to as the "global best". r 1 and r 2 denote two random numbers in the range (0, 1) and c 1 and c 2 denote the coefficients of acceleration, which regulate how far a particle will reach during a single generation. The inertia weight w k determines the effects of the previous particle velocity on its current velocity. k max defines the maximum number of iterations. Generally, the weight of inertia is decreased in a linear sense from 0.9 to 0.4 during the progression of the exploration. This is done to appropriately balance the global and local searching of the swarm capabilities. The corresponding equation for inertia weight w k can be specified as: In Equation (39), w max and w min have values of 0.9 and 0.4, respectively. Iteration max refers to the maximum iterations permitted.

PSO Stability
The conditions necessary and sufficient for swarm stability were derived in [34] and are stated below and Such conditions ensure that the system converges to a stable equilibrium. Nevertheless, it cannot be said with absolute certainty that the proposed solution would be the global optimum.

Computation of Fitness Function
The fitness function computation is the primary operation that needs to be conducted while executing a PSO algorithm. Accordingly, the evaluation of the particle quality is done using the volume of the domain of attraction. Ω G (i) = (γ(G (i) )) n det(G (i) ) is the objective function that needs to be maximised. Given the position G (i) for the ith particle, the domain of attraction γ G (i) , and the matrix G (i) , the process specified in Algorithm 1 may be used to evaluate the fitness function for the ith particle. (27);

4.
Run an GEVP feasibility test to compute the domain of attraction γ i G (i) .

Fundamentals of Chaotic PSO (C-PSO)
In the PSO technique, the varying parameters w k , r 1 and r 2 are the crucial features influencing the characteristics of the algorithm convergence [37]. The weight of inertia regulates the steadiness between local searchability and global examination. While a larger inertia weight favours global exploration, a smaller weight of inertia favours local searchability. Therefore, the search process typically involves decreasing linearly the weight of inertia, starting at 0.9 and ending at 0.4. Given the fact that the logistic maps are the often-used chaotic behaviour maps, chaotic sequences may be rapidly created and easily saved. Additionally, there is no requirement to store long sequences [38]. The CPSO technique comprises sequences produced by the logistic maps which update randomised parameters r 1 and r 2 in PSO.
The logistic map modifies the parameters r 1 and r 2 as per the specified Equation (37), such that r 1 will be substituted by Cr k where: r 2 will be substituted by: At the value of α = 4, Equation (42) outputs a chaotic sequence, having a value between 0 and 1. In Equation (42), for every independent run, Cr 0 is generated randomly. Figure 1 depicts the chaotic value Cr generated through a logistic map with Cr 0 = 0.001 and 300 iterations.
For CPSO, the velocity update formula may be specified as: The following pseudo-code describes the computational process for calculating the domain of attraction in a polynomial system by using the proposed Chaotic-PSO-based algorithm. Initialise the particles swarm in a random way 03: Create Cr 0 randomly 04: while (the ending criterion is not obtained or iteration number is performed) 05: Assess the particle swarm fitness criterion 06: for p = 1 to particles' number 07: Determine G (i) pk 08: Determine G (i) gk 09: for k = 1 till particle dimension number 10: Refresh the chaotic Cr value as described by Equation (42)  11: Refresh the particles position as defined by Equations (37) and (38)  12: Increment k 13: Increment p 14: Refresh the value of the inertia weight given by Equation (39)  The flowchart describing the developed strategy is described by Figure 2.

Numerical Examples
In this section, we consider certain examples from previous research for which DA has been determined using Lyapunov quadratic functions to illustrate the effectiveness of the presented approach [10,26,39,40].

Example 1
Consider the Van der Pol equation's state-space representation, given below [26].
The state vector here is represented by x = [x 1 x 2 ] T . We employ the GEVP optimisation technique to estimate the initial DA of the steady equilibrium.
It is necessary to obtain a quadratic approximation of the asymptotic stability area around the origin. To achieve this, we examine the GEVP structure in (34) when the assigned Lyapunov function with degree η v = 1, which defines that the LEDA shape is chosen as This implies that V(x; G) can be written using the SMR with regards to the vectors To verify Theorem (2), let us use the polynomial It is noteworthy that D f (α, G), indicating that the CSMR matrix of d f (x) is the corresponding vector of free variables.
where the free variables are computed with the values (8) and (12). Let us note that the degree of We can choose the degree of e(x) as We select η e = 1 e(x) = e 11 x 2 1 + 2e 12 x 1 x 2 + e 22 To calculate Q(E, G) (Equation (31)) the SMR of the polynomial Q(x, e(x), V(x; G)) , we choose the matrix K as Lastly, we compute the Q(E, G) SMR of the product V(x; G)e(x) For this purpose, we begin by encoding the variables to be determined as a position of every particle as per following matrix: The fitness value can be given by the maximum γ * for which there is a viable solution of the GEVP optimisation method as explained in the pseudo-code in Section 4.
The variables' global optimisation G (i) k is carried out by the CPSO operators, the location and the velocities.
The size of swarm is p = 50 and the stopping condition is the highest number of iterations, k max = 100.
To finish this step, by solving (34) with µ = 0.1 we get: After solving this optimisation issue with G * , we get the optimum domain of attraction as given below: The real value of the volume is: Ω(G * ) = 2.1118 (61) Figure 3 depicts the evolution of the CPSO process for 100 iterations. It can be clearly seen that the convergence of the algorithm takes place in less than 10 iterations. This shows the high performance of the developed method in computing the exact theoretical volume value Ω(G * ) = 2.1118. The above result, which provides the volume value Ω(G * ) = 2.1118, helps obtain the largest estimated DA. This result is given by Figure 3. The approximated DA of the equilibrium with γ(G * ) = 1.8770 is represented by the blue ellipsoid in Figure 4. However, the dashed red line, determines the boundary of the region where . V(x, G * ) < 0. Herein, one can accurately check the validity of the computed domain throughout the Lyapunov theory and its specified conditions.

Example 2
This instance is taken from Hahn [39]. It is necessary to find the "best" quadratic approximation of the stability area around the system's origin. The system is as follows: .
The precise stability area is known to be x 1 x 2 < 1. In this area, the origin is depicted asymptotically. The LF describing the LEDA shape is chosen as As the degree η d of .
V(x; G) amounts to 4, we can choose η e = 1 which entails η r = 2. Vectors x {η e } , x {η v } and x {η r } are chosen as Thus, To compute Q(E, G) (Equation (31)), the SMR of the polynomial Q(x, e(x), V(x; G)). We select the matrix K as We obtain Finally, we compute the Q(E, G) SMR of the product V(x; G)e(x).
(θ 11 e 12 − θ 12 e 11 ) 0 0 0 (θ 11 e 12 − θ 12 e 11 ) (θ 11 e 22 + 4θ 12 e 12 + θ 22 e 11 ) (θ 22 e 12 + θ 12 e 22 ) 0 0 0 (θ 22 e 12 + θ 12 e 22 ) θ 22 e 22 The dynamics of the synthesised optimisation algorithm is showed in Figure 5. Compared to the results presented in [39], this demonstrates the accuracy and swiftness of the CPSO technique established in this paper. Figure 6 describes the maximised DA of the equilibrium. The blue ellipsoid represents the DA estimate with γ(G * ) = 5.5715, the dashed red line determines the boundary of the light blue area and represents the region in where · V(x) < 0. A high accuracy can be concluded regarding the obtained performance.
To consider the variables of concern for the matrix G, we reformulate the polynomial with its SMR as As the degree η d of .

(87)
The SMR of the product Q(E, G) In Figures 7 and 8, it is straightforward to check the direct convergence of the algorithm. Indeed, no further fluctuation can be observed, which makes it possible to minimise the convergence time. The steady-state dynamic of the estimated volume is therefore reached after around 65 iterations. A small number of approaches in the literature have studied this kind of system. The technique developed in this paper shows its superiority in providing the maximal volume value in shorter computation time. Figure 8 shows the validity of the obtained DA. Nevertheless, it is theoretically possible to further increase the estimated DA. As a direct explanation for this fact, it can be noted that the obtained domain is not maximal due to the selected LF. The optimal computing of the LF constitutes one of the most advantageous aspects of this work.

Results Analysis and Discussion
This section is dedicated to an evaluative comparative analysis between the synthesised CPSO estimation strategy in this paper and another peer reviewed technique [26]. Table 1 provides the estimated DA features for three dynamical nonlinear systems with quadratic LF(s) taken from the literature [39,40]. The selected examples are presented in their nonlinear polynomial form. Examples E1 [39] and E2 [40] are second-order systems. However, example E3 [39] is a third-order nonlinear polynomial model. Please note that the main evaluation criterion is the domain volume. The obtained results are quite satisfactory in terms of the value of the domain volume. Moreover, the designed strategy provides a complete solution starting from defining the optimal LF and finishing with providing the maximal domain volume. This result is achieved for all studied examples with an easy implementation concept and low consuming time.
A second step was performed to evaluate the performance of the CPSO strategy, consisting of investigating the domain volume for the three examples based on the same LF that was reported in [26,39,40].
For each example, the maximum possible value of Ω(G * ) is computed by the CPSO method and compared with the outcomes obtained by a peer optimisation-based technique, reported in the literature [26]. It appears obvious in Table 2 that for the three investigated examples, the estimated volumes of the DA(s) attained by the CPSO technique are significantly better than the estimates resulting from optimisation-based approaches. For the case of example E3, the performance of the CPSO scheme is even larger and more accurate. E1 [39,40] .
V(x; G) = 0.2742x 2 1 + 0.3634x 1 x 2 +0.4578x 2 2 + 0.1820x 1 x 3 +0.1892x 2 x 3 + 0.1077x 2 3 ≤ 1 ≤ 12.0849 1.0246 12.5328 Likewise, it can be concluded that the developed CPSO procedure is appropriate for estimating enlarged DA of polynomial nonlinear systems. As a matter of fact, it is computationally operational and converges to the optimal DA estimate in a significantly reduced time. Despite the benefit of offering more accurate DA in less time, CPSO can be very useful for real-time control problems. It is likewise valuable for control strategies that use online sequential composition formalism as described in [26]. Finally, we should emphasise the fact that the work presented in [26] is based on a random procedure that characterises the evolution of the designed algorithm. This means that an intensive simulation study should be performed each time to obtain the optimal solutions. This fact is the main disadvantage for the described approach. For the CPSO strategy, the final estimated domain is the same independent of the simulation settings.

Conclusions
This paper studies the quadratic Lyapunov function calculation, which expands the DA volume estimate for systems with polynomials. The designed scheme presented in this paper has the main objective of computing an optimal LF so as to search the largest sublevel set of the DA. Iterative methods based on CPSO were established to obtain the LEDA lower bound. By leveraging CPSO capacity, a new enhanced PSO by using chaotic maps for global optimisation and encoding the variables of the quadratic Lyapunov function to be determined as particle positions was presented. The main advantage of the established algorithm is that one can evaluate the necessary and sufficient conditions as stated in Lyapunov theory with respect to initial conditions selected for the system state variables. The outcomes of the simulations substantiated that the parameters chosen using CPSO possibly leads to the biggest DA. Moreover, no trade-off between the algorithm computational cost and convergence speed is imposed. The convergence rate is therefore reduced, favouring the online implementation for trajectory tracking purposes of complex nonlinear dynamical systems or hybrid nonlinear systems.
Motivated by the excellent outcomes obtained for the instances that were employed, the designed methodology can be considered for other kinds of nonlinear systems and for real processes, which by definition possess nonlinearities. Consideration of the case of rational LF is also a beneficial aspect to this work.