A Hybrid Approach for Solving Systems of Nonlinear Equations Using Harris Hawks Optimization and Newton’s Method

Systems of nonlinear equations are known as the basis for many models of engineering and data science, and their accurate solutions are very critical in achieving progress in these fields. However, solving a system with multiple nonlinear equations, usually, is not an easy task. Consequently, finding a robust and accurate solution can be a very challenging problem in complex systems. In this work, a novel hybrid method namely Newton-Harris hawks optimization (NHHO) for solving systems of nonlinear equations is proposed. The proposed NHHO combines Newton’s method, with a second-order convergence where the correct digits roughly double in every step, and the Harris hawks optimization (HHO) to enhance the search mechanism, avoid local optima, improve convergence speed, and find more accurate solutions. We tested a group of six well-known benchmark systems of nonlinear equations to evaluate the efficiency of NHHO. Further, comparisons between NHHO and other optimization algorithms, including the original HHO algorithm, Particle Swarm Optimization (PSO), Ant Lion Optimizer (ALO), Butterfly Optimization Algorithm (BOA), and Equilibrium Optimization (EO) were performed. The norm of the equation system was calculated as a fitness function to measure the optimization algorithms’ performance. A solution with less fitness value is considered a better solution. Furthermore, the experimental results confirmed the superiority of NHHO over the other optimization algorithms, in the comparisons, in different aspects, including best solution, average fitness value, and convergence speed. Accordingly, the proposed NHHO is powerful and more effective in all benchmark problems in solving systems of nonlinear equations compared to the other optimization algorithms. Finally, NHHO overcomes the limitations of Newton’s method, including selecting the initial point and divergence problems.


