A Chaotic Krill Herd Optimization Algorithm for Global Numerical Estimation of the Attraction Domain for Nonlinear Systems

Abstract

Nowadays, solving constrained engineering problems related to optimization approaches is an attractive research topic. The chaotic krill herd approach is considered as one of most advanced optimization techniques. An advanced hybrid technique is exploited in this paper to solve the challenging problem of estimating the largest domain of attraction for nonlinear systems. Indeed, an intelligent methodology for the estimation of the largest stable equilibrium domain of attraction established on quadratic Lyapunov functions is developed. The designed technique aims at computing and characterizing a largest level set of a Lyapunov function that is included in a particular region, satisfying some hard and delicate algebraic constraints. The formulated optimization problem searches to solve a tangency constraint between the LF derivative sign and constraints on the level sets. Such formulation avoids possible dummy solutions for the nonlinear optimization solver. The analytical development of the solution exploits the Chebyshev chaotic map function that ensures high search space capabilities. The accuracy and efficiency of the chaotic krill herd technique has been evaluated by benchmark models of nonlinear systems. The optimization solution shows that the chaotic krill herd approach is effective in determining the largest estimate of the attraction domain. Moreover, since global optimality is needed for proper estimation, a bound type meta-heuristic optimization solver is implemented. In contrast to existing strategies, the synthesized technique can be exploited for both rational and polynomial Lyapunov functions. Moreover, it permits the exploitation of a chaotic operative optimization algorithm which guarantees converging to an expanded domain of attraction in an essentially restricted running time. The synthesized methodology is discussed, with several examples to illustrate the advantageous aspects of the designed approach.

1. Introduction

The main problem of theory of optimization consists in selecting an optimal vector within a given space of search. The attained vector can minimize or maximize an objective function which is expected to offer the optimal solution. Mostly, advanced intelligent techniques are developed to deal with the design of these kinds of optimization processes. In view of their methodology, optimization techniques can be classified into several main categories: (1) analytical or deterministic, (2) heuristic or random advanced methods and (3) Multi-Objective or single objective optimization [1]. The first class, based mainly on the gradient theory, are known to be “strict step” methods. Indeed, identical solutions are obtained whenever the initial starting conditions of the algorithm are the same. However, heuristic algorithms are founded on a random walks principle. As a result, the optimization approaches cannot be reiterated under any initializing conditions. Nonetheless, both classes mostly provide satisfactory solutions leading to similar final optimal results. Lately, bio-inspired algorithms demonstrate an impressive and valuable performance in solving challenging nonlinear optimization problems [2].

All meta-heuristic methods, to some degree, manage to identify a balance between local and global searches (intensification and randomization) [3]. To solve high-dimension nonlinear optimization problems including task-resource assignment, these flexible nature-inspired meta-heuristic methods are used. These approaches can fully use the entire data pertaining to the population in order to identify precise solutions. Evolution theory has received a lot of attention to date. Genetic algorithms (GAs) [4,5] are based on the gradient-free approach that mimics evolution. Since then, several meta-heuristic nature-inspirated techniques have been formulated; these include particle swarm optimization (PSO) [6,7,8], evolutionary strategy (ES) [9], firefly algorithm (FA) [10], ant colony optimization (ACO) [11], differential evolution (DE) [12], probability-based incremental learning (PBIL) [13], big bang–big crunch algorithm [14], biogeography-based optimization (BBO) [15], harmony search (HS) [16], cuckoo search (CS) [17], animal migration optimization (AMO) [18], krill herd method (KH) [19], bat algorithm (BA) [20], teaching–learning-based optimization (TLBO) [21], and charged system search (CSS) [22]. KH is a modern swarm intelligence optimization technique inspired by krill herding behaviour [23]. In KH, the krill movement objective function is determined by the krill distance, krill swarm density, and food density. Krill location comprises three mechanisms: (i) movement due to other krill, (ii) foraging movement, and (iii) physical diffusion. The KH algorithm is beneficial because of its simplicity, facilitating easy implementation in parallel computing [24].

Additionally, KH is a new population-dependent swarm intelligence algorithm [24] considering the Lagrangian and evolutionary activities of krill in their natural environment. This algorithm can help understand and use krill behaviour for optimizing real-world problems. Although the KH based algorithm is convenient for many optimization problems, it cannot avoid local optima and therefore is inappropriate to determine a global optimal solution [25]. Several techniques are used in the literature to enhance the fundamental KH technique and improve its performance [26].

In order to address this concern, Wang et al. [27] suggested using an enhanced algorithm based on chaotic patterns. Chaotic maps have been implemented for preventing local optima and directly reaching the globally optimal solution [24]. Furthermore, this technique requires less computation time to reach a solution. Hence, it is feasible to avoid local optima and provide rapid convergence, which makes this technique usable in real-world applications for addressing optimization problems with distinct conditions.

This study introduces the Chaotic KH-based technique (CKH) for increasing the KH convergence rate. Numerous one-dimensional chaotic maps replace KH parameters. To generate chaotic sequences efficiently and rapidly, Chebyshev maps has been implemented in the CKH algorithm. The suggested technique is tested on benchmark problems in order to estimate the domain of attraction (DA). It is expected that CKH will provide better performance than the techniques like those developed in [28,29], specified in the literature. The primary reason for performance increase is the use of deterministic chaotic patterns instead of linearly declining values.

Stability is the most critical factor in the study of nonlinear dynamic systems. For most applications, knowing the asymptotic stability associated with an equilibrium position is typically inadequate because it is essential to consider the maximum distance between the initial state and equilibrium position to ensure asymptotic stability. Consequently, the principle of the domain of attraction was born. Since the publication of the stability theorem in 1892, the challenge of precisely evaluating or estimating the DA has remained unresolved. In general, evaluating the DA is very complex [30,31]; hence, the problem concerning the computation of the most significant invariant DA subset is also taken into account. Lyapunov and non-Lyapunov methods are the two types of DA estimation methods [30].

The maximal Lyapunov function (LF) can be used to calculate the precise DA of a nonlinear system. A rational LF can be used to approximate such an expression [31,32]. Lyapunov quadratic functions typically compute stability area approximations conservatively. Regardless, quadratic LFs are widely employed for stabilization [33]. Moreover, these functions are used to build asymmetric LF used for forecasting DA associated with linear systems concerning asymmetrical extrinsic action [34].

In contrast, research indicates that random LFs fulfilling Lyapunov constraints offer DA estimates concerning control systems of the nonlinear polynomial type [35]. The set of LFs appropriate for computing DA was augmented using a proposition proved in [36]. A novel technique for computing LFs was formulated; its degrees of freedom were higher than the traditional technique [37]. Imprecise polynomial expressions concerning non-polynomial systems must be estimated in order to apply Lyapunov polynomial functions for computing the DA specific to non-polynomial systems [38]. Several approximation techniques that rely on the multi-parameter polynomial method to assess residual expression were formulated in [39]. The DA can also be computed using piecewise affine functions used as LFs [40,41,42].

Presently, nonlinear system stability and formulation of new DA estimates typically rely on the Pupov criterion, The Zubov technique [43], LaSalle’s theorem, the circular criterion, and the framework of central manifolds [29].