I. INTRODUCTION
Systems of nonlinear equations are the basis of many fields such as engineering models and data sciences. However, finding accurate solutions to systems of nonlinear equations is a challenging issue. Let F(x) be a nonlinear multivariable The associate editor coordinating the review of this manuscript and approving it for publication was Easter Selvan Suviseshamuthu . function such that F(x) : R n − > R n . Finding the solution α = (α 1 , α 2 , . . . , α n ) t of F (X ) = 0 is one of the prevalent problems in applied sciences. This type of equation can occur directly in many sciences including biology, physics, chemistry, and many engineering branches. Besides, systems of nonlinear equations can occur indirectly in many nonlinear problems that come in the form of differential equations or integral equations. These problems can be converted into VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ systems of nonlinear equations using some numerical techniques such as the finite difference method [1]. Furthermore, finding the exact solution to nonlinear problems is an intricate process. For instance, there is no algebraic formula for the general solution of a quintic polynomial (Abel's theorem) [2]. Many equations have higher orders and may include specific hyperbolic, exponential, and trigonometric terms, or they may be entirely transcendental equations that are difficult to be solved analytically to achieve the exact solution. As a consequence, the field of finding approximate solutions to nonlinear problems has become very common. Population-based approaches (PBAs) have been used as alternatives for finding the solution(s) of nonlinear systems of equations. We can define PBAs simply as a group activity that has arisen from social insects operating under very few laws. Usually, the popular algorithms of the population-based approaches result from animals' lifestyles. Many population-based approaches have been proposed to solve optimization problems in the literature, for instance, Harris hawks optimization (HHO) [3], equilibrium optimization (EO) [4], particle swarm optimization (PSO) [5]- [7], Slime Mould algorithm (SMA) [8], ant lion optimization (ALO) [9], grasshopper optimization algorithm (GOA) [10], genetic algorithm (GA) [11], firefly algorithm (FA) [12], Cuckoo search algorithm (CSA) [13], artificial bee colony (ABC) [14], krill herd algorithm (KHA) [15], ant colony optimization (ACO) [16], and glowworm swarm optimization (GSO) [17]. Many other examples of optimization algorithms can be found in [18].
Nonlinear equation systems have been solved using several optimization algorithms such as FA, CSA, GA, ABC, and PSO. To boost the canonical optimization algorithms performance, some researchers modified the optimization algorithm to better suit the targeted problem for more accurate solutions. For example, Ariyaratne et al. [19] presented a modified FA to solve nonlinear systems of equations in an optimization problem, such that the continuity, differentiation, and initial assumptions were not required, and it can provide multiple root approximations simultaneously. Another approach was proposed by Zhou and Li [20] to solve systems of nonlinear equations by presenting a modified CSA variant that provides a harmonious solution without modifying the cuckoo's evolution. Similarly, Chang [21] solved systems of nonlinear equations using a suitable modification to the GA to estimate the parameters of the problems.
By presenting the systems as a multi-objective optimization problem, and then solving it using GA, Grosan and Abraham [22] solved the complex systems of nonlinear equations. In contrast, Ren et al. [23] used efficient GA with harmonious and symmetric individuals to solve nonlinear systems of equations. Also, Jaberipour et al. [24] solved systems of nonlinear equations using a modified PSO algorithm. The modification aimed at overcoming the fundamental PSO drawbacks such as trapping at local minimums and the slow convergence. Then as well, Mo and Liu [25] solved a system of nonlinear equations by an improved PSO algorithm where the improvement incorporated the conjugate direction method (CDM) into the original algorithm. The CDM increased the algorithm efficiency in optimizing high-dimensional problems and overcame the local minima incident [26].
Other researchers went further in searching for better solution approximations to the nonlinear systems of equations by combining two PBAs. As a result, the hybrid will inherit the advantages and alleviates the weaknesses of both algorithms [27]. Some examples of hybridizing PBAs are the hybrid of the ABC and PSO [28], hybrid GA [29], [30], hybrid FA [31], hybrid PSO [32], hybrid KHA [33], hybrid ABC [34], and many others [36]- [40].
Mathematicians usually employ mathematical methods that depend on successive approximations, called iterative methods, to find an approximate solution to nonlinear systems of equations. The most widespread iterative method is the one proposed by Isaak Newton, known as Newton's method or Newton-Raphson method. Newton's method has been widely studied and modified to increase both the accuracy of the approximate solutions and the convergence rate to the desired solution [40]- [42] Most of the proposed optimization algorithms in the literature deal with the systems of nonlinear problems using a single optimization algorithm, or by using a hybridized algorithm of two different algorithms. In the proposed work, we the systems of nonlinear problems using a combination of the well-known HHO optimizer and Newton's iterative technique to benefit from the fast convergence of Newton's method towards the desired solution. However, few types of research were proposed based on the idea of hybridizing Newton's method with an optimization algorithm. Karr et al. [43] presented a hybrid technique that combines GA and Newton's method for solving systems of nonlinear equations. The GA was used to find proper initial solutions which will be later fed to Newton's method. Besides, another hybrid algorithm that combines the GA to Powell algorithm together with Newton's method was proposed by Luo et al. [44] for solving systems of nonlinear equations. Finally, Luo et al. [45] presented a hybrid scheme for solving nonlinear systems using the quasi-Newton method and chaos optimization. However, all of these works have focused on a specific issue or a special type of problem and did not apply their proposed algorithms to arbitrary nonlinear systems of equations.
It is worth mentioning that optimization algorithms model natural phenomena where the local search mechanism in an optimization algorithm is critical to obtain better solutions. However, there is no guarantee that the search always converges to the optimal solution in the search space especially if the search space contains many candidate solutions. Further, in the case of solving nonlinear systems, optimization algorithms attempt to avoid premature convergence that can lead to local optima [18]. However, obtaining more accurate solutions in less time, or iterations, is considered an advantage for the optimization algorithm in solving nonlinear equations.
As mentioned earlier, hybridized algorithms are created by combining two optimization algorithms to overcome the demerits of one algorithm and taking advantage of the merits of the other. However, most of the hybridized algorithms still suffer from premature convergence depending on the original optimization algorithms mechanism [46]. The key point is to find a good match between the selected algorithms while creating the hybrid to obtain a more efficient variant. The following paragraph highlights the selected hybrid poles, The HHO and Newton's method, coupled with reasons for this candidature.
On one hand, Harris Hawk optimization (HHO) is a fast, powerful, and efficient PBA. The algorithm imitates the style of Harris Hawk birds in looking and pursuing the prey in nature. A rabbit symbolizes the prey in the algorithm, which is treated as the best solution [46]. HHO algorithm has gone beyond other well-known algorithms, such as GOA, PSO, EO, ALO, and BOA [3]. However, HHO has some constraints as most optimization algorithms. Firstly, Harris hawks' solutions may end up at local optima instead of the optimal solution. Besides, to specify the type of search that the algorithm will follow, HHO depends on the rabbit energy. The energy of the rabbit begins at approximately 2 and progressively decreases to 0 through iterations. Therefore, when the energy of the rabbit is greater than 1 (E ≥ 1) in the first half of iterations, the global search is performed. However, when the energy of the rabbit is less than 1 in the second half of iterations, HHO does not conduct a global search, although the area might not be the optimal one [47]. Hence, instead of converging to a global solution, HHO may converge to a local solution which causes a premature convergence [48].
Several researchers have modified the original HHO to accelerate the global search phase and to prevent it from falling into local search space. Dhawale and Kamboj [49] proposes a hybrid meta-heuristic optimization technique that combines the benefits of HHO and Improved Grey wolf optimization (IGWO) to develop a more efficient method for reaching an exploration to the exploitation phase. In addition, Kamboj et al. [50] developed a hybrid variant of Harris Hawks optimizer using Sine-Cosine algorithm. Further, an improved version of the HHO algorithm that combines HHO with Canis lupus inspire grey wolf optimizer (hHHO-GWO) has been proposed by Nandi and Kamboj [51] to find the solution to various optimization problems. Recently, Krishna et al [52] proposed a hybrid Harris hawks pattern search algorithm (hHHO-PS). In their next work, Nandi and Kamboj [53] hybridized the HHO with sine-cosine algorithm (SCA) using memetic algorithm approach to solve the photovoltaic constrained UCP of an electric power system. Other works have been proposed by improving different optimization algorithms, including SMA algorithm, to solve different numerical and engineering design challenges such as in [54] and [55].
On the other hand, Newton's method is limited to initial point selection and divergence problems caused by improper guesses of the initial points. In other words, it performs a narrow and local search that entirely depends on the initial given values as a start with the incapability of performing a global search.
The aforementioned has motivated the proposal of an HHO hybrid for solving nonlinear equations systems. To improve the search mechanism of the HHO, Newton's method is applied within the algorithm as a local search assistant. Newton's iterations enhance the distribution of the initial solutions in the search space. In comparison to the random distribution used for the original HHO, Newton's technique improves the algorithm's computational accuracy and speeds up its convergence rate.
Therefore, this work is intended to improve the accuracy of the solution of nonlinear equations systems by combining HHO Newton's method, as a search assistant, and HHO. The key contributions to this paper are summarized as follows: 1) To propose an improved hybrid HHO variant, namely Newton Harris hawks optimization (NHHO), that relies on Newton's method to improve the local search in HHO and accelerate its convergence rate. 2) To apply the proposed NHHO algorithm for solving nonlinear systems of equations. 3) To evaluate the efficiency of the proposed NHHO with several well-known optimization algorithms, including HHO, PSO, EO, ALO, and BOA in terms of accuracy, fitness value, and satisfying the convergence criterion. The rest of the work is arranged as follows: The stateof-the-art Newton's method and HHO algorithm are introduced in Section II. Section III describes the details of the proposed NHHO. The details of the experiments, the benchmark systems, and the results are presented in Section IV. A discussion and analysis of the results are presented in Section V. Eventually, Section VI states the conclusion of the paper.

II. STATE-OF-THE-ART A. NEWTON'S METHOD
The method which aims to find an approximate solution using successive approximations is called an iterative method. In general, iterative methods do not attain the exact solutions; because of that, researchers usually select a tolerance limit for the solutions obtained by the iterative methods, which is the difference between the exact and the approximate solutions. The most widespread iterative method is the one proposed by Isaak Newton, known as Newton's method or Newton-Raphson method. The method is given by where F (X n ) is the Jacobian of F(X ). Newton's method is of the second order of convergence and can be applied easily to a wide range of nonlinear algebraic problems. For this reason, mathematical programs such as Mathematica and MATLAB use built-in functions based on Newton's method to find the roots of nonlinear equations. The principle of Newton's method is that continuous and differentiable functions can VOLUME 9, 2021 be approximated by a straight tangent line to them. That is, starting by initial solutionx 0 , the approximation x 1 is the x-intercept of the tangent line to the graph of f at (x 0 , f (x 0 )).
In the same manner, we can find the next approximation x 2 by finding the x-intercept of the tangent line to the graph of (x 1 , f (x 1 )) and so on, as illustrated in Figure 1 [56]. Newton's method has been widely studied and modified to increase the accuracy of the approximate solutions of the nonlinear problems and increase the order of convergence that reflects consecutively on the speed of obtaining the desired solution [40]- [42].

B. HARRIS HAWKS OPTIMIZATION (HHO)
Harris hawks optimization (HHO) is a recent optimization algorithm introduced by Heidari et al. in 2019. HHO is considered as a meta-heuristic, population-based, and nature-inspired algorithm that imitates the Harris hawks birds' style in chasing rabbits. The hawks' natural strategies are applied in the algorithm, such as exploring the area looking for prey, preaching, surprise pounce, and predation strategies. As meta-heuristic algorithms, HHO is composed of exploration and exploitation phases for searching on the prey. More precisely, HHO consists of two exploration phases that make the hawks explore different areas in search of the rabbit and four stages of exploitation in which the hawks search for the rabbit in the neighborhood, as shown in Figure 2.

1) EXPLORATION PHASE
In the exploration phase, the Harris hawks search for the rabbit, representing the best location in the algorithm. The Harris hawks have sharp eyes that help them to monitor the area to find the rabbit. Practically, each Harris hawk is a candidate solution, and its fitness value is calculated based on its location from the rabbit. However, it is not easy sometimes to detect the prey. In this case, the Harris hawks wait on a high tree for a while hoping to observe the prey based on the following equation: where X (t + 1) represents the position of the Harris hawk in the next iterationt, X rabbit (t) is the location of the rabbit, and X (t) is the current location of the hawk. The average of the current Harris hawks locations is represented by X m (t).
While the variables r 1 , r 2 , r 3 , r 4 , and q are random variables between [0, 1], LB and UB refer to the lower and upper limits of the variables. To calculate the average locations of the Harris hawks X m (t), HHO uses a simple way as in the following equation: where N represents the maximum number of Harris Hawks (solutions), and X i (t) is the location of the Harris hawks in iteration t.

2) THE TRANSITION FROM EXPLORATION TO EXPLOITATION
Indicating the right balance to shift from the exploration and exploitation phases is critical to the algorithm's performance. HHO assumes that the rabbit's escaping energy (E) is reduced while escaping from the hawks. Therefore, it relies on the value of (E) to perform the transitions between the two phases as shown in the following equation: where E refers to the rabbit's escaping energy, E 0 is random initial energy between (−1, 1), and T is the maximum number of iterations. In HHO, the exploration phase occurs when the value of |E| ≥ 1, and the exploitation phase when |E| < 1, as shown in Figure 3.