In a majority of the situations, the problem concerning DA estimation is simplified to a convex or non-convex optimization problem [44]. This is addressed using optimization techniques such as SOS [45], LMI [46], intelligent optimization techniques [47], integration of the genetic algorithm and LMI [4]. In comparison to optimization approaches, a swift and computationally efficient technique formulated for real-time applications relies on an algorithm that identifies the most-extensive sublevel set pertaining to the LF under the constraint its time derivative is negative [28].

Predicting DA pertaining to inexact scenarios is highly challenging, specifically for situations comprising large scale systems [48]. An internal DA approximation pertaining to a type of inexact nonlinear systems can be determined using a novel technique that LFs independently of parameters [49]; moreover, the branch and bound technique is also used [50]. Furthermore, researchers [51] proposed an updated technique to compute sliding mode control stability under uncertain conditions. Filling measures have been widely used for stochastic control applications; these measures are required for computing the DA [52].

DA approximation is also performed using the non-Lyapunov trajectory inversion technique that relies on topological aspects. However, computational complexity is one challenge associated with this technique; this can be addressed by integrating the approach with Lyapunov techniques [53,54]. Formulating reachable sets is a different approach sometimes used for evaluation of DA. Formulating reachable sets using an approximate calculation method is specified in [55].

This paper suggests a metaheuristic scheme for Lyapunov-based methods to determine DA corresponding to several nonlinear systems. This technique is intended to have better computational effectiveness and to be helpful for real-time implementations. When a potential LF is selected, an evolutionary algorithm explores the biggest sublevel LF set under the constraint that the derivative with respect to time is negative for all values in the sublevel set. A heuristic optimization scheme is set up, comprising a tangency restriction concerning level set and LF sign requirements. These restrictions facilitate skipping numerous potential local solutions specific to the nonlinear optimization process.

In [56] a hybridized monarch butterfly optimization method has been proposed. Two state-of-the-art swarm intelligence approaches have been implemented and the designed technique was enhanced to be applied to a convolutional neural network scheme problem. In [57], the authors created a general model, the adaptive network-based fuzzy inference system code, that has been applied for estimating nanofluids’ relative viscosity. A statistical analysis has proved the precision of the synthesized model.

The main contribution of this paper consists in combining the CKH algorithm with position updating equations in a chaotic map to compute the best solution related to the DA estimation problem. The implemented CKH algorithm in this work will guarantee a high convergence speed in computing the radius of the maximal estimated DA. Moreover, the coordinates of the tangency point between the LF and its derivative will be computed and efficiently determined. Along with this valuable performance, a minimal optimization routine time budget will be insured, favoring trajectory tracking engineering control problems. To obtain the maximal expanded DA, a particular optimization problematic is investigated in this work. This later includes a delicate tangency constraint related to the LF sign and its derivative. The challenging objective consists in avoiding important potential dummy and local solutions for the applied optimization routine. The designed strategy will permit the tangency point to be fixed accurately and converge to the solution in a reduced running time. The assumption behind the investigated research is that the designed method could tackle problems better at hand than peer methods; indeed, no efficient procedure has been formulated to determine this particular point.

This paper discusses a novel approach for modifying the shape and increasing the size of the DA specific to nonlinear systems. Section 2 details how iterative techniques are set up to formulate the LF that optimizes DA. Section 3 details the CKH-specific parameter choice. Section 4 details simulation results of prevalent numerical problems. Section 5 presents a comparison of the proposed algorithm with other techniques; additionally, there are analyses of the results and discussion concerning algorithmic efficiency assessment. Summarizing conclusion and perspectives are provided in Section 6. An appendix comprising an overall discussion concerning the class of systems and several basic notions has been provided at the end of the paper.

2. Estimation of the Domain of Attraction

2.1. Fundamentals of the Domain of Attraction

The constraints defined in Theorem 1 (see Appendix A) will be used with a Lyapunov to compute the inner approximations corresponding to this set. The entire section is based on the assumption that $V:R^n\to R$ a is positive and differentiable function satisfying $V\left(x_0\right)=0$ and $V\left(x\right)>0$ when $x\ne x_0$. DA approximations obtained using Lyapunov-based schemes typically use statements like [30], if for positive $c$ and $\epsilon$.

$$V\left(x\right)\le c\Rightarrow \frac{\partial V\left(x\right)}{\partial x}f\left(x\right)\le -\epsilon V\left(x\right)$$

(1)

The following expression defines the set correlated to the LF:

$$\Omega =\left\{ x\in R^n:V\left(x\right)\le c\right\}$$

(2)

The negative time derivative area is specified as:

$$\mathrm{\Upsilon }=\left\{ x\in R^n:\dot{V}\left(x\right)<0\right\}\cup \left\{0\right\}$$

(3)

Then, $\Omega$ is an approximation of the DA corresponding to the origin, such that $\Omega \subseteq \mathrm{\Upsilon }$.

An unsophisticated initial technique for broadening this analysis requires being disjoint from the area specified by:

$${\displaystyle \sum =}\left\{ x\in R^n:\dot{V}\left(x\right)\ge 0\right\}$$

(4)

We try to determine an area specified using a set of smooth “Derivative LF” $\dot{V}\left(x\right):R^n\to R$ where no failure areas are considered. The points corresponding to $\dot{V}\left(x\right)=0$ are known as “barriers”.

Hence, we can determine a solution that is true for all x in the following expression:

$$x\in R^n\in \Omega \cap {\displaystyle \sum =\left\{x\left|V\left(x\right)\right.\le c, \dot{V}\left(x\right)\ge 0\right\}}$$

(5)

this is a part of the area of attraction and is also positively invariant.

Put differently, these constraints evaluate the “safe set” ${\displaystyle \sum }$ created using the intersection of functions meeting $\dot{V}\left(x\right)=0$, while the c-sublevel set corresponding to $V\left(x\right)$ is the latest DA inner approximation.

The present study is based on a function $V\left(x\right)$ that is quadratic with respect to the state $x$.

$$V\left(x\right)=x^{T}Px, P=P^T\in R^{n\times n}, P>0$$

(6)

We assume that $P$ has two dimensions, which helps simplify the presentation. The outcomes can be generalised for bigger matrices. It is supposed that:

$$P=\left[\begin{array}{cc}{p}_1& p_2\\ p_2& p_3\end{array}\right]$$

(7)

Then

$$V\left(x\right)=p_1{x}_1^2+2p_2{x}_1{x}_2+p_3{x}_2^2$$

(8)

This form of a LF is used for specifying a general ellipsoid present in the $(x_1,x_2)$ plane, as described by the equation:

$$\Omega =\left\{x\in R^n:V\left(x\right)=p_1{x}_1^2+2p_2{x}_1{x}_2+p_3{x}_2^2=c, c>0\right\}$$

(9)

The DA is the region that comprises the convergence of all system trajectories to a particular equilibrium position, as specified using Equation (9).

2.2. Concept Illustration for Enlarging the DA

The concept objective is to determine the guaranteed large value of DA specific to a particular nonlinear polynomial system corresponding to a LF and research the state space using a meta-heuristic technique. This concept is illustrated in Figure 1.