3) EXPLOITATION PHASE
In the exploitation phase, the hawks identify the rabbit location in the previous stage and are ready to attack. However, the rabbit always attempts to escape from the attack. Therefore, the four exploitation phases are designed based on the absolute value of rabbit's escaping energy |E| and the probability of escaping (r) as illustrated in Figure 3. A summary of each attack strategy is introduced in the following.

4) SOFT BESIEGE
In the soft besiege, the rabbit has some energy with a high escaping probability as |E| ≥ 0.5 and |r| ≥ 0.5. In this case, the hawk softly encircles around the hawks as in the following equation: where X (t) refers to the difference between the rabbit location and hawk locations in iteration t, and J is the strength of the random rabbit jump with J = 2(1−r 5 ) as r 5 is a random variable in the interval [0, 1].

5) HARD BESIEGE
In the hard besiege, the rabbit has lost most of its energy, but it still has some escaping probability as |E| < 0.5 and |r| ≥ 0.5. Consequently, the Harris hawk hardly encircles around the rabbit to prepare for the final surprise pounce using the following mathematical equation:

6) SOFT BESIEGE WITH PROGRESSIVE RAPID DIVES
In this situation with |E| ≥ 0.5 and |r| < 0.5, the rabbit can still escape an attack. Accordingly, the hawk performs some dives and smart zigzag movements around the rabbit. First, the hawk tries to approach the rabbit based on the following equation: Second, a decision is made based on the last dive and the expected result to dive again or not. If not, the hawk begins to perform smart and irregular dive based on Levy flight concept (LF) as in the following equation: where dim refers to the dimension of the optimization problem and S is a vector randomly generated of size 1×dim. The LF function is calculated as: where u and v are random variables between [0, 1], and β = 1.5 is a constant. Therefore, updating the locations of Harris hawks in soft besiege can be formulated as follows: where F represents the calculated fitness value for Y and Z .

7) HARD BESIEGE PROGRESSIVE RAPID DIVES
In the case of |E| < 0.5 and r < 0.5, the rabbit has low energy and less chance to escape. The hawk attempt to get closer to the rabbit by performing rapid dives before performing the surprise pounce as formulated in the following equation: where Finally, the fitness function value is calculated based on the Euclidean norm, which is also called norm-2 or square norm. For this research, the solution with less norm is more accurate than a higher norm solution. Hence, a solution with norm= 0 is a precise solution. This norm represents the ordinary distance from the origin to the vector F (x) = (f 1 , f 2 , . . . , f n ) and can be defined as: Aside from all of the advantages, the HHO, like any other stochastic population-based optimizer, has some drawbacks in terms of its original ability. For instance, When the HHO population becomes stuck in local optima for a complex problem, a proper modification is needed to speed up the algorithm's exploration. Another shortcoming is that HHO does not always maintain the proper balance between the cores, or it does not always find genuine moments, or the transition is not always smooth, particularly when dealing with very complex attribute spaces. Most of the HHO-based research works have been pointed to these two reasons to improve the convergence speed or local optima avoidance [57].

III. NEWTON-HARRIS HAWKS ALGORITHM (NHHO)
The original HHO is a robust optimization algorithm that can be used to solve many problems. However, based on the no free lunch theorem (NFL) [58], no algorithm is perfect for solving all problems. The proposed NHHO aims to improve the performance of HHO in solving systems of nonlinear equations by adopting Newton's method. Newton's method in NHHO works as a local search to enhance the search mechanism in HHO as illustrated in Figure 4.
Newton's method is applied to the rabbit location in every iteration and then the fitness value of the rabbit location is compared to the fitness of the location calculated by Newton's method, which is represented by (X n+1 ) in Figure 4. Consequently, the rabbit location equals the location with a minimum fitness value. The pseudocode of the proposed NHHO algorithm is presented in Algorithm 1.
The algorithm shows the different stages of the HHO algorithm, including initialization, exploration, and exploitation stages. The transition from exploration to exploitation depends on the rabbit energy |E|. Exploitation is performed when the value of |E| ≥ 1 and one of the four exploitation phases is selected when |E| < 1. Newton's method, in the red box, is applied at the end of each iteration. The fitness value is calculated for the new location provided by Newton's method and compared to the X rabbit location provided by HHO. Finally, the location that has better fitness is selected.

IV. EXPERIMENTS
To evaluate the proposed NHHO performance, comparisons between NHHO and other optimization algorithms on solving Algorithm 1 Newton Harris Hawks Optimization (NHHO) Input: N: population size, T: maximum iterations. Output: The rabbit location (X rabbit ), and the rabbit fitness Initialize population randomly X i (i = 1, 2, . . . , N ) While (less than maximum iteration (t < T )) do Evaluate upper and lower boundaries and the hawks' fitness Set the hawks best location to X rabbit For every hawk X i do Update initial rabbit energy E 0 Update rabbit energy E using Equation (4) If (|E| ≥ 1) then Update the hawks' location using Equation (2) If (|E| 1) then If |E| ≥ 0.5&&r ≥ 0.5 then %Soft besiege Update the hawks' location using Equation (5) ElseIf |E| < 0.5&&r ≥ 0.5 then %Hard besiege Update the hawks' location using Equation (7) ElseIf |E| ≥ 0.5&&r < 0.5 then %Soft besiege with progressive rapid dives Update the hawks' location using Equation (11) ElseIf |E| < 0.5&&r < 0.5 then %Hard besiege with progressive rapid dives Update the hawks' location using Equation (12) Calculate Newton's location X n+1 using Equation (1) Calculate the fitness of X n+1 and X rabbit using Equation (15) If Fitness (X n+1 ) < Fitness (X rabbit ) X rabbit = X n+1 Return the rabbit location (X rabbit ) six systems of nonlinear equations were performed. The six systems are general problems and known benchmarks that other researchers studied. Each problem was solved using NHHO and other well-known optimization algorithms including the original HHO, PSO, ALO, BOA, and EO algorithms.
Further, to assure fair comparisons to all optimization algorithms, the comparisons were made based on the fitness value of F(x) presented in Equation (15). The same settings have been used in all experiments, such as search agents (population size) was set to 20 and maximum iteration to 100. Furthermore, it is well-known that the performance of optimization algorithms is improved by fine parameter tuning. In this work, the parameter setting of the utilized algorithms is closed based on the original work as shown in Table 1.
The optimization algorithms were implemented using MATLAB software R2020a. The experiments were performed on an Intel i5 processor PC with 2.2 GHz, 8 GB of RAM, and Windows 8 for the operating system.
To perform the experiments six problems (benchmarks), which are commonly used in the literature, are selected. See for instance the following references [18], [22], [23], [55]. The selected problems are of different dimensions and contains several types of transcendental functions, for instance: polynomials and trigonometric functions.

A. PROBLEM 1
This problem is represented by a system of two nonlinear equations as following [23]: The best solution obtained by optimization algorithms after running the algorithms 30 times is shown in Table 2. The performance of the provided solution is measured by calculating its fitness value (F).
The proposed NHHO obtained the best solution with x = 2 and y = 1, which is one of the exact solutions for the equations in problem 1. The fitness value of NHHO (0 in the case of exact solution) outperformed all other algorithms. Running the experiments for 30 runs can, sometimes, give alternative solutions with the same fitness, such as x = 1 and y = 2 for NHHO, or with lower fitness values. However, for this research, we consider the best solution and the average fitness value of solutions. It is worth mentioning that none of the problems, except for problem 1, has two best solutions with the same fitness value.
Graphical representation can illustrate the solutions of nonlinear systems. Problems 1, 3, and 4 are of nonlinear systems of dimension two. Thus, we are able to illustrate the equations on the coordinate plane. However, this is not possible for the remaining problems as dimension n ≥ 3. Note that the solution(s) of the systems is the point(s) of the intersection between the two curves. The graphical representation for benchmark problem 1 is illustrated in Figure 5.

B. PROBLEM 2
In this problem, the system is of three nonlinear algebraic equations as following [23]: where The best solution obtained by the optimization algorithms for problem 2 is presented in Table 3. With fitness value (F) = 2.2204E −16, the proposed NHHO obtained the best solution compared to the other algorithms. The solutions in the tables are rounded to 12 digits. That is why NHHO and PSO have obtained the same solution in the table but with different fitness values. For instance, the NHHO solution rounded to 17 digits is as following x = −0.03275902177937206, y = 1.264628719362605, and z = 1.400643891245019 while the solution of PSO is x = −0.03275902177936415, y = 1.2646287193626, and z = 1.400643891245025.

C. PROBLEM 3
The system in this problem consists of two nonlinear algebraic equations; such that f is not differentiable as follows [23]:    Here, the best solution for problem 3 was found by the proposed NHHO and PSO with the same fitness value (F) = 3.3566E − 16. Besides, both of the optimization algorithms found the same solution with x = 1.33981159986596 and y = 2.832851967564935. In contrast, BOA has given an approximate solution, but not as accurate as NHHO and PSO. The result of the best solution attained by running each algorithm 30 times is shown in Table 4. The graphical representation for problem 3 is illustrated in Figure 6.

D. PROBLEM 4
In this benchmark problem, there are two nonlinear algebraic equations as following [62]: The best solution obtained by the optimization algorithms for problem 4 is presented in Table 5. With fitness value (F) = 1.1102E − 16 the proposed NHHO has achieved the best solution and outperformed the other optimization algorithms. The solutions of all algorithms, except for EO, are close to each other in the first few digits. However, the proposed NHHO has proved to be more accurate in solving nonlinear equations due to the powerful local search mechanism NHHO adopts in searching for the best solution in the search space. The graphical representation for problem 4 is illustrated in Figure 7.