2.2.1. Research on Maximum DA

As per (9), a greater level set value c corresponds to relatively precise DA estimation. The computation of the LF’s maximum level set relates to DA estimation. It can be determined by solving the pseudo-optimization as specified below:

$$\left\{\begin{array}{l}\underset{c, x}{\max } c\\ s.t.\left\{x belong to level set V\left(x\right)-c=0\right\}\end{array}\right.$$

(10a)

$$\left\{x is any point belonging to region \Omega contained in D\right\}$$

(10b)

The basis of problem (10) is to determine the maximum level set corresponding to $V\left(x\right)$ (10a), which lies entirely within the definitive negative region of $\dot{V}\left(x\right)$ (10b). It may be obtained by determining the globally optimal solution for the following problem:

$$\left\{\begin{array}{l}\underset{c, x}{\min } c\\ s.t. V\left(x\right)-c=0\\ \dot{V}\left(x\right)=0\\ c>0\end{array}\right.$$

(11)

The intent behind using problem (11) is to determine the minimum level set corresponding to function $V\left(x\right)$, as contained in the level set $\dot{V}\left(x\right)=0$. The required solution is one point located in the state space that correlates to the level sets $V\left(x\right)=c$ and $\dot{V}\left(x\right)=0$ coming in contact. The nonlinear nature of model (11) is highlighted; consequently, it might have several local solutions. To be able to converge to an acceptable estimate, global optimality must be attained. This concept is graphically presented using Figure 2. The optimisation problem defined in (11) is nonlinear, and several solutions can be identified. Nevertheless, most solutions are inappropriate for estimating the DA. The following schematic version is assessed for additional evaluation of this concept.

Figure 2 presents the global optimization problem qualitatively.

The expression mentioned below indicates a hypersurface:

$$dV\left(x\right)/dt=0, x\ne 0$$

(12)

This is related to the area boundary where $\dot{V}\left(x\right)$ is negative definite, where we intend to determine the guaranteed estimation $\Omega$. Such expressions corresponding to QLF denote the inside of the ellipsoid specified using (9).

Distinct points on this figure can be derived as optimization problem solutions (11). The requirements for optimisation (11) apply for points $x_5$ and $x_7$. The LF and its differentiated version verify the sign change corresponding to a particular area in the state space. However, other paths cross the predicted DA transversally. These solutions correspond to local minima and are considered dummy; hence, these are not considered. Furthermore, coordinates $x_1, x_2, x_4, x_6,$ and $x_8$ validate this inference. Considering the intended DA, several transversal intersections can be found.

Nevertheless, it can be observed that $x_3$ is considered as a global minimum in Figure 2. It should be noted that all system paths intersect the computed DA tangentially. It is crucial to verify the below-mentioned constraint in order to specify the optimisation problem (11) in a refined way.

Constraint 1: 

$$\left\{\begin{array}{l}V\left(x\right)=c\\ \dot{V}\left(x\right)<0\\ c>0\end{array}\right.$$

(13)

Constraint 2: 

The optimal solution must correspond to a state-space area prior to a sign change corresponding to $V\left(x\right)$ and its derivative tackles as $\left|x\right|$ approaches infinity.

Constraint 3: 

The level set specified by:

$$V\left(x\right)=c$$

(14)

corresponds to a global minimum.

It is necessary to ensure that these three constraints are met so that “dummy solutions” can be avoided and the determined globally optimal formulation (11) outputs tangential coordinates. To have detailed examples related to this concept, the reader can consult the following reference [58].

2.2.2. Algorithm of Computing the DA

We now specify algorithmic aspects that enhance DA estimates iteratively using numerous calculations of $V\left(x\right)$ that help determine the state $x$ where the first constraint is fulfilled. These cycles enhance the approximation of the radius of the contained ellipse (9).

Before algorithm execution begins, the origin is in an asymptotically stable state. Hence, as specified by the linear theory, there exists a locally applicable QLF:

$$V\left(x\right)=x^{T}Px$$

(15)

where $P$ is a positive definite, symmetric matrix having $n\times n$ dimension. The below mentioned Lyapunov expression is solved to compute $P$.

$$PA+A^{T}P=-Q$$

(16)

corresponding to a positive definite matrix $Q$ having $n\times n$ dimension where $A$ denotes the Jacobian of $f\left(x\right)$ at equilibrium.

The derivative of the LF is computed as:

$$\dot{V}\left(x\right)=\frac{dV\left(x\right)}{dx}f\left(x\right)$$

(17)

State-space sampling facilitates an instantaneous assessment of $V\left(x\right)$ and $\dot{V}\left(x\right)$.

This technique changes lower and upper bounds corresponding to $c$ denoted $c_{\min }$ and $c_{\max }$, respectively. When combined, these values provide a more precise DA estimate. The lower bound $c_{\min }$ is zero, while the upper bound $c_{\max }$ is infinite when the algorithm begins execution. Considering a randomly selected $x^{\left(i\right)}$ and conditions $\dot{V}\left(x^{\left(i\right)}\right)<0$ and $c_{\min }<V\left(x^{\left(i\right)}\right)<c_{\max }$ are met, the lower bound is changed to the LF‘s value, i.e., $c_{\min }=V\left(x^{\left(i\right)}\right)$. On contrast, if $\dot{V}\left(x^{\left(i\right)}\right)\ge 0$ and $V\left(x^{\left(i\right)}\right)<c_{\max }$ conditions are met, then the upper bound is changed to $V\left(x^{\left(i\right)}\right)$. When the sampling process has processed numerous samples, the lower bound increases; however, the increase is not always monotonic. The final convergence point $c$ and most significant sublevel set $\Omega \left(c_{\min }\right)$ is determined. Furthermore, $c_{\max }$ exhibits a monotonic decrease and eventually converges down to $c$.

When the requirements for Theorem 1 are met corresponding to state $x^{\left(i\right)}$, then $V\left(x^{\left(i\right)}\right)$ represents a probable estimate of $c$ therefore, it is retained in an array. This condition ensures that the approximate DA values calculated using $c_{\min }$ satisfy the conditions corresponding to Theorem 1, thereby providing more precise estimates. Array $\Im$ initially contains zero.

Considering that conditions $\dot{V}\left(x^{\left(i\right)}\right)<0$ and $V\left(x^{\left(i\right)}\right)<c_{\max }$ are met, $\Im$ stores the value of $V\left(x^{\left(i\right)}\right)$, i.e., $\Omega \left(V\left(x^{\left(i\right)}\right)\right)$, which is a probable DA estimate. If conditions $\dot{V}\left(x^{\left(i\right)}\right)\ge 0$, $V\left(x^{\left(i\right)}\right)<c_{\max }$ and $c_{\min }\ge c_{\max }$ are met, then the method tries to select a new lower bound from the values stored in the array. Condition $c_{\min }<c_{\max }$ is considered when trying to identify the maximum value from the array. Choosing an existing lower bound value meets the constraint $\dot{V}\left(x^{\left(i\right)}\right)<0$ for the identified sublevel set $\Omega \left(c_{\min }\right)$. The lower bound is zero, considering the worst-case situation.

The sampling technique used in this study is conservative; however, there could be a substantial overestimation of DA. One example is the region $\dot{V}\left(x^{\left(i\right)}\right)<0$ not having a simple connection. Here, the technique might not disregard tiny holes inside the area where $\dot{V}\left(x^{\left(i\right)}\right)\ge 0$. Presently, there is no assurance that this algorithm always converges; however, empirical data obtained using numerous simulations and real-world tests indicate that this method converges to the precise level set if the sample count is large.

Algorithm 1 presents a series of steps for identifying the DA corresponding to the origin of (A1). Ensuring that Lyapunov theory concepts are met, the following steps are followed:

Algorithm 1 DA estimation

Initialization: Initialize the parameters $c$, $c_{\min }, c_{\max }$, $\Im =\left\{0\right\}$

Liniarization: The Jacobian of system (A1) at the origin and determine the matrix LF          $A\leftarrow \partial f∕\partial x\left(x\right){\parallel }_{x=0}$          $P\leftarrow solution of P^{T}A+AP=-I$

$V\left(x\right)=x^{T}Px$

for $i=1:n$ do

Evaluate $V\left(x^{\left(i\right)}\right)$, $\dot{V}\left(x^{\left(i\right)}\right)$ for each state

if $(V\left(x^{\left(i\right)}\right)=c_{\max })$ && $(\dot{V}\left(x^{\left(i\right)}\right)=0)$

$\Im =\{V(x^{(i)})\}$

else if $(V\left(x^{\left(i\right)}\right)<c_{\max })$ && $(\dot{V}\left(x^{\left(i\right)}\right)\ge 0)$

$c_{\max }=V(x^{(i)})$

If $c_{\min }\ge c_{\max }$

$c_{\max }=\max \{c\in \Im ,c\le c_{\max }\}$

end

end

end for

Result: Return the best value of DA

Draw the graphic Domain definite by $\mathrm{\Omega }=\left\{V\left(x^{(i)}\right)=c_{\max }\right\}$;

End.

Our technique relies on ascertaining state $x$ in space where DA radius is maximised, corresponding to a particular ellipsoid having boundaries situated within the negative time derivative area. This technique is not feasible for higher dimensions; moreover, state–space discretisation and several simulations would be needed.

3. Proposed Krill Herd Algorithm for State Assessment

3.1. The Motivation for Using the Krill Herd Method

Extensive research has been conducted to determine the mechanisms that result in the formation of non-random patterns by marine animal populations [18]. Several mechanisms have been understood: safety from predators, feeding, environmental characteristics, and better reproduction [23].

The Antarctic krill is one of the most commonly researched marine species [24]. In fact, the krill herd is characterized by several uncertainties regarding its representative distribution [25]. Various conceptual frameworks have been proposed to explain the pattern of krill herds [26]. The results suggest that krill swarms are the fundamental organizational unit of this species.

Marine predators such as sea birds or penguins attack krill by leading individual krill to an area with lesser krill density. Following the predatory attack, krill herd formation has two primary objectives: (1) enhance krill density and (2) access to food. Krill behaviour to enhance the density and locate food are identified as the objective function. Subsequently, herding is observed around local minima. Individual krill movement is such that the best solution can be found in this search for food and increased density.

3.2. Krill Herd

Gandomi and Alavi first proposed the Krill Herd (KH) algorithm in 2012. KH [27] is a novel meta-heuristic technique for addressing optimisation problems. Krill swarm patterns during food search are the motivation for this algorithm. The position of individual krill is affected by three aspects specified below [27]:

i.

Krill movement due to other individuals;

ii.

Foraging behaviour;

iii.

Random diffusion;

A simple Lagrangian model [27] can be formulated for the three aspects of KH. Equation (18) specifies this model:

$$\frac{d{X}_k}{dt}=N_k+F_k+D_k$$

(18)

where

$N_k$: krill motion due to other individuals;

$F_k$: foraging-specific movement;

$D_k$: is the physical diffusion corresponding to the k^th individual krill;

$k$: denotes a single krill;

$t$: generation count;

3.2.1. Krill Movement Due to Other Individuals

The first movement direction is approximated using three characteristics: repulsive effect, target effect, and local effect. The movement of an individual krill i^th is modelled as:

$$N_k^{next}=N^{\max }{\alpha }_k+{\omega }_d{N}_k^{present}$$

(19)

where

$${\alpha }_k={\alpha }_k^{local}+{\alpha }_k^{t\arg et}$$

(20)

and

$N^{\max }$ denotes the maximum induced speed,

${\omega }_d$ denotes inertia weight in the range [0,1],

$N_k^{present}$ denotes the previous induced motion of the i^th individual krill

${\alpha }_k^{t\arg et}$ and ${\alpha }_k^{local}$ denote the target and local effects, respectively.

3.2.2. Foraging Movement

Two primary aspects determine the second movement: location of food, and prior experience concerning the location of food. This motion can be modelled for the i^th krill as specified below [27]:

$$F_k^{next}=V_f{\beta }_k+{\omega }_f{F}_k^{previous}$$

(21)

where

$${\beta }_k={\beta }_k^{food}+{\beta }_k^{best}$$

(22)

and

$V_f$ denotes foraging speed; for the present study, $V_f$ is set to 0.02 [27].

${\omega }_f$ denotes inertia weight in the range [0,1],

$F_k^{previous}$ denotes previous foraging movement,

${\beta }_k^{food}$ denotes food attraction

${\beta }_k^{best}$ is the influence of best fitness

3.2.3. Random Diffusion

The third motion is fundamentally random. Two factors are used for expressing the model corresponding to this motion: random vector and highest diffusion speed. The model expression is specified below [27]:

$$D_k=D^{\max }\delta V_f$$

(23)

where,

$D^{\max }$ is the maximum induced speed;

$\delta$ denotes the random directional vector [0,1]

3.3. Movement in KH

Typically, the motions specified in KH change often with respect to krill position to attain the best position. The two subsequent motions (2nd and 3rd) caused by other krill comprise two local and two global techniques used concurrently to make KH an efficacious and potent technique. The following equations specify the location of individual krill i in the $\left[t, t+\epsilon \right]$ time interval:

$$X_k(t+\epsilon )=X_k\left(t\right)+\epsilon \frac{d{X}_k}{dt}$$

(24)

where $\epsilon$ is a crucial constant (speed vector scaling factor).

Equation (25) may be used to determine the value of $\epsilon$

$$\epsilon =P_{cst}{\displaystyle \sum _{j=1}^{N_v}(U_j-L_J)}$$

(25)

where,

$N_v$: denotes the total variable count

$U_j, L_J$: denotes the upper and lower bounds corresponding to the jth variables

$P_{cst}$: position constant in the [0,2] range.

To achieve better optimisation and enhanced convergence speed, multiplication and crossover aspects specific to the algorithm are integrated with KH [27].

3.3.1. Crossover

The crossover probability (Cr) parameter regulates the crossover process. The position of every individual krill is changed based on its interaction with other krill. The jth component corresponding to the ith krill may be specified as [27]:

$$X_{i.j}=\left\{\begin{array}{l}{X}_{m,j} if rand\le Cr\\ X_{i,j} if rand> Cr\end{array}\right.$$

(26)

where, $m=1, 2, \dots , N$ and $N$ denotes krill population and $z_i^{best}$: is the best historical meeting point for krill i.

3.3.2. Mutation

The mutation probability (Mr) parameter regulates the mutation process; it may be specified as [27]:

$$X_{i.j}=X_{best,j}+\theta (X_{s,j}-X_{t,j})$$

(27)

where,

$X_{best,j}$: globally best vector,

$X_{s,j}, X_{t,j}$: two random vectors,

$\theta$: scalar in the [0,1] range.

The updated $X_{i.j}$ value is expressed as [27]:

$$X_{i, j}^{\mathrm{mod}}=\left\{\begin{array}{l}{X}_{i, j}^{new} if rand\le Mr\\ X_{i, j} if rand> Mr\end{array}\right.$$

(28)

3.4. Computation of the Fitness Function

Execution of the KH algorithm comprises the determination of the fitness function as its primary step. Correspondingly, particle quality assessment is made using the domain of attraction. $\mathrm{\Omega }=\left\{ x\in R^n:V\left(x^{\left(i\right)}\right)\le c\right\}, c>0$ is the objective function expression that must be maximised. Considering the i^th particle position $x^{\left(i\right)}$, the third step of Algorithm 2 needs to be verified. States $x^{\left(i\right)}$ and the domain of attraction $\mathrm{\Omega }\left(V\left(x^{\left(i\right)}\right)\right)$ can be used as per the below-mentioned process to determine the fitness function corresponding to the ith particle.

3.5. KH Computation Procedure

Algorithm 2 lists the computational steps for this process.

Algorithm 2 The KH computional process

Begin

Initialization. Set the initial values for the population size (Np), the maximum induced speed ($N_{obj}^{\max }$), the foraging speed  ($N_{obj}^{\max }$), the maximum diffusion speed $N_{obj}^{\max }$), the generation counter (m) and the maximum number of fitness function evaluations ($N_{obj}^{\max }$).

Fitness evaluation. Evaluate the fitness function and randomly set up the position n x(i);i = 1,2,…,Np set of each krill individual then assess the fitness function value for all individuals.

Whilem < ($N_{obj}^{\max }$) or the stopping criteria are not accomplished do

Organize the population initiating from the best to worst.

Keep the best krill individual in “STORE”

fori = Np do

Compute the:      a) Physical diffusion      b) Foraging motion      c) Induced motion

Implement genetic operators

Update the krill individual position in the search domain

Assess each krill individual based on its updated position

end for

Organize the population from best to worst and fix the current best one. Increment the counter m

m = m + 1

end while

Results: Presentation of the results and visualition;

End.

3.6. Fundamentals of Chaotic KH (CKH)

All metaheuristic algorithms share a similar search scheme divided into two parts: exploration and exploitation. During the exploration phase, the search space is evaluated extensively to generate additional solutions. In contrast, the exploitation phase comprises a rapid convergence of the algorithm to the global optimum. These two steps are sophisticated aspects required to determine the global optimum; here, an accurate solution is typically balanced between these two factors. However, a significantly complex optimisation such as a DA problem with an unknown search space presents challenges. Hence, the most appropriate trade-off between exploitation and exploration cannot be ascertained. In order to focus on these two aspects in the search space, we have considered several adaptive and random variables along with the KH.

Chaos is among the mathematical techniques that have recently been used for enhancing exploitation and exploration [27]. Chaos theory provides the leading method for enhancing evolutionary algorithm performance concerning the ignorance of local optima and enhancing convergence rate; moreover, the chaos technique does not rely on random parameters. Additionally, except for chaotic maps, another fundamental enhancement is the integration of the elitism approach to the CKH framework. Like other optimisation schemes, we add elitism in some form to retain the optimal krill population.

The elitism scheme which is implemented in the CKH algorithm is considered as one of the most important enhancements offered by this. As with peer population-based optimization strategies, some kinds of elitism are typically combined so as to retain the best solutions in the investigated population. In the main cycle of the CKH, to start with, the KEEP best population results are saved in a variable KEEPKRILL. Generally, the KEEP worst solutions are eventually replaced by the KEEP best solutions. This later elitism process can ensure that the whole population cannot be declined to the population with worse fitness than the former. This avoids the best population from being abandoned by the dynamic of the motion calculation operator. Note that the elitism strategy is used to save the feature of the krill that has the best fitness in the CKH process, so even though motion calculation operation degrades its matching krill, this later has been stored and can be recovered by its previous good status if required.

Inertia weights ${\omega }_d, {\omega }_f$ are the primary parameters that regulate the convergence rate of KH. To enhance search efficacy and ensure convergence to the best solution, the chaos concept is integrated with the KH framework for this study’s scope. Considering that Chebyshev maps are the most-extensively used chaotic behavioural maps, it is possible to generate chaotic sequences rapidly and efficiently. Furthermore, there is no need to retain long sequences. Chebyshev maps used during the CKH method change the value of randomised parameters ${\omega }^{chebychev}=\left({\omega }_d, {\omega }_f\right)$ in KH.

Chebyshev maps update parameters ω_d and ω_f in accordance with the Equations (19) and (21):

$${\omega }_j^{chebychev}=\cos \left(j\ast {\cos }^{-1}\left({\omega }_{j-1}^{chebychev}\right)\right)$$

(29)

Equation (29) produces a chaotic sequence in the [0,1] range. For every independent execution of equation (29), ${\omega }_0^{chebychev}$ is output randomly. The chaotic value ${\omega }_j^{chebychev}$ produced using a logistical map with 300 runs and ${\omega }_0^{chebychev}=0.001$ is depicted in Figure 3.

In the case of CKH, Equations (19) and (21) can be written as:

$$N_k^{next}=N^{\max }{\alpha }_k+{\omega }_{d, j}^{chebychev}{N}_k^{present}$$

(30)

$$F_k^{next}=V_f{\beta }_k+{\omega }_{f, j}^{chebychev}{F}_k^{previous}$$

(31)

Algorithm 3 presents the pseudocode for the CKH computation process.

Algorithm 3 The CKH computional process

Begin

Initialization. Initialize the CKH’s parameters m, (Np), ($N_{obj}^{\max }$), ($N_{obj}^{\max }$), ($N_{obj}^{\max }$) and ($N_{obj}^{\max }$)

Set up the chaotic maps value in random way and inertia weights ωchebychev = (ωd, ωf) as well.

STORE: number of the best krill swarm to retain from one generation to the next.

Fitness evaluation. Evaluate the fitness function and randomly set up the position x(i);i = 1,2,…,Np set of each krill individual, then assess the fitness function value for all individuals.

Whilem < $N_{obj}^{\max }$ or the stopping criteria are not accomplished do

Organize the population initiating from the best to worst.

Keep the best krill individual in “STORE”

for i = 1: Np do

Compute the:      a) Physical diffusion      b) Foraging motion      c) Induced motion