E. PROBLEM 5
The system in this problem contains ten nonlinear algebraic equations [22]: Similarly, the proposed NHHO has been approved to be efficient in solving a complex set of nonlinear equations as shown in Table 6. NHHO algorithm outperformed all of the other compared optimization algorithms with a significant difference. While NHHO came first with a fitness value of 5.2147E − 17, PSO achieved a fitness value of 6.1466E − 06. The original HHO came last with fitness equals 2.7272E −01.
This system of nonlinear equations is a complex system that consists of ten equations and ten variables. Usually, such a system has many solutions. In other words, the possibility is higher for optimization algorithms to be trapped in local optima. There are three different solutions in Table 5 with different fitness values. However, the hybrid NHHO can improve its solution in every iteration as it has the power of both search mechanisms; the local search of the HHO and Newton's method.
The best solution obtained by the optimization algorithms for problem 6 is presented in Table 7. With fitness value (F) = 6.9871E − 21, the proposed NHHO obtained the best solution compared to the other algorithms. PSO and EO came second with fitness values equal to 1.5098E − 11 and 7.2487E − 10, respectively. From problem 5 and problem 6, it can be noticed that the proposed NHHO is flexible and can work efficiently even in a wide interval [−10, 10].
Eventually, the obtained results by the proposed hybrid NHHO on all the previous problems confirm the theoretical hypothesis discussed earlier, i.e., the hybrid will inherit the powerful merits from both algorithms. This is obvious in case the comparison has been made separately between NHHO and the HHO where the proposed hybrid surpasses the original HHO in all obtained fitness values. This is imputed to the employed local search (Newton's method) which enhanced the hybrid capability to avoid the local optima in Problem 1 VOLUME 9, 2021  (since the fitness is 0) and significantly improve the obtained fitness values of all other algorithms in the remaining problems (which will be considered as the global optima until better solutions could be found). Furthermore, the obtained standard deviation values by the proposed hybrid indicate more stable and repeatable results in the majority of the problems. The aforementioned observations are considered as a piece of strong evidence that the proposed method has successfully avoided the entrapment at least in the local optima that all other algorithms have been in the tested problems.

V. RESULTS AND ANALYSIS
To analyze the efficiency and the repeatability of the proposed NHHO algorithm results, the average of the 30 runs should be considered and not only the best solution, which represents the best case. Therefore, the average fitness values for the previous six problems are calculated. The average fitness results are presented in Table 8.
Based on the results, NHHO has ranked first, outperforming all the other compared algorithms in an average of 30 runs. First, there is a significant difference between the average fitness of the proposed NHHO and the original HHO. Accordingly, NHHO has more flexibility to find the best solution and avoid falling into local optima. Second, comparing NHHO with PSO, ALO, BOA, and EO optimization algorithms, the significant improvement can be seen clearly in all problems especially in problem 1 and problem 5 Further, it is crucial to maintain consistency in the 30 runs. In other words, to achieve consistency the average of solutions has to be close, if not the same, to the best solution. The consistency is maintained in five problems in NHHO except for problem 6.
On the other hand, convergence speed is another critical issue in optimization algorithms. The proposed NHHO improves the convergence speed of HHO by avoiding premature convergence and having the right balance between the exploration and exploitation phases. A comparison between the NHHO algorithm and other optimization algorithms is shown in Figure 8.
It can be noticed from the results that NHHO can reach the best solution faster than all the other optimization algorithms, which indicates the superiority of NHHO. Further, the NHHO outperformed all the optimization algorithms in all problems except for problem 3, where NHHO has almost the same fitness as PSO. However, based on Figure 8, NHHO reached the solution much faster than PSO. The average of reaching the best solution in 30 runs for NHHO is 17 iterations, while the same fitness value needed an average of 72 iterations for PSO.
The convergence of the new hybrid algorithm is due to the combination of Newton's method, as a local search, and HHO algorithm. On one hand, Newton's method is designed for solving nonlinear systems. On the other hand, HHO guarantees a very good selection of initial guesses for Newton's method. Further, it is well-known that Newton's method is considered as a second-order convergence, which means that the number of correct digits roughly doubles in every step [41]. This advantage does not exist in any optimization algorithm, even in the HHO algorithm without utilizing Newton's method.
However, each optimization algorithm suffers from some limitations that make it difficult for the algorithm to locate the best solution. PSO has some limitations, such as poor population diversity and balancing between local optima. BOA gets stuck into local optima and it lacks exploitation of local values [48]. The limitations of ALO include the trapping of ants in antlions' trip which prevents the algorithm from identifying the global optimum solution. In contrast, the EO algorithm does not work effectively for large-scale problems [63].
Newton's method has been approved to be very effective when combines with HHO algorithm. Figure 9 presents a comparison between NHHO and other optimization algorithms in the six problems' fitness value in an average of 30 runs. The y-axis in Figure 9 must be set to a logarithmic scale as the gap between the values is too big. The average fitness of NHHO is 4.2267E − 10, while it is 0.07935 for the original HHO, 0.048372 for PSO, 0.33967 for ALO, 0.167762 for BOA, and 0.001294 for EO.
It is worth mentioning that all calculations have been performed using variables with 16 digits of precisions. The 16 digits variable is the default for MATLAB software and it is time-saving comparing to more significant digits. However, the 16 digits may limit the calculation and affect the results in some cases. The MATLAB function 'vpa' can be used to increase the number of digits. Increasing the number of digits may enhance the experimental results, but it is undoubtedly a time-consuming process. A comparison of applying NHHO to the six benchmark problems with other optimization algorithms in terms of fitness value and computational time based on variables of 16, 32, and 64 digits is shown in Table 9.
The results show a significant improvement in the fitness for all tested problems. For example, a fitness of 3.357E-16 using 16 digits to 1.504E-39 using 32 digits and to 1.082E-71 using 64 digits for problem 3. However, the consumed time increased dramatically from 1.4 seconds using variables of 16 digits to 8.8 seconds and 9.1 seconds, when using variables of 32 digits and 64 digits, respectively.
The standard deviation has been calculated for all problems to examine the stability of the proposed NHHO in finding the best solution in 30 runs. As can be noticed in Tables 2-6, the value of the standard deviation is zero in Table 2 and close to zero in Table 3-6. The results of NHHO outperformed the other algorithms which confirm the stability of NHHO in finding the best solutions in the benchmark. In benchmark problem 6, Table 7, NHHO came second with a standard deviation of ±0.04470 slightly higher than the original HHO that calculated ±0.02190.
However, every improvement comes with a cost. Despite the efficiency of the proposed hybrid, the overall complexity is quite higher if compared to pure algorithms. This is a common cause within hybrids where more than one algorithm complexity is added to another. While optimizing a problem solution, there should be a tradeoff between gain and loss factors. The goal of proposing such a hybrid was to enhance these three factors (best and average fitness and the convergence speed) at the expense of a slight increase in the required time and the complexity in finding the optimal solution. VOLUME 9, 2021

A. COMPARISON WITH NEWTON'S METHOD
To clarify the efficiency of NHHO on the accuracy of the obtained solutions, and its ability to avoid local optima compared to Newton's method, we applied both techniques to problems 1-4. A comparison between the fitness values obtained by NHHO and Newton's method using three arbitrary initial points is shown in Tables 10-13. For comparison purposes, we tested both schemes at the 5 th , 7 th , and 10 th iteration (Iter.). Besides, the best solution achieved by NHHO, in 10 runs and 1200 digits variables, was selected for each problem.
Tables [10][11][12][13] show that NHHO outperformed Newton's method in all selected problems. The main critical issue while applying Newton's method is the choice of the   initial point because the method is very sensitive to the selected initial point. Choosing an improper initial point may affect Newton's method by slowing its convergence (see Table 11 -12), or even make the method diverge (see Table 13), and sometimes the selected initial point may cause a singularity to the Jacobian matrix of F, therefore, the inverse of the Jacobian does not exist. Consequently, Newton's method cannot be applied (see Table 10).
The proposed NHHO overcomes these problems by selecting random multiple initial points, which are called search agents. NHHO tests the fitness of all search agents and the best search agent, which achieves the lowest fitness, is selected as an initial point. Accordingly, the selection of the best search agents plays an important role in overcoming the problem of choosing a proper initial point and eventually the solution's accuracy.    The results from Tables 10-13 also confirmed the dominance of NHHO over Newton's method. The reason for the significant improvement in the results is the ability of the hybrid NHHO to overcome the limitation of the initial VOLUME 9, 2021 guess in Newton's method by relying on 20 search agents at the initial stage instead of one as Newton's method does.
In these experiments, we fixed the variable to 1200 digits, and that is the reason why the fitness results for NHHO in iteration (Iter.) 10 are equal E − 1207 and E − 1208. Accordingly, the results may be improved as we increase the number of digits as mentioned previously.

B. COMPARISON WITH RELATED WORK
In this section, NHHO is compared with the related work of El-Shorbagy [18]. The authors in the related work proposed a hybrid technique (GOA-GA) that combines grasshopper optimization algorithm with genetic algorithm to solve systems of nonlinear equations. Hybrid-GOA-GA take the merits of both GOA and GA; where GOA's exploitability and GOA's exploration potential are combined.
Further, this work performed the experiments on the same problems (benchmarks) used in El-Shorbagy's work. Therefore, a comparison between the results of both algorithms for the six problems based on the best solution achieved in 30 runs was performed, and the achieved fitness values are presented in Table 14.
As noted from the results, NHHO outperformed the hybrid GOA-GA in all benchmark problems. The consistency of GOA-GA is reduced as the initial interval is expanded. For example, the fitness value increased from 5.3607E-09 to 2.7607E-04 when expanding the initial interval to the length of 10000 in problem 1. The increase in fitness value is noticed in all other problems as well, which affected the stability of the approach to be generalized to other benchmark problems. However, NHHO is more consistent as explained in the next section.

C. CONSISTENCY ON DIFFERENT INITIAL INTERVALS
It is critical to study the performance of the proposed NHHO over different initial intervals. Therefore, NHHO was run on every benchmark problem for 10 times at the following initial intervals  Table 15 shows the best fitness achieved for all problems. It is very clear, based on the results, that NHHO shows high consistency by achieving the same fitness values in all initial intervals.
The hybrid approach has a very strong search mechanism that applies Newton's method to search for the best solution even in a large search space such as in applying the initial interval of [−10000, 10000]. Although the initial intervals have widely expanded, NHHO managed to get the same best solutions mentioned in Tables 2-7.
However, the number of times NHHO reached the best solutions were not the same; while the best solutions were achieved almost in all of the 10 runs in the initial interval of [−10, 10], they were only achieved few times in the interval [−10000, 10000]. Accordingly, the proposed NHHO can avoid trapping into local optima even in large search spaces and that approves the superiority of NHHO. The complexity of NHHO is calculated based on adding the complexity of HHO to the complexity of Newton's method. It is worth mentioning here that the complexity of Newton's method is considered too high compared to the complexity of optimization algorithms. Therefore, it is expected that combining Newton's method with any optimization algorithm will increase its complexity by O(T ×W ) as Newton's method performs an additional local search. The computational times of HHO, Newton's method, and NHHO based on 100 iterations (T = 100) are shown in Table 16. However, to have a fair comparison, the experiments on NHHO performed on a stopping criterion that the achieved fitness value is less than 1 × E-15 (ε ≤ 1 × E-15) as it is well known that Newton's method can reach a significant fitness value with few iterations. Note that an N by N system of linear equations must be solved at each iteration. Thus, such derivative-based methods can be time-consuming because any computation of the Jacobian needs n^2 evaluations of scalar functions.
Combining HHO with Newton's method to develop NHHO did not add significant computational time compared to Newton's method as can be noticed from Table 16. However, NHHO can overcome the limitations of Newton's method, as discussed previously, such as selecting the initial points and divergence problems. Therefore, it is more effective in solving nonlinear equation systems.

VI. CONCLUSION
In this work, a hybrid between a powerful optimization algorithm (HHO) and an iterative approximation method (Newton's method) was proposed to solve systems of nonlinear equations. The role of the HHO was to guarantee a good selection of initial guesses for Newton's method to widen its applicability to different problems where Newton's method was employed as a local search mechanism with the advantage of being a second-order convergence method.
In the experiments, we utilized six well-known problems that consist of complex systems of nonlinear equations to evaluate the performance of NHHO. Besides, NHHO was compared to well-known optimization algorithms such as the original HHO, PSO, EO, ALO, and BOA. The experiments were evaluated based on each optimization algorithm's fitness value, which was calculated based on the norm of the equations. The results have confirmed the superiority of NHHO over other algorithms in solving systems of nonlinear equations in terms of the best solution, fitness value, and convergence speed. Besides, the results of running the algorithms 30 times confirmed the consistency of NHHO.
Moreover, the nonlinear system's solution can be significantly improved by increasing the number of digits for the variables used by NHHO but it is time-consuming. The consistency of NHHO has been also approved on different initial intervals. Finally, NHHO overcomes the limitations of Newton's method including initial point selection and divergence problems.
The scope of this manuscript is to evaluate the ability of the proposed NHHO algorithm in solving systems of nonlinear equations and hypothetically, all the solved problems were considered as global optimization problems since no studies have provided a clear categorization of them. The performance of the proposed algorithm in solving standard optimization benchmarks may be further investigated in the future.
The motivation of this work is to solve systems of nonlinear equations. However, the effect of NHHO on other types of problems, such as unimodal and multimodal problems will be considered for future work.
We believe that the local search of Newton's method can be applied to other optimization algorithms to improve their search mechanisms for future work. Besides, NHHO can be applied to real-life problems such as data science and engineering problems. KHAIRUDDIN OMAR received the bachelor's and master's degrees in computer science from Universiti Kebangsaan Malaysia, in 1986 and 1989, respectively, and the Ph.D. degree in philosophy from Universiti Putra Malaysia, in 2000. He is currently a Professor with the Faculty of Information Science and Technology, Universiti Kebangsaan Malaysia. His research interests include artificial intelligence (pattern recognition in decision making with uncertainty-Bayesian reasoning; neural networks; fuzzy logic; and fuzzy neural networks, image processing-2D and 3D, edge detection, thinning, segmentation, feature extraction, image improvement, texture, resolution, image transforms, e.g., Trace, Fourier, and Wavelet) with applications to Jawi/Arabic manuscripts and biometric authentication.
KHAIRUL AKRAM ZAINOL ARIFFIN (Member, IEEE) received the bachelor's and master's degrees (Hons.) in system engineering with computer engineering from the University of Warwick, U.K., in 2008 and 2009, respectively. He later joined Universiti Teknologi PETRONAS (UTP), in 2010, to pursue his journey towards academic research and teaching courses to earn the Ph.D. degree in information systems. During his time at UTP, many journal articles and conference papers have been produced and published internationally. Then, he was appointed as a Researcher with the Digital Forensic Department, CyberSecurity Malaysia, and had been entrusted with the research on embedded systems and live forensic. He is currently a member of the Faculty Information Science and Technology, Universiti Kebangsaan Malaysia, to pursue his passion for research toward cybersecurity, digital forensics, algorithms, and embedded systems. He is GCFA certified and a member of IET.