Implement genetic operators

Update the krill individual position in the search domain

Assess each krill individual based on its updated position

end for

Update the “STORE” with the new best value;

Organize the population from best to worst and fix the current best one Increment the counter m: m = m + 1

end while

Results: Presentation of the results and visualition;

End.

The flowchart describing the developed strategy is described in Figure 4.

4. Illustrative Examples

4.1. Example 1

The nonlinear dynamic expression below might be assumed as a model for a pendulum [28].

$$\left\{\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-\sin \left(x_1\right)-0.5x_2\end{array}\right.$$

(32)

where the pendulum, having angular velocity $x_2$, is at an angle $x_1$ with the vertical axis. The state vector is specified as: $x={\left[x_1, x_2\right]}^T$. We use the proposed approach with several chaotic map distributions to find an approximate DA corresponding to the stable equilibrium $x=\left(0, 0\right)$. To calculate a potential LF, we linearise the dynamic model for a pendulum. Hence, the objective is to identify an acceptable stability domain:

$$\mathrm{\Omega }=\left\{ x\in R^n:V\left(x\right)\le c\right\}, c>0$$

(33)

where c is a real non-zero constant and $V\left(x\right)$ denotes the positive definite QLF expression specified below:

$$V\left(x\right)=x^{T}P x, P^T=P$$

(34)

P denotes a symmetric and positive matrix that is ascertained by determining the solution of the Lyapunov expression.

$$A^{T}P+PA=-Q$$

(35)

$Q$ is an identity matrix as specified below:

$$Q=I$$

(36)

In this case,

$$P=\left[\begin{array}{cc}2.25& 0.5\\ 0.5& 1\end{array}\right]$$

(37)

The computed LF is:

$$V\left(x\right)=x^T\left[\begin{array}{cc}2.25& 0.5\\ 0.5& 1\end{array}\right]x$$

(38)

which can be expressed as:

$$V\left(x\right)=2.25x_1^2+x_1{x}_2+2x_2^2$$

(39)

In the case the derivative of $V\left(x\right)$, as defined by the expression below, is negative definite, and the nonlinear model has assured stability.

$$\dot{V}\left(x\right)={\dot{x}}^{T}Px+x^{T}P\dot{x}$$

(40)

This equation may be specified as:

$$\dot{V}\left(x\right)=4x_1{x}_2-x_2^2-x_1\sin \left(x_1\right)-4x_2\sin \left(x_1\right)$$

(41)

In this regard, we start encoding variables where particle position is defined according to the below-mentioned vector:

$$x^{\left(i\right)}=\left(\begin{array}{l}{x}_{{}_1}^{\left(i\right)}\\ x_2^{\left(i\right)}\end{array}\right); i=1,\dots ,N_p$$

(42)

Initially, we consider that the LF has quadratic nature.

Consider

$$V\left(x^{\left(i\right)}\right)=2.25x_1^{2\left(i\right)}+x_1^{\left(i\right)}{x}_2^{\left(i\right)}+2x_2^{2\left(i\right)}; i=1,\dots , N_p$$

(43)

whose time derivate is specified below.

$$\dot{V}\left(x^{\left(i\right)}\right)=4x_1^{\left(i\right)}{x}_2^{\left(i\right)}-x_2^{2^{\left(i\right)}}-x_1^{\left(i\right)}\sin \left(x_1^{\left(i\right)}\right)-4x_2^{\left(i\right)}\sin \left(x_1^{\left(i\right)}\right); i=1,\dots ,N_p$$

(44)

Fitness value can be determined using maximum $\dot{V}\left(x\right)=4x_1{x}_2-x_2^2-x_1\sin \left(x_1\right)-4x_2\sin \left(x_1\right)$ for which the algorithm has an acceptable solution, as mentioned in the pseudocode specified in Section 3.

Subsequently, the proposed approach is used to ascertain $V\left(x^{\left(i\right)}\right)$, $\dot{V}\left(x^{\left(i\right)}\right)$ and $c$.

The settings used for the CKH technique are specified as follows. The scenario is based on a population size of 25, maximum diffusion speed of 0.005, foraging speed of 0.02, running for 150 iterations and a krill number of 25. The CKH algorithm was executed 25 times, and evolution was recorded in Figure 5a–d. These figures depict how the objective function evolved. After evaluating the optimisation outcome, the following was selected:

$$x^{\left(i\right)}=\left(\begin{array}{l}{x}_1^{\left(i\right)}\\ x_2^{\left(i\right)}\end{array}\right)=\left(\begin{array}{l}1.746\\ 0.747\end{array}\right)$$

(45)

while the domain is:

$$\mathrm{\Omega }=\left\{x\in R^n:V\left(x\right)=2.25x_1^2+x_1{x}_2+2x_2^2=9.2870\right\}$$

(46)

The dynamic of the state variables initialized from the tangency point state locus is represented in Figure 6. This is just to confirm the reliability of the obtained solution; indeed, the state variables converges asymptotically to the stable equilibrium point.

4.2. Example 2

The below-mentioned expression is the state-space equivalent of the Van der Pool equation [28].

$$\left\{\begin{array}{l}{\dot{x}}_1=-x_2\\ {\dot{x}}_2=x_1-x_2+x_1^2{x}_2\end{array}\right.$$

(47)

Here, $x={\left[x_1, x_2\right]}^T$ denotes the state vector. We use our approach comprising several chaotic maps for estimating the general DA corresponding to a stable equilibrium $x=(0,0)$. Initially, the dynamic model of Van der Pol is linearised in proximity to the equilibrium to determine a potential LF.

$$A={\frac{\partial f}{\partial x}‖}_{(0,0)}=\left[\begin{array}{cc}0& -1\\ 1& -1\end{array}\right]$$

(48)

The symmetric matrix $P$ is positive definite, and it is computed by solving the Lyapunov expression:

$$A^{T}P+PA=-Q$$

(49)

Q is an identity matrix.

$$Q=I$$

(50)

In this case,

$$P=\left[\begin{array}{cc}1.5& -0.5\\ -0.5& 1\end{array}\right]$$

(51)

The potential LF is specified as:

$$V\left(x\right)=x^T\left[\begin{array}{cc}1.5& -0.5\\ -0.5& 1\end{array}\right]x$$

(52)

which can be expressed as:

$$V\left(x\right)=1.5x_1^2-x_1{x}_2+x_2^2$$

(53)

The stable region corresponding to the nonlinear model is assured if the derivate of $V\left(x\right)$

$$\dot{V}\left(x\right)={\dot{x}}^{T}Px+x^{T}P\dot{x}$$

(54)

is negative definite.

$$\dot{V}\left(x\right)=-x_1^2-x_2^2-x_1^3{x}_2+2x_1^2{x}_2^2$$

(55)

In this regard, the following vector is initially employed for variable encoding related to particle position.

$$x^{\left(i\right)}=\left(\begin{array}{l}{x}_1^{\left(i\right)}\\ x_2^{\left(i\right)}\end{array}\right); i=1,\dots ,N_p$$

(56)

We initially assume that the LF is quadratic.

Let

$$V\left(x^{\left(i\right)}\right)=1.5x_1^{2\left(i\right)}-x_1^{\left(i\right)}{x}_2^{\left(i\right)}+x_2^{2\left(i\right)}; i=1,..., N_p$$

(57)

Then its time derivative can be specified as:

$$\dot{V}\left(x^{\left(i\right)}\right)=-x_1^{2\left(i\right)}-x_2^{2\left(i\right)}-x_1^{3\left(i\right)}{x}_2^{\left(i\right)}+2x_1^{2\left(i\right)}{x}_2^{2\left(i\right)}; i=1,...,N_p$$

(58)

Fitness value can be specified using the maximum $c^{*}$ where a feasible solution exists for the algorithm specified as pseudocode in Section 3.

Subsequently, use the proposed approach to determine $V\left(x^{\left(i\right)}\right)$, $\dot{V}\left(x^{\left(i\right)}\right)$ and $c$.

The CKH algorithm was configured as: 0.02 foraging speed, 150 iteration count, 0.005 maximum diffusion speed, 25 population size. The CKH algorithm was executed 20 times. Objective function evolution history is specified in Figure 7a and the DA show in Figure 7b. Considering the optimisation outcomes, we identify:

$$x^{\left(i\right)}=\left(\begin{array}{l}{x}_1^{\left(i\right)}\\ x_2^{\left(i\right)}\end{array}\right)=\left(\begin{array}{l}-0.8538\\ 0.7534\end{array}\right)$$

(59)

while the more extensive domain is:

$$\mathrm{\Omega }=\left\{x\in R^n:V\left(x\right)=1.5x_1^2-x_1{x}_2+x_2^2=2.3045\right\}$$

(60)

The dynamic of the state variables initialized from the tangency point state locus is represented in Figure 8. This is just to confirm the reliability of the obtained solution; indeed, the state variables converges asymptotically to the stable equilibrium point.

Unlike existing techniques, the developed method can be applied for both polynomial and rational LF. Indeed, it allows the development of a chaotic operative optimization algorithm which ensures the convergence to an enlarged DA in a significantly limited running time. Despite the performance, offering a larger accurate DA at high speed converging time, CKH is valuable for real-time trajectory tracking control problems

5. Discussion

This section comprises an in-depth discussion of a comparative analysis and assessment of the formulated CKH technique and peer-reviewed methods [28,29].

Table 1. Estimated DA features for CKH, ref. [28,29] methods.

Example	Systems Dynamics	Lyapunov Function	Estimated DA			Optimization Algorithm Convergence Time (ms)			Tangency Point (Best Solution)
			[29]	[28]	CKH	[29]	[28]	CKH	[29]	[28]	CKH
E1	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=-2x_1+x_1{x}_2\\ {\stackrel{·}{x}}_2=-x_2+x_1{x}_2\end{array}\right.$	$V\left(x\right)=x_1^2+x_2^2$	4.0804	4.112	4.0955	-	6.6	0.350	-	-	$\left\{\begin{array}{l}{x}_1=1.223\\ x_2=1.611\end{array}\right.$
E2	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=-x_2\\ {\stackrel{·}{x}}_2=x_1-x_2+x_1^2{x}_2\end{array}\right.$	$V\left(x\right)=1.5x_1^2-x_1{x}_2+x_2^2$	2.09	2.318	2.3045	-	6.7	0.542	-	-	$\left\{\begin{array}{l}{x}_1=0.852\\ x_2=-0.755\end{array}\right.$
E3	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=x_2\\ {\stackrel{·}{x}}_2=-0.2x_2+0.81\sin \left(x_1\right)\cos \left(x_1\right)-\sin \left(x_1\right)\end{array}\right.$	$V\left(x\right)=x_1^2+x_1{x}_2+4x_2^2$	0.699	0.708	0.6993	-	7.2	0.236	-	-	$\left\{\begin{array}{l}{x}_1=0.741\\ x_2=-0.307\end{array}\right.$
E4	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=-0.25x_1+\ln \left(1+x_2\right)\\ {\stackrel{·}{x}}_2=-0.3750x_1-0.2x_1{x}_2+\left(0.125x_1-x_2\right)\cos \left(x_1\right)\end{array}\right.$	$V\left(x\right)=x_1^2+x_2^2$	0.2737	0.278	0.2737	-	8.3	0.258	-	-	$\left\{\begin{array}{l}{x}_1=0.443\\ x_2=-0.277\end{array}\right.$
E5	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=-x_1+x_2{x}_3^2\\ {\stackrel{·}{x}}_2=-x_2+x_1{x}_2\\ {\stackrel{·}{x}}_3=-x_3\end{array}\right.$	$V\left(x\right)=x_1^2+x_2^2+x_3^2$	4.9188	4.971	4.969	-	8.4	0.77	-	-	$\left\{\begin{array}{l}{x}_1=1.200\\ x_2=1.452\\ x_3=-1.190\end{array}\right.$
E6	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=1+x_3+0.125x_3^2-{\mathrm{e}}^{\left(x_1\right)}\\ {\stackrel{·}{x}}_2=-x_2-x_3\\ {\stackrel{·}{x}}_3=-x_2-2x_3-0.5x_1^2\end{array}\right.$	$V\left(x\right)=x_1^2+x_2^2+x_3^2$	2.655	2.887	2.6617	-	8.6	0.813	-	-	$\left\{\begin{array}{l}{x}_1=-1.325\\ x_2=0.611\\ x_3=-0.728\end{array}\right.$

lists the approximate DA features corresponding to six dynamic nonlinear benchmark models having a quadratic LF, computed using the literature [28,29]. The nonlinear polynomial and nonpolynomial variants of the investigated examples are presented. E1, E2, E3 and E4 are second-order systems. Nevertheless, E5 and E6 are nonlinear polynomial and nonpolynomial third-order models. The DA remains the primary assessment criteria. The observed outcomes are better, concerning DA values. Figure 9 shows the approximated DAs for the origins in examples E1–E6 described in Table 1 using the CKH method.

Furthermore, the formulated technique offers a comprehensive solution that comprises the definition of a LF, concluding with an enhanced DA. This outcome is obtained for every example assessed in this study. Implementation was straightforward and less time-intensive.

Every example uses the CKH technique to maximise the permissible domain value. These values are contrasted against the results identified using a peer optimisation-specific method mentioned in the literature [28,29]. It is evident from Table 1 that the six assessed samples yield better domains of attraction when identified using the CKH method as compared to optimisation-based techniques. Formulation (11) comprises nonlinear optimisation for a problem with numerous solutions having improper DA estimations corresponding to the system being studied. For instance, consider the many solutions of problem (11) corresponding to systems in E1 and E2, as indicated in Figure 10 and Figure 11.

Figure 10 and Figure 11 depict LF level sets that are indicated using a solid line. Dashed lines are used to depict the time derivative level set at level zero. It is also noteworthy that, like the previous instance, E1 has a unique solution corresponding to $V\left(x^{\left(i\right)}\right)=c$. The discussion ahead specifies a unique solution which does not satisfy the tangency constraint for the method described in [28]. The later weakness is well clarified in Figure 10a where a clear intersection between the LF and its derivative is shown. Consequently the value of solution $V\left(x^{\left(i\right)}\right)=c=4.112$ provided in Table 1 is not realistic. The solution corresponding to E2 solves the problem (11), where c is less than the required estimate in Figure 11a. This solution correlates to a transversal intersection between level sets $V\left(x^{\left(i\right)}\right)=c$ and $\stackrel{·}{V}\left(x^{\left(i\right)}\right)=0$. It may be observed that these level sets intersect at a significant distance from the origin. Hence, a minor fraction of the level set $V\left(x^{\left(i\right)}\right)$ provides an actual DA estimate (the small circle close to the origin).

For the Van der Pool model, all techniques specified in the literature provide different results. Using the proposed method, we found $V\left(x\right)=2.3045$, while for the technique in [28], the corresponding value is 2.318, which corresponds to an explicit intersection between the LF and its derivative and should be rejected as a solution. At the same time, Hachiho’s technique yields a value of 2.09. All these methods were assessed using identical LF: $V\left(x\right)=1.5x_1^2-x_1{x}_2+x_2^2$. Figure 11 depicts the numerous solutions corresponding to the system (A1). The solid line indicates the distinct levels sets corresponding to the LFs. The time derivative for zero value is also depicted using dashed lines.

The solution presented in Figure 11a addresses the problem (11), with c more than the required estimate. When the solution point is determined, two solution points are observed specific to the transversal intersection of level sets described by $V\left(x\right)=2.318$ and $\stackrel{·}{V}\left(x\right)=0$. Nevertheless, neither of these two are precise DA estimates. Considering that all solutions depicted in Figure 11a are applicable to the problem (11), but none are appropriate, these solutions are therefore considered “dummy”.

In Figure 11c, we have indicated that $V\left(x\right)=2.09$ and level set $\stackrel{·}{V}\left(x\right)=0$ do not intersect. The solution presented in the figure corresponds to the problem (11) having a c value less than the required DA estimation. Nevertheless, DA estimation size is insufficient since the LF level set and that corresponding to its derivative do not intersect.

Considering Figure 11b, it can be observed that the LF intersects its derivative tangentially at two coordinates owing to symmetry. One result is reported because the two solutions are similar concerning DA estimation.

An important observation related to the time cost shown by the compared methods involves the superiority of the CKH hybrid method developed in this paper. Indeed, in [29] no further data have been reported regarding the time cost. In [28] authors reported more than 5 ms time cost overall for the studied examples. It is interesting to notice that CKH involves less than 0.5 ms as a consuming time in the overall investigated example. This is to emphasize the fact that CKH is the better and more appropriate technique to be used in online trajectory tracking problems.

Table 2. Summary comparative key performance of the investigated methods [28,29] in contrast with the designed CKH method.

Method	Pros	Cons
CKH	- High accuracy in defining the solution characterizing the tangency point coordinates and the DA radius. - Reliable meta-heuristic technique, high convergence speed due to the induced chaotic movement parameter with strong robustness. - Global enhanced search ability. - Local optimum avoidance. - Dummy DA solution avoidance. - No random space research is needed. - Interesting feasibility for a wider range of benchmark models. - Performs easily and simply, and no need of parameter tuning. - Efficient for both polynomial and rational function - Reasonable amount of computational requirements allowing to solve control tracking problems. - Smaller number of parameters	- Suitable chaotic map should be selected. - Suitable algorithm parameter values should be initialized.
Sampling technique	- Perform easily and simply. - Efficient for both polynomial and rational Lyapunov function.	- May diverge - Needs intensive numerical simulationto fix the solution. - Risk of stagnation in local optimum.
Analytical method	- A step by step recursive algorithm is well defined analytically. - No dummy solution can be found.	- May diverge - Delicate algebraic and analytical developments are indispensable. - Risk of stagnation in local optimum - Polynomial function approximation is needed.

summarises qualitatively the analytical technique given in [29], the sampling method given in [28] and the CKH optimizing method developed in this work. The comparison criterion focuses mainly on the ability to implement the related algorithms easily and to avoid local optimizing solutions, while ensuring convergence [24].

An advanced metaheuristic strategy pointing to maximizing the DA of stable equilibriums using Lyapunov theory is established in this work. The estimated region is fully bounded by the domain of the LF negative definiteness and its time derivative. The CKH optimizer is designed such that the requirements on the sign of the LF and the tangency constraint between the level sets are satisfied. As a result, numerous potential local solutions have been avoided by such constraints.

Unlike existing techniques, the developed method can be applied for both polynomial and rational LF. Indeed, it allows the development of a chaotic operative optimization algorithm which ensures the convergence to an enlarged DA in a significantly limited running time. Despite the performance, offering a larger accurate DA at high speed converging time, CKH is valuable for real-time trajectory tracking control problems.

To highlight the performance of the CKH technique a comparative analysis is performed with a swarm-intelligence based approach. The study compares the CKH with the CPSO technique while considering the performance criteria as the DA radius, the running time (ms) and the capability to provide an accurate description of the tangency point. The result of this study is summarized in

Table 3. Estimated DA features for CKH and CPSO methods.

Example	Systems Dynamic	Lyapunov Function	Domain Radius		Running Time (ms)		Tangency Point
			CPSO	CKH	CPSO	CKH	CPSO	CKH
E2 [28,29]	$\left\{\begin{array}{l}{\stackrel{·}{x}}_1=-x_2\\ {\stackrel{·}{x}}_2=x_1-x_2+x_1^2{x}_2\end{array}\right.$	$V\left(x\right)=1.5x_1^2-0.5x_1{x}_2+x_2^2$	2.3045	2.3045	1.821	0.542	$\left\{\begin{array}{l}{x}_1=-0.8575\\ x_2=0.7483\end{array}\right.$	$\left\{\begin{array}{l}{x}_1=0.852\\ x_2=-0.755\end{array}\right.$

.

It appears in this table that both methods are highly comparable in term of provided performance. However, the CKH showed a superiority in the operating consumed time. This concludes that the CKH is an efficient, fast and accurate optimization technique.

6. Conclusions

The CKH algorithm provided the largest solution for the DA optimization. It was implemented in this work to compute the QLF, and allows expanding of the DA for a large class of nonlinear systems. The main objective of computing the largest sublevel set of the DA was perfectly achieved. The CKH-based iterative technique was developed to calculate DA lower bound by exploiting CKH capacity and implementing a novel enhanced KH. The synthesized design establishes chaotic maps for global optimization analysis and encodes the systems’ states that will be determined as particle positions.

The observations and conclusions of this study are summarized below:

• The CKH technique suggested in the present study is straightforward, flexible, and easily understood. There are no challenges concerning the timing of parameter tuning; hence, it may be implemented for complex computational optimizations.
• The outcomes of the CKH algorithm, when benchmarked against the results of other popular optimization techniques, indicate the flexibility of the suggested approach.
• The proposed CKH technique requires less CPU time compared to other techniques because of the integration of chaos maps that provide chaotic maps benefits.

Hence, the benefits associated with the proposed technique used for estimating the DA showed that the concurrent use of KH with chaotic maps leads to rapid convergence and provides superior performance and better results. The chaotic maps produce a search space different from that produced by the KH process.

Motivated by the efficient results attained for the studied benchmark examples, the developed approach can be applied for real time tracking control problems with delicate time budget constraints. Theoretical concerns related to rational LF are a challenging aspect and a key perspective promoting this work.