# Survey of Methods for Solving Systems of Nonlinear Equations, Part II:
Optimization Based Approaches

Ilias S. Kotsireas
<EMAIL_ADDRESS>
Wilfrid Laurier University
Canada  Panos M. Pardalos
<EMAIL_ADDRESS>
Toxeus Systems LLC
Orlando, Florida, USA  Alexander Semenov
<EMAIL_ADDRESS>
University of Florida
Gainesville, Florida, USA  William T. Trevena
<EMAIL_ADDRESS>
University of Florida
Gainesville, Florida, USA  Michael N. Vrahatis
<EMAIL_ADDRESS>
University of Patras
Patras, Greece

###### Abstract

This paper presents a comprehensive survey of methods which can be utilized to
search for solutions to systems of nonlinear equations (SNEs). Our objectives
with this survey are to synthesize pertinent literature in this field by
presenting a thorough description and analysis of the known methods capable of
finding one or many solutions to SNEs, and to assist interested readers
seeking to identify solution techniques which are well suited for solving the
various classes of SNEs which one may encounter in real world applications.

To accomplish these objectives, we present a multi-part survey. In part one,
we focused on root-finding approaches which can be used to search for
solutions to a SNE without transforming it into an optimization problem. In
part two, we introduce the various transformations which have been utilized to
transform a SNE into an optimization problem, and we discuss optimization
algorithms which can then be used to search for solutions. We emphasize the
important characteristics of each method, and we discuss promising directions
for future research. In part three, we will present a robust quantitative
comparative analysis of methods capable of searching for solutions to SNEs.

Keywords: systems of nonlinear equations, global optimization, localization of
zeros, computation of roots, metaheuristics, topological degree, total number
of solutions and extrema

##  1 Introduction

This is the second part of a survey on methods for finding one or many real
solutions to a system of nonlinear equations (SNE):

${F}_{m}(x)={\mathit{\Theta}}_{m}\equiv(0,0,\ldots,0)^{\top}\quad\Longleftrightarrow\quad\left\\{\begin{aligned}
&f_{1}(x_{1},x_{2},\ldots,x_{n})=0,\\\\[2.84544pt]
&f_{2}(x_{1},x_{2},\ldots,x_{n})=0,\\\ &\kern 41.25641pt\vdots\\\
&f_{m}(x_{1},x_{2},\ldots,x_{n})=0,\end{aligned}\right.$ (1)

where
${F}_{m}=(f_{1},f_{2},\ldots,f_{m}):{\mathcal{D}}_{n}\subset{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$,
where $f_{1},f_{2},\ldots,f_{m}$ are real-valued continuous or continuously
differentiable functions on the domain ${\mathcal{D}}_{n}$, and where at least
one of $f_{1},f_{2},\ldots,f_{m}$ is nonlinear. For example, consider the
system of transcendental equations

${F}_{2}(x)={\mathit{\Theta}}_{2}\equiv(0,0)^{\top}\quad\Longleftrightarrow\quad\left\\{\begin{aligned}
f_{1}(x_{1},x_{2})&=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0,\\\\[2.84544pt]
f_{2}(x_{1},x_{2})&=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0,\end{aligned}\right.$
(2)

which is comprised of two transcendental equations of two unknowns (See Figure
1).

Figure 1: An example of a SNE with two transcendental equations of two
unknowns as introduced by Eq. (2):
(Blue):
$f_{1}(x_{1},x_{2})=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0$;
(Red):
$f_{2}(x_{1},x_{2})=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0$.
Solutions to this SNE are defined as the points where the blue and red
contours intersect. Finding all of the points within a certain region which
satisfy both equations is a challenging task.

Finding one or more solutions to a SNE is a challenging and ubiquitous task
faced in many fields including chemistry [1, 2, 3], chemical engineering [4],
automotive steering [5], power flow [6, 7], large-scale integrated circuit
designs [8], climate modeling [9], materials engineering [10], robotics [11,
12, 13, 14], nuclear engineering [15], image restoration [16], protein
interaction networks [8], neurophysiology [17], economics [18], finance [19],
applied mathematics [20], physics [21], finding string vacua [22], machine
learning [23, 24], and geodesy [25, 26] among others. The problem of solving
even a system of polynomial equations has been proven to be NP-hard [27].
Furthermore, it has also been proven [28] that no general algorithm exists for
determining whether an integer solution exists for a polynomial equation with
a finite number of unknowns and only integer coefficients. The latter has been
known as Hilbert’s 10th problem.

###  1.1 Notation / Scientific Style

Throughout this paper, we utilize
$x=(x_{1},x_{2},\ldots,x_{n})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$ to
denote a real vector within the bounded domain $\mathcal{D}$. Furthermore, we
utilize
$x^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$
to denote a real solution to a SNE such that all equations in the SNE are
satisfied ($F_{m}(x^{*})=0$). In an iterative method, we utilize
$x^{k}=(x_{1}^{k},x_{2}^{k},\ldots,x_{i}^{k},\ldots
x_{n}^{k})^{\top}\in\mathcal{D}\subset\mathbb{R}^{n}$ for $k=0,1,\ldots$ to
denote the vector found during the $k-$th iteration of the iterative method.
Here, $x^{k}_{i}$ denotes the $i-$th coordinate of the vector $x^{k}$. When we
discuss optimization methods, we denote an objective function as
$\varphi:\mathbb{R}^{n}\rightarrow\mathbb{R}$, and we denote the corresponding
gradient as $\nabla\varphi(x)$.

###  1.2 Terminology

Although we refer to Eq. (2) as a system of nonlinear equations (SNE), such
systems have been referred to in a variety of different ways in literature.
For example, articles [29, 30, 31, 32, 33, 34] utilize the abbreviation “SNLE”
to refer to a system of nonlinear equations, and article [35] uses the
abbreviation “SoNE”. Other papers refer to Eq. (2) as a nonlinear system of
equations, and use the abbreviations “NSE” [36, 37] and “NLS” [38]. Eq. (2)
has also been referred to as a nonlinear equation system (NES) [39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]. The survey in article [54]
uses the terminology Nonlinear Equations (NEs) to refer to a system of one or
more nonlinear equations.

When $m>n$, a SNE can be referred to as an overdetermined SNE, and when $n>m$,
a SNE can be referred to as an underdetermined SNE. When $m=n$, a SNE can be
referred to as a square SNE [55]. Furthermore, a SNE is considered to be
consistent if a solution exists which satisfies all equations [56].

For a square SNE, a solution
$x^{*}=(x_{1}^{*},x_{2}^{*},\ldots,x_{n}^{*})^{\top}$ of the SNE
${F}_{n}(x)={\mathit{\Theta}}_{n}$ or equivalently a zero $x^{*}$ of the
function ${F}_{n}(x)$ or a root $x^{*}$ of the function ${F}_{n}(x)$ is called
simple if for the determinant of the corresponding Jacobian matrix:

$J_{F_{n}}(x)\equiv
F^{\prime}_{n}(x)_{ij}\equiv\left\\{\frac{\displaystyle\partial
f_{i}(x)}{\displaystyle\partial
x_{j}}\right\\}_{ij}\equiv\left[\begin{array}[]{cccc}\vspace*{0.1cm}{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{1}}}&{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{2}}}&\cdots&{\frac{\displaystyle\partial
f_{1}(x)}{\displaystyle\partial x_{n}}}\\\\[5.69046pt]
{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{1}}}&{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{2}}}&\cdots&{\frac{\displaystyle\partial f_{2}(x)}{\displaystyle\partial
x_{n}}}\\\ \vdots&\vdots&\ddots&\vdots\\\\[2.84544pt]
{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{1}}}&{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{2}}}&\cdots&{\frac{\displaystyle\partial f_{n}(x)}{\displaystyle\partial
x_{n}}}\end{array}\right],\vspace*{0.1cm}$ (3)

at $x^{*}$ it holds that ${\rm det}\,J_{F_{n}}(x^{*})\neq 0$, otherwise it is
called multiple. The problem of conservation and decomposition of a multiple
root into simple roots in the case of systems of homogeneous algebraic
equations has been tackled in [57]. This approach can be applied to high
dimensional CAD where it is sometimes required to compute the intersection of
several hypersurfaces that are a perturbation of a set of original unperturbed
hypersurfaces.

When $F_{n}$ satisfies the monotonicity condition:

$\bigl{(}F_{n}(x)-F_{n}(y)\bigr{)}^{\top}(x-y)\geqslant 0,\kern
8.5359pt\forall\,x,y\in\mathbb{R}^{n},$ (4)

the corresponding SNE can be referred to as a system of monotone nonlinear
equations [58]. Furthermore, $F_{n}$ is considered to be Lipschitz continuous
if there exists $L>0$ such that

$\|F_{n}(x)-F_{n}(y)\|_{2}\leqslant L\|x-y\|_{2},\kern
5.69046pt\forall\,x,y\in\mathbb{R}^{n}.$ (5)

###  1.3 Comparison to other surveys

Other surveys discussing solution techniques for SNEs include [59] and [54].
We have decided to conduct this comprehensive literature review because many
new solution techniques for SNEs have been introduced since the publication of
[59] in 1994, and because the recent survey presented in [54] focuses mainly
on methods which first convert a SNE into an optimization problem, and then
search for multiple solutions to the optimization problem using Intelligent
Optimization Algorithms (IOAs). The IOAs discussed in article [54] are
primarily metaheuristics for global optimization. Although the survey in
article [54] provides a very nice discussion of IOAs for solving SNEs
reformulated as optimization problems, many of the IOAs they discuss are only
introduced at a very high level, only eight IOAs were tested in their
computational study, and the IOAs were evaluated on SNEs comprised of 20
equations or less. Also, article [54] only briefly mentions methods which can
be used to search for solutions to SNEs without transforming them into
optimization problems.

We would like to present a broader survey which covers in detail the large set
of methods which can be used to solve a SNE without transforming it into an
optimization problem (i.e. homotopy and symbolic computation methods). These
methods were our main focus in part one of this survey. In part two, we expand
upon article [54] by introducing additional reformulation techniques and
optimization algorithms which have been used to solve SNEs, and by discussing
in much more detail many optimization algorithms which were only briefly
introduced in article [54]. This will allow us to appropriately set the stage
for the comprehensive empirical study we will present in part three of this
survey. Furthermore, we believe it is imperative to introduce the reader to a
technique for determining the number of solutions to a SNE that exist within a
bounded domain. Such techniques are of critical practical importance for those
interested in finding all solutions to a SNE that exist within a domain of
interest.

###  1.4 Organization of this survey

Before we discuss optimization methods which can be used to search for
solutions to SNEs, we first introduce a method which can be used to determine
the total number of solutions to a SNE that exist within a given bounded
domain. If one can determine the number of solutions to a SNE which exist
within a bounded domain of interest, one can confidently define an appropriate
stopping criteria for algorithms which can be used to search for solutions.
Next, we will discuss the various ways that a SNE can be transformed into an
optimization problem, and we will introduce techniques that can be utilized to
search for solutions to the global optimization problem that arises when the
most common reformulation is performed. We will conclude our paper by
highlighting promising areas for future research.

##  2 Determining the number of solutions to a SNE in a bounded domain

The knowledge of all the solutions of a system of nonlinear equations and/or
all the extrema of a function is of major importance in various fields. The
total number of the solutions of a system of nonlinear equations can be
obtained by computing the topological degree. Suppose that the function
$F_{n}=(f_{1},f_{2},\ldots,f_{n})\colon{\mathcal{D}}_{n}\subset$
${\mathbb{R}}^{n}$ $\to$ ${\mathbb{R}}^{n}$ is defined and is two times
continuously differentiable in a bounded domain ${\mathcal{D}}_{n}$ of
${\mathbb{R}}^{n}$ with boundary $b({\mathcal{D}}_{n})$. Suppose further that
the solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ are not located on
$b({\mathcal{D}}_{n})$, and that they are simple (that the determinant of the
Jacobian of $F_{n}$ at these solutions is non-zero). Then the topological
degree of $F_{n}$ at ${\mathit{\Theta}}_{n}$ relative to ${\mathcal{D}}_{n}$
is denoted by ${\rm deg}[F_{n},{\mathcal{D}}_{n},$ ${{\mathit{\Theta}}_{n}}]$
and can be defined by the following relation:

${\rm
deg}[F_{n},{\mathcal{D}}_{n},{{\mathit{\Theta}}_{n}}]=\displaystyle\sum_{{\scriptstyle
x\in F_{n}^{-1}({{\mathit{\Theta}}_{n}})}}{\rm sgn}\,{\rm det}\,J_{F_{n}}(x),$

where ${\rm det}\,J_{F_{n}}(x)$ denotes the determinant of the Jacobian matrix
and ${\rm sgn}$ defines the three-valued sign function. The above definition
can be generalized when $F_{n}$ is only continuous [60].

It is evident that, since ${\rm deg}[F_{n},{\mathcal{D}}_{n},$
${{\mathit{\Theta}}_{n}}]$ is equal to the number of simple solutions of
$F_{n}(x)={\mathit{\Theta}}_{n}$ which give positive determinant of the
Jacobian matrix, minus the number of simple solutions which give negative
determinant of the Jacobian matrix, then the total number $N^{s}$ of simple
solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ can be obtained by the value of
${\rm deg}[F_{n},{\mathcal{D}}_{n},$ ${\mathit{\Theta}}_{n}]$ if all these
solutions have the same sign of the determinant of the Jacobian matrix. Thus,
Picard considered the following extensions of the function $F_{n}$ and the
domain ${\mathcal{D}}_{n}$ [61, 62]:

$F_{n+1}=(f_{1},f_{2},\ldots,f_{n},f_{n+1})\colon{\mathcal{D}}_{n+1}\subset{\mathbb{R}}^{n+1}\to{\mathbb{R}}^{n+1},$
(6)

where $f_{n+1}=y\,\,{\rm det}\,J_{F_{n}}$,
${\mathbb{R}}^{n+1}:x_{1},x_{2},\ldots,x_{n},y$, and ${\mathcal{D}}_{n+1}$ is
the direct product of the domain ${\mathcal{D}}_{n}$ with an arbitrary
interval of the real $y$-axis containing the point $y=0$. Then the solutions
of the following system of equations:

$\begin{array}[]{l}f_{i}(x_{1},x_{2},\ldots,x_{n})=0,\quad
i=1,2,\ldots,n,\\\\[2.84544pt] y\,\,{\rm
det}\,J_{F_{n}}(x_{1},x_{2},\ldots,x_{n})=0,\end{array}$

are the same simple solutions of $F_{n}(x)={\mathit{\Theta}}_{n}$ provided
that $y=0$. Obviously, the determinant of the Jacobian matrix obtained for the
function (6) is equal to $({\rm det}\,{J_{F_{n}}}(x))^{2}$ which is always
positive. Thus, the total number $N^{s}$ of the solutions of the system
$F_{n}(x)={\mathit{\Theta}}_{n}$ can be obtained by the following value of the
topological degree:

$N^{s}={\rm deg}[F_{n+1},{\mathcal{D}}_{n+1},{\mathit{\Theta}}_{n+1}].$

For example, in the one dimensional case, using the above Picard’s extensions
it is proved that the total number of simple solutions $N^{s}$ of the equation
$f(x)=0$, where $f\colon(a,b)\subset$ ${\mathbb{R}}$ $\to$ ${\mathbb{R}}$ is
twice continuously differentiable in a predetermined interval $(a,b)$, is
given by the following relation [61, 62]:

${N^{s}}=-\frac{1}{\pi}\left[\varepsilon\int_{a}^{b}\frac{f(x)\,f^{\prime\prime}(x)-{f^{\prime}}^{2}(x)}{f^{2}(x)+{\varepsilon}^{2}{f^{\prime}}^{2}(x)}\,dx+\arctan\left(\frac{\varepsilon
f^{\prime}(b)}{f(b)}\right)-\arctan\left(\frac{\varepsilon
f^{\prime}(a)}{f(a)}\right)\right],$ (7)

where $\varepsilon$ is a small positive constant. Note that ${N^{s}}$ was
shown to be independent of the value of $\varepsilon$. Also, the above
approach can be applied for computing the number of multiple solutions.
Obviously, the total number $N^{e}$ of the extrema of $f\in C^{3}$ i.e.
$x\in(a,b)$ such that $f^{\prime}(x)=0$ can be obtained using the above
formula (7). For details of the topological degree we refer the interested
reader to the books [63, 64, 65, 66, 60, 67]. Details of the computation of
the value of the topological degree and its usefulness as well as some
applications and issues related to the number of zeros can be found for
example in [68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82].

Article [83] discusses the robust satisfiability of SNEs, and discusses
utilizing the topological degree to determine the existence of robust
solutions to a SNE.

##  3 Transforming a SNE into an optimization problem

Although many of the root finding methods described in part one of this survey
are guaranteed to converge to a solution on certain classes of SNEs (such as
SNEs that satisfy the conditions of monotonicity (cf. Eq. (4)) and Lipschitz
continuity (cf. Eq. (5))), many of the root finding methods described in part
one of this survey are incapable of reliably finding solutions to SNEs which
do not satisfy such conditions. As a result, many approaches for solving
general classes of SNEs first reformulate a SNE as an optimization problem,
and then utilize optimization algorithms to attempt to find solutions.
However, it is important to note that the resulting optimization problem may
be multimodal and potentially nonsmooth. Therefore, global optimization
techniques are often utilized to attempt to find solutions. For the interested
reader, the books [84, 85] discuss feasibility and infeasibility in
optimization.

From the best of our knowledge, this section discusses the various techniques
which have been presented in literature for reformulating a SNE as an
optimization problem. In literature, techniques have been presented for
transforming a SNE into a single objective global optimization problem
(Section 3.1), into a multiobjective optimization problem (Section 3.3), into
a constrained optimization problem (Section 3.2), and into a nonlinear
programming problem (Section 3.4). Article [54] provides a nice analysis and
comparison of many of the reformulation techniques described in this section,
especially in regards to the techniques they use to support efforts aimed at
finding multiple solutions to a SNE.

###  3.1 Transforming a SNE into a single-objective global optimization
problem

A nonlinear system of equations (Eq. (1)) may be transformed into a single-
objective global optimization problem by introducing the objective function

$\varphi^{0}(x)=\alpha\sum_{i=1}^{m}|f_{i}(x)|^{p}$ (8)

where $\alpha>0$ and $p>0$ (See Figure 2). Alternate formulations of
$\varphi^{0}(x)$ have also been proposed in literature, for example, article
[86] proposes utilizing $\varphi^{0}(x)=\sqrt{\sum_{i=1}^{m}f_{i}^{2}(x)}$. In
order to find solutions to a SNE (Eq. (1)), the following global optimization
problem should be solved:

$\min_{x\in{\mathcal{D}}}\varphi^{0}(x).$ (9)

That is, to find the value of the minimizer $x$:

$\arg\min_{x\in{\mathcal{D}}}\varphi^{0}(x),$

for which $\varphi^{0}(x)$ attains its minimum value.

Although almost all papers which reformulate SNEs as optimization problems
reformulate them as minimization problems, SNEs can also be transformed into
maximization problems. Article [87] proposes transforming SNEs into
maximization problems by maximizing over an objective function of the form
$\varphi(x)={1}/({1+\sum_{i=1}^{m}|f_{i}(x)|})$. In this reformulation,
solutions to a SNE are represented by points with an objective function value
of 1. Article [87] is the only work that we are aware of that reformulates
SNEs as maximization problems. For the purposes of this paper, when we discuss
reformulating a SNE as an optimization problem, we are referring to the
standard case of reformulating a SNE as a minimization problem such as Eq.
(9).

(a) $\,\,\alpha=1,\,p=1$ (b) $\,\,\alpha=1,\,p=2$

Figure 2: The objective function surface
$\varphi^{0}(x)=\alpha\sum_{i=1}^{m}|f_{i}(x)|^{p}$ when $\alpha=1$, $p=1$ and
$\alpha=1$, $p=2$ for the example SNE introduced in Eq. (2) with two
transcendental equations and two unknowns:
$f_{1}(x_{1},x_{2})=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0$;
$f_{2}(x_{1},x_{2})=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0$.

Since systems of nonlinear equations may have more than one solution, the
reformulated optimization problem Eq. (9) may have more than one global
minimizer, and each global minimizer with an objective function value of 0
(global minimum) will correspond to a solution to the SNE. Although the
problem shown by Eq. (9) can be repeatedly solved to find many solutions to a
SNE, depending on the solution technique used, the same solution(s) to a SNE
may be found repeatedly during the search for new unfound solutions. To
address this, a plethora of alternate reformulations have been utilized in
literature which are specifically designed to help algorithms find further
solutions to a SNE. These reformulations typically introduce a penalty (also
referred to as a “repulsion” [46] or “polarization” [5] technique) into the
objective function to repel the search for new solutions away from solutions
which have already been found. These reformulations are often of the form:

$\varphi(x)=\varphi^{0}(x)+\varphi^{k}(x),$ (10)

where $\varphi^{k}(x)$ creates a penalty region around the $k$ already found
solutions (See Figure 3). Article [37] proposed to represent $\varphi^{k}(x)$
as:

$\varphi^{k}(x)=\sum_{j=0}^{k}{a\,\xi\left(\frac{\|x-x^{j}\|_{2}}{\rho}\right)},$
(11)

where $x^{j}$ are already found solutions, $k$ is the number of solutions
found so far, $a>0$ and $\rho>0$ are user defined parameters, and
$\xi(\delta)$ is a continuous real function that is positive for $|\delta|<1$
and zero elsewhere. Article [37] proposes to utilize the following variants of
$\xi(\delta)$:

  * (a)

$\xi(\delta)=1-|\delta|$ for $|\delta|<1$, otherwise 0,

  * (b)

$\xi(\delta)=(1-\delta^{2})^{2}$ for $|\delta|<1$, otherwise 0,

  * (c)

$\xi(\delta)=(1+\cos(\pi\delta))/2$ for $|\delta|<1$, otherwise 0,

  * (d)

$\xi(\delta)=\exp\bigl{(}-{\delta^{2}}/(1-\delta^{2})\bigr{)}$ for
$|\delta|<1$, otherwise 0.

(a) $\,\,\alpha=1,\,p=1,\,a=1,\,\rho=1$ (b)
$\,\,\alpha=1,\,p=2,\,a=4,\,\rho=1$

Figure 3: The penalized objective function
$\varphi(x)=\varphi^{0}(x)+\varphi^{k}(x)=\alpha\sum_{i=1}^{m}|f_{i}(x)|^{p}+\sum_{j=1}^{k}{a\,\xi\left({\|x-x^{j}\|_{2}}/{\rho}\right)}$
(where $\xi(\delta)=(2-\delta^{2})^{2}$ for $|\delta|<\sqrt{2}$ and $0$
otherwise) with a penalty $\varphi^{k}(x)$ for a single found solution ($k=1$)
$x^{1}$ at $(0,0)$ when $\alpha=1,p=1,a=1,\rho=1$ and
$\alpha=1,p=2,a=4,\rho=1$ for the example SNE introduced in Eq. (2):
$f_{1}(x_{1},x_{2})=x_{1}-x_{1}\sin(x_{1}+5x_{2})-x_{2}\cos(5x_{1}-x_{2})=0$;
$f_{2}(x_{1},x_{2})=x_{2}-x_{2}\sin(5x_{1}-3x_{2})+x_{1}\cos(3x_{1}+5x_{2})=0$.
Notice the sharply elevated region surrounding the point $(0,0)$ on both
objective function surfaces.

Notice that if there are no solutions found yet, then Eq. (10) reduces to the
original optimization problem Eq. (9). A similar approach is suggested in the
papers [88, 89] where

$\varphi^{k}(x)=\beta\sum_{j=0}^{k}{\exp\bigl{(}-\|x-x^{j}\|_{2})\,\chi(\|x-x^{j}\|_{2}\bigr{)}}.$
(12)

Here, $\beta$ is a large constant, $\varrho$ is a small constant,
$\chi(\delta)=1$ when $\delta<\varrho$ and otherwise $\chi(\delta)=0$. For the
rest of the reformulations introduced in Section 3.1, we assume that at least
one solution has been found so far ($k\geqslant 1)$.

Article [5] proposes minimizing an objective function of the form

$\varphi(x)=\dfrac{\varphi^{0}(x)}{\Pi_{j=1}^{k}\arctan\|x-x^{j}\|_{2}},$ (13)

while article [90] proposes minimizing an objective function of the form

$\varphi(x)=\bigl{(}\varphi^{0}(x)+\varepsilon\bigr{)}\,\Pi_{j=1}^{k}\bigl{|}\coth(\alpha\|x-x^{j}\|_{2})\bigr{|},$
(14)

where $\varepsilon>0$ is a small user defined constant and $\alpha\geqslant 1$
is a user defined parameter which is utilized to adjust the radius of the
penalty region.

Article [91] proposes minimizing an objective function of the form

$\varphi(x)=\varphi^{0}(x)\,\Pi_{j=1}^{k}\xi(x,x^{j},\alpha,\rho),$ (15)

where

$\displaystyle\xi(x,x^{j},\alpha,\rho)=|\,\text{erf}(\alpha\|x-x^{j}\|_{2})|^{-1}\,\,\text{
if }\,\,\|x-x^{j}\|_{2}\leqslant\rho,\text{ and }\,0\text{ otherwise}.$

Here, $\rho>0$ and $\alpha>0$ are user defined parameters which adjust the
radius of the penalty region and the magnitude of the penalty respectively,
and $\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}\exp(-t^{2})dt$.

Article [44] proposes a similar reformulation to Eq. (15) with a dynamic
repulsion radius that is updated within a evolutionary algorithm, and article
[92] proposes minimizing an objective function of the form:

$\varphi(x)=\varphi^{0}(x)+\sum_{j=1}^{k}\xi(x,x^{j},\alpha,\rho),$ (16)

where

$\displaystyle\xi(x,x^{j},\alpha,\rho)=\alpha\bigl{(}1-\text{erf}(\|x-x^{j}\|_{2})\bigr{)}\,\text{
if }\,\|x-x^{j}\|_{2}\leqslant\rho,\text{ and }\,0\text{ otherwise},$

where $\rho>0$ and $\alpha>0$ are user defined parameters which adjust the
radius of the penalty region and the magnitude of the penalty respectively,
and where $\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}\exp(-t^{2})dt$.

Other single objective reformulations have been utilized for particular
applications. For example, articles [93, 94] propose minimizing over an
objective function of the form:

$\varphi(x)=\varphi^{0}(x)\,\Pi_{j=1}^{k}\left(\frac{\Pi_{i=1}^{n}((b_{i})^{k\alpha})}{\Pi_{i=1}^{n}|x_{i}-x^{k}_{i}|^{\alpha}}\right),$
(17)

where $a\leqslant x\leqslant b$, $x^{j}$ for $j=1,2,\ldots,k$ are already
found solutions ($k$ is the number of solutions found so far), $x_{i}$ is the
$i-$th coordinate of the solution $x$, and $\alpha>0$ is a user defined
parameter. Article [94] also proposes a slightly different objective function,
they propose minimizing over an objective function of the form

$\varphi(x)=\varphi^{0}(x)\,\Pi_{j=1}^{k}\exp\left(\frac{\Pi_{i=1}^{n}((b_{i})^{k\alpha})}{\Pi_{i=1}^{n}|x_{i}-x^{k}_{i}|^{\alpha}}\right).$
(18)

It is important to note that for the reformulations in Eq. (17) and Eq. (18),
if any of the coordinates of a new solution are the same as that of a previous
found solution, the objective function value will be infinity. By using this
reformulation, one can ensure that for each new solution that is found, every
dimension will be different from that of all previously found solutions.

###  3.2 Transforming a SNE into a single-objective constrained optimization
problem

In addition to the non-constrained reformulations we have described so far,
constrained reformulations include:

$\min\varphi(x)\,\,\text{ for }\,\,\varphi(x)=\sum_{i=1}^{m}f_{i}(x),$ (19)

subject to

$f_{i}(x)\geqslant 0,\kern 8.5359pt\forall\,i=1,2,\ldots,m,$

as presented in [34]. Alternatively, article [95] first partitioned the SNE
into two sets of equations $S_{1}=\\{1,2,\ldots,m\\}$ and
$S_{2}=\\{1,2,\ldots,m\\}/S_{1}$. Based on these two sets of equations,
article [95] then utilized the reformulation:

$\min\varphi(x)\,\,\text{ for }\,\,\varphi(x)=\sum_{i\in S_{1}}f_{i}^{2}(x),$
(20)

subject to

$f_{j}(x)=0,\kern 8.5359pt\forall\,j\in S_{2}.$

Constrained optimization based transformation techniques such as those
presented above require a constraint handling method to solve.

###  3.3 Transforming a SNE into a multiobjective optimization problem

Instead of transforming a SNE into a single-objective optimization problem, a
number of research papers (see for example [18, 40]) suggest transforming a
SNE (Eq. (1)) into a multiobjective optimization problem of the form

$\left\\{\begin{aligned} &\text{min
}|f_{1}(x_{1},x_{2},\ldots,x_{n})|,\\\\[5.69046pt] &\text{min
}|f_{2}(x_{1},x_{2},\ldots,x_{n})|,\\\ &\kern 48.36958pt\vdots\\\ &\text{min
}|f_{m}(x_{1},x_{2},\ldots,x_{n})|,\end{aligned}\right.$ (21)

which can be solved via multiobjective optimization methods. Alternatively,
instead of having an objective function for each equation, article [96]
randomly splits the $m$ equations into two subsets $S_{1}$ and $S_{2}$ such
that $S_{1}\cap S_{2}=\emptyset$ and $S_{1}\cup S_{2}$ is equivalent to the
SNE (1). In terms of these two sets, article [96] then transforms a SNE into a
biobjective optimization problem of the form:

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\sum_{i=1}^{|S_{1}|}|f_{S_{1}(i)}(x)|,\text{ s.t.
}f_{S_{1}(i)}(x)=0,\\\\[5.69046pt] &\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=\sum_{i=1}^{|S_{2}|}|f_{S_{2}(i)}(x)|,\text{ s.t.
}f_{S_{2}(i)}(x)=0,\\\ \end{aligned}\right.$ (22)

or

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\sum_{i=1}^{|S_{1}|}|f_{S_{1}(i)}(x)|,\text{ s.t.
}f_{S_{1}(i)}(x)\geqslant 0,\\\\[5.69046pt] &\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=\sum_{i=1}^{|S_{2}|}|f_{S_{2}(i)}(x)|,\text{ s.t.
}f_{S_{2}(i)}(x)\geqslant 0.\end{aligned}\right.$ (23)

Article [40] reformulates a SNE into a biobjective optimization problem of the
form:

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=x_{1}+\sum_{i=1}^{m}|f_{i}(x)|,\\\\[5.69046pt]
&\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=1-x_{1}+m\,\max|f_{i}(x)|,\\\ \end{aligned}\right.$ (24)

where $i=1,2,\ldots,m$, and $x_{1}$ represents the first decision variable of
$x$. In this reformulation, $x_{1}$ is utilized to ensure that the two
objective functions conflict. Note that when using this reformulation, it may
be difficult to detect multiple solutions that have the same value of $x_{1}$
[54]. Expanding upon this approach, article [39] presented a biobjective
reformulation of the form:

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\frac{\sum_{i=1}^{n}w_{i}\,x_{i}}{\sum_{i=1}^{n}w_{i}}+\sum_{j=1}^{m}|f_{j}(x)|,\\\\[5.69046pt]
&\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=1-\frac{\sum_{i=1}^{n}w_{i}\,x_{i}}{\sum_{i=1}^{n}w_{i}}+\sum_{j=1}^{m}|f_{j}(x)|,\\\
\end{aligned}\right.$ (25)

where $i=1,2,\ldots,n$, $w_{i}$ is a randomly generated weight for each
variable such that $\sum_{i=1}^{n}w_{i}=1$, $n$ is the number of decision
variables in $x$, and $m$ is the number of equations. By utilizing this
reformulation, the risk of not being able to identify two solutions with the
same value of $x_{1}$ inherent in reformulation (24) is decreased.

Also expanding upon reformulation (24), article [41] transforms a SNE into a
multiobjective optimization problem with $n+1$ objectives where $n$ is the
number of decision variables in the SNE. The reformulated multiobjective
optimization problem is of the form:

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\frac{x_{1}}{n}+\frac{x_{2}}{n-1}+\cdots+\frac{x_{n-1}}{2}+\frac{x_{n}}{1}+C\,\sum_{i=1}^{m}|f_{i}(x)|\,\ln(n+2),\\\
&\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=\frac{x_{1}}{n}+\frac{x_{2}}{n-1}+\cdots+\frac{x_{n-1}}{2}+(1-x_{n})+C\,\sum_{i=1}^{m}|f_{i}(x)|\,\ln(n+1),\\\
&\min\varphi^{3}(x)\,\,\text{ for
}\,\,\varphi^{3}(x)=\frac{x_{1}}{n}+\frac{x_{2}}{n-1}+\cdots+\frac{x_{n-2}}{3}+(1-x_{n-1})+C\,\sum_{i=1}^{m}|f_{i}(x)|\,\ln(n),\\\
&\kern 159.90434pt\vdots\\\ &\min\varphi^{n}(x)\,\,\text{ for
}\,\,\varphi^{n}(x)=\frac{x_{1}}{n}+(1-x_{2})+C\,\sum_{i=1}^{m}|f_{i}(x)|\,\ln(3),\\\
&\min\varphi^{n+1}(x)\,\,\text{ for
}\,\,\varphi^{n+1}(x)=(1-x_{1})+C\,\sum_{i=1}^{m}|f_{i}(x)|\,\ln(2).\end{aligned}\right.$
(26)

In article [47], a SNE is transformed into a biobjective optimization problem
of the form:

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\sum_{i=1}^{m}[f_{i}(x)]^{2},\\\\[5.69046pt]
&\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=-\frac{1}{S_{n}}{\sum_{j=1}^{S_{n}}K(x,x^{j})},\\\
\end{aligned}\right.$ (27)

where $S_{n}$ is the size of the population utilized in the article’s proposed
evolutionary algorithm, $x^{j}$ for $j=1,\ldots,S_{n}$ are the elements of the
current population,
$K(x,x^{j})=\exp\bigl{(}-\|x-x^{j}\|_{2}/{(2\sigma^{2})}\bigr{)}$ is the
Gaussian kernel of $x$ and $x^{j}$, and $\sigma>0$ is the kernel radius. In
[47], the authors set $\sigma=2$. Here, the first objective focuses on
exploitation while the second focuses on exploration. There is no guarantee
that the two objectives $\varphi^{1}(x)$ and $\varphi^{2}(x)$ will conflict.

Article [97] reformulates a SNE as a dynamic tri-objective optimization
problem, and introduces a self-adaptive ranking multi-objective differential
evolution algorithm which is utilized to solve the dynamic tri-objective
optimization problem. The proposed tri-objective optimization problem is of
the form

$\left\\{\begin{aligned} &\min\varphi^{1}(x)\,\,\text{ for
}\,\,\varphi^{1}(x)=\frac{\sum_{j=1}^{S_{n}}\|x-r^{j}\|_{2}}{\max_{t=1,\ldots,P_{s}}\sum_{j=1}^{S_{n}}\|x^{t}-r^{j}\|_{2}}+\frac{9ng^{3}}{\gamma^{3}}\left(\sum_{i=1}^{m}|f_{i}(x)|+\max_{i=1,\ldots,m}|f_{i}(x)|\right),\\\\[8.5359pt]
&\min\varphi^{2}(x)\,\,\text{ for
}\,\,\varphi^{2}(x)=1-\frac{\sum_{j=1}^{S_{n}}\|x-r^{j}\|_{2}}{\max_{t=1,\ldots,P_{s}}\sum_{j=1}^{S_{n}}\|x^{t}-r^{j}\|_{2}}+\frac{9ng^{3}}{\gamma^{3}}\left(\sum_{i=1}^{m}|f_{i}(x)|+\max_{i=1,\ldots,m}|f_{i}(x)|\right),\\\\[8.5359pt]
&\min\varphi^{3}(x)\,\,\text{ for
}\,\,\varphi^{3}(x)=\sum_{j=1}^{S_{n}}\max\left(1-\frac{\|x-x^{j}\|}{\frac{\sqrt{n}}{\sqrt[n]{S_{n}}}(1-\frac{g}{\gamma})},0\right)+\frac{9ng^{3}}{\gamma^{3}}\left(\sum_{i=1}^{m}|f_{i}(x)|+\max_{i=1,\ldots,m}|f_{i}(x)|\right),\\\
\end{aligned}\right.$ (28)

where $S_{n}$ is the population size utilized in the proposed evolutionary
algorithm, $r^{j}$ for $j=1,2,\ldots,S_{n}$ are randomly generated reference
points at the initial stage, $x^{t}$ for $t=1,2,\ldots,S_{n}$ are the current
points in the population, $n$ is the dimension of the solution vector
$x\in\mathbb{R}^{n}$, $g$ is the index of the current evolutionary generation,
and $\gamma$ is the maximum number of generations over which the evolutionary
algorithm is allowed to evolve.

Since the proposed bi-objective and tri-objective optimization problems
introduced in Eq. (27) and Eq. (28) contain parameters which depend on the
configuration of algorithms proposed in the original papers, we refer the
interested reader to the source articles [97, 47] for additional details.

###  3.4 Transforming a SNE into a nonlinear programming problem

SNEs can also be explicitly reformulated as nonlinear programming problems.
Article [38] reformulates a square SNE (such that $m=n$) into a nonlinear
programming problem of the form:

$\min\,\bigl{\|}\varepsilon_{i}-\varepsilon_{i}^{*}\big{\|}_{\infty},$ (29)

subject to

$\left\\{\begin{aligned} \ \ &\\\\[-12.80365pt] \ \
&\bigl{|}f_{i}(x)\bigr{|}-\varepsilon_{i}\leqslant 0,\\\\[8.5359pt] \ \
&\varepsilon_{i}\geqslant 0,\\\\[1.42271pt] \end{aligned}\right.$ (30)

where $i=1,2,\ldots,n$ and $\varepsilon^{*}\in\mathbb{R}^{n}$ is the precision
of the desired solution (meaning that $\varepsilon_{i}^{*}=0$ corresponds to a
true solution to a SNE).

##  4 Optimization techniques to find solutions to a SNE

By formulating SNEs of interest as optimization problems, many techniques
traditionally utilized for finding one or more minima of a function can be
applied to attempt to find solutions to SNEs. There are many well known
optimization algorithms which can be used to attempt to find one or more
minima of an objective function, and the performance of these methods can vary
considerably depending on the characteristics of the objective function.
However, some optimization algorithms can only be applied to objective
functions with certain characteristics (of differentiability class $C^{2}$ for
example). Since some of the reformulations described in Section 3 can result
in non-smooth and multimodal objective function surfaces, it is important that
the characteristics of the reformulated objective function are considered when
a solution technique is selected.

In cases where the characteristics of the reformulated objective function are
difficult to determine, one should consider utilizing an optimization
algorithm which can be applied to non-smooth and multimodal objective
functions. If one can assume that the reformulated objective function is
unimodal and smooth (neither of which is true in the general case for SNEs),
one could simply use a local search method to attempt to find a solution to a
SNE via minimization. If one can only assume that the reformulated objective
function is unimodal, then a gradient-free local search technique should be
used. If one can not assume that the objective function is unimodal, then a
local search method alone is unlikely to be sufficient for finding a solution
to a SNE. In the multimodal case, one should attempt to find solutions to a
SNE by using an algorithm for global optimization, such as those described in
Section 4.2. Although the scope of our paper is focused on discussing
techniques which have been directly utilized in literature to search for
solutions to SNEs, there are many optimization algorithms which have yet to be
utilized to search for solutions to SNEs. In addition to learning about the
techniques explicitly described in this paper, we suggest that the interested
reader learn more about the taxonomy of global optimization algorithms by
reviewing articles such as [98].

In this paper, we have classified the optimization algorithms which have been
utilized to search for solutions to SNEs into two categories: (1) local search
methods, and (2) global search methods. Section 4.1 describes a few of the
many local search techniques which can be used to attempt to find a point
$x^{\prime}$ sufficiently close to a critical point on the objective function
surface. Therefore, a solution to a SNE can be found if a local search
algorithm terminates at a point $x^{\prime}$ which has an objective function
value sufficiently close to zero. Since the reformulated objective function of
an arbitrary SNE may be highly multimodal, local search techniques are often
utilized within hybrid metaheuristics and memetic algorithms to find solutions
to SNEs via global optimization. Section 4.2 focuses on introducing
metaheuristics which have been utilized to find solutions to SNEs via global
optimization. The book [99] provides a nice introduction to heuristics for
global optimization, and the books [100, 101] discuss strategies for
developing portfolios of such optimization algorithms. Unless otherwise noted,
in Section 4 we will describe techniques which can be used to solve SNEs
reformulated into single objective non-constrained global optimization
problems such as Eq. (9).

###  4.1 Local search methods

Since Newton’s method is capable of finding the roots of a differentiable
function, Newton’s method can be applied to the gradient $\nabla\varphi(x)$ of
a twice differentiable function $\varphi(x)$ to find critical points.
Therefore, if an objective function such as Eq. (10) is twice differentiable,
Newton’s method can be applied from an arbitrary starting point
$x^{0}\in\mathcal{D}\subset\mathbb{R}^{n}$ to find a nearby local minimum or
saddle point on the objective function surface. Once Newton’s method
terminates at a point $x^{\prime}\in\mathcal{D}\subset\mathbb{R}^{n}$, one can
then check whether or not the objective function value at $x^{\prime}$ is
sufficiently close to zero to determine whether or not $x^{\prime}$ is a
solution to the NSE of interest (since every solution to a SNE corresponds to
a point with an objective function value of zero).

If the objective function is differentiable, then nonlinear gradient methods
and nonlinear conjugate gradient methods may be used to attempt to find one if
its critical points. Gradient descent methods are based on the observation
that the function $\varphi(x)$ decreases the most in the direction opposite to
its gradient $\nabla\varphi(x)$. Hence, starting from an initial solution
$x^{0}\in\mathcal{D}\subset\mathbb{R}^{n}$, the gradient descent method
iteratively finds a sequence of new solutions by calculating
$x^{i+1}=x^{i}-\alpha^{i}s^{i}$ where $s^{i}=\nabla\varphi(x^{i})$ is the
search direction and $\alpha^{i}$ is the step size. There are multiple methods
for determining the step size $\alpha^{i}$, including the well known Barzilai-
Borwein method [102]. Article [103] provides a nice introduction to the
Barzilai-Borwein method and compares the Barzilai-Borwein method against
alternative methods. Overview of gradient descent optimization methods may be
found in [104]. Popular extension of the gradient descent method is called
stochastic gradient descent, that can be efficient when the function
$\varphi(x)$ can be represented as a sum of large number of functions
$\varphi(x)=\sum_{i}^{K}\varphi_{i}(x)$. Then, the gradient is estimated on a
random subset of $\varphi_{i}(x)$. Widely used stochastic gradient descent
algorithms are ADAM [105], AdaGrad[106], and others.

Conjugate gradient methods differ from gradient descent methods by using the
conjugate search direction $s^{i}=-\nabla\varphi(x^{i})+\beta^{i}s^{i-1}$.
Multiple methods for calculating the parameter $\beta^{i}$ of the conjugate
search direction have been proposed in the literature. Commonly utilized
methods for calculating $\beta^{i}$ were introduced by Fletcher–Reeves [107],
Polak–Ribière [108], Hestenes-Stiefel [109], and Dai-Yuan [110] among others.
Many variants based on conjugate gradient methods have been proposed in
literature. For example, articles [16, 111, 55, 112] propose new methods which
combine hyperplane projection techniques with the conjugate gradient method to
solve systems of monotone nonlinear equations (SNEs which satisfy the
inequality in Eq. (4)). Article [113] applies the conjugate gradient method to
solve SNEs, and it suggests combining the Fletcher–Reeves and Polak–Ribière
methods with the Quasi-Newton update. The paper first reformulates a SNE into
a unconstrained single objective optimization problem, and then searches for a
solution. Article [113] numerically evaluates the performance of the proposed
method on 20 test problems with up to 10k dimensions, and most of the problems
were solved by the proposed method in under 3 seconds and with less than 30
iterations of the algorithm when coded in MATLAB and run on a computer with a
2.13 GHz CPU processor and 2GB of RAM. Additional works which discuss using
conjugate gradient based methods for solving SNEs reformulated as single
objective optimization problems include [114, 115].

In addition to Quasi-Newton methods, other commonly utilized derivative-free
direct search methods worth mentioning include the Rosenbrock search method
[116], the Hooke and Jeeves pattern search method [117], Brent’s method [118],
and the Nelder-Mead simplex method [119]. The Nelder-Mead simplex method in
particular has been frequently utilized as a local search method within hybrid
metaheuristics and memetic algorithms for finding solutions to SNEs via global
optimization [37, 92, 120, 2]. The books [121, 122] provide a nice
introduction to methods for derivative-free optimization.

###  4.2 Global search / global optimization methods

Global optimization is a well developed field, and There are many ways to
classify global search methods as shown in article [98], but in this paper we
decide to classify them into three main categories: (1) exact methods, (2)
single-solution-based methods, and (3) population-based methods. We consider
exact methods to be methods which are guaranteed to find a solution via
exhaustive search, we consider single-solution-based methods to be methods
which search the solution space with a single agent at a time
($\text{population}=1$), and we consider population-based methods to be
methods which search the solution space simultaneously with many agents
($\text{population}\gg 1$). Since much of the recent literature discussing
solution techniques for SNEs has focused on metaheuristcs, we focus our
attention on the two main classes of metaheuristics: the single-solution-based
and population-based methods. Hybrid methods, or methods which combine
metaheuristics with local search algorithms, are not explicitly separated from
the single-solution based methods and the population-based methods in this
survey as many of the metaheuristics mentioned in these sections can easily be
hybridized if desired. For example, see article [37] which hybridizes the
Continuous Variable Neighborhood Search (C-VNS) metaheuristic with the Nelder-
Mead method for local search. For readers who are interested in learning more
about global optimization methods and applications, we refer the interested
reader to the books [123, 124, 125, 126, 127].

####  4.2.1 Exact methods

Exact methods for global optimization guarantee that a solution will be found
within a specific tolerance (if one exists) by exhaustively searching the
solution space. Exact methods for global optimization are typically too
computationally expensive to be applied in practice, especially in cases where
the objective function is highly multimodal or the number of dimensions is
very large. The Branch-and-bound algorithm is a well known exact method for
global optimization [128].

####  4.2.2 Single-solution-based methods

Single-solution-based methods explore the solution space with a single agent
who makes decisions about which path to take based upon a predefined set of
rules typically based upon the objective function values at the points which
the agent has explored so far. This section introduces single-solution-based
methods which have been utilized in literature to search for solutions to SNEs
through global optimization.

Continuous Greedy Randomized Adaptive Search Procedure (C-GRASP): In [129],
Hirsch et al. introduced the Continuous Greedy Randomized Adaptive Search
Procedure (C-GRASP) which extends Feo and Resende’s Greedy Randomized Adaptive
Search Procedure (GRASP) [130] for discrete optimization to the domain of
continuous global optimization. The C-GRASP metaheuristic utilizes multistart
and a derivative-free stochastic local search method to perform global
optimization. In [88], Hirsch et al. utilize an enhanced version of the
original C-GRASP metaheuristic proposed in [129] to find multiple solutions to
Eq. (9) via global optimization. Each iteration of C-GRASP is initialized with
an initial coarse discretization of the search space and an initial random
solution vector $x$. Each iteration of C-GRASP consists of a series of
construction-local improvement cycles where the output of the construction
phase is used as the input to the local search phase, and where the output
from the local search phase is used as the input for the next iteration’s
construction phase. As the iterations progress, the discretization of the
search space becomes increasingly dense as C-GRASP converges to a minimum.
After reaching a maximum number of multi-start iterations, the C-GRASP
metaheuristic returns all the solutions to the SNE that it found.

During C-GRASP’s construction phases, C-GRASP utilizes greediness and
randomization to create a diverse set of quality initial solutions from which
to start local search. Each construction phase takes an initial feasible
solution $x$ as an input, and starts by letting all coordinates of $x$ be
unfixed. Then, in an iterative fashion, C-GRASP’s construction phase will
first conduct a line search in each unfixed coordinate direction of $x$ (while
holding all other $n-1$ coordinates at their current values) to find the
coordinate value along each dimension (using the current discretization) that
results in the best objective function value. After conducting this line
search, the best and worst objective function values found across all
coordinates are saved in the variables $min$ and $max$. Then, each unfixed
coordinate of $x$ is set to the best value found during the coordinate-wise
line search. From here, a restricted candidate list of unfixed coordinates is
formed which is comprised of the coordinates whose best objective function
values found during line search were less than or equal to
$\alpha\max+(1-\alpha)\min$ where $\alpha\in[0,1]$ is randomly selected at the
beginning of each construction phase. Then, a dimension in the restricted
candidate list is selected at random, and the corresponding coordinate value
of $x$ is fixed and will no longer be moved throughout the current
construction phase. By selecting a coordinate at random from the restricted
candidate list (meaning the not necessarily the best coordinate is selected to
be fixed) the authors preserve some randomness which they argue helps the
C-GRASP metaheuristic explore more of the search space. This process repeats
until all coordinates of $x$ are fixed which ends the construction phase. The
final feasible point $x$ found by the construction phase will then serve as
the starting point for the next derivative-free local search phase. Here, even
though $x$ may be a good solution, it may not be optimal as the line search
was performed along each coordinate direction. As a result, local search is
then conducted starting at the feasible point $x$ found by the construction
phase.

During each derivative-free local search phase, C-GRASP searches in the
neighborhood of the current solution for a better solution, where the
neighborhood of the current solution $x$ is defined as the projection of all
the current grid points onto the hyper-sphere centered at $x$ with radius
equal to the current discretization level $h$. C-GRASP then determines at
which points in the neighborhood of $x$, if any, the objective function
improves. Then, if an improving point is identified, it is made the new
current point and local search continues from this new point. This local
search process continues until an approximate local minimum $x^{\prime}$ is
reached such that taking a step in any evaluated direction would result in a
worse objective function value. The $x^{\prime}$ returned by each local
improvement phase then serves as the starting point for the next construction
phase. Specifically, the first construction phase (before any local search has
been conducted) starts at a random feasible point, every construction phase
after the first starts at the approximate local minimum $x^{\prime}$ found by
the previous local improvement phase, and each local improvement phase starts
at the point generated by the previous construction phase.

The original version of C-GRASP [129] was evaluated on three SNEs: (a) a SNE
arising in robot kinematics that is comprised of eight nonlinear equations in
eight unknowns and which has 16 known roots [131, 124], (b) a SNE from the
area of chemical equilibria which has ten equations and ten unknowns and which
has one physical solution where all variables are positive [132], and (c) a
transformed version of the original SNE from the area of chemical equilibria
which has five equations, fine unknowns, and one solution where all variables
are positive [132]. The original version of C-GRASP was tested ten times on
the SNE arising in robot kinematics, and successfully found all 16 solutions
each time with an average runtime of 3,048s. On the original version of the
SNE arising in chemical interaction and equilibrium modeling, C-GRASP
successfully found the solution on eight out of ten runs with an average run
time of 201.58s. On the transformed version of the SNE arising in chemical
interaction and equilibrium modeling, C-GRASP successfully found the solution
on all ten runs with an average runtime of 37.53s. The improved version of
C-GRASP [88] was evaluated on four SNEs, three of which had two variables and
two unknowns, and one of which had three solutions and three unknowns. Over
thousands of trials, C-GRASP found all solutions to each system almost every
time, and the longest average runtime for any of the problems tested was
5.06s. Other systems tested were solved in fractions of a second.

Continuous Variable Neighborhood Search (C-VNS): Continuous Variable
Neighborhood Search (C-VNS) is a continuous extension of the variable
neighborhood search (VNS) metaheuristic proposed by [133], a popular method
for solving discrete optimization problems. C-VNS systematically explores the
solution space by first defining a set of neighborhood structures, and then
repeatedly performing local search starting from a point within the
corresponding neighborhoods which increase in size when improving solutions
are not found. C-VNS was utilized to search for solutions to SNEs in [37]. The
paper reformulates a SNE as an optimization problem and performs minimization
of the resulting objective function using C-VNS. The paper defines the set of
neighborhood structures to be a set of hyperspherical shells with increasing
internal and external radii, so that the union of the shells covers the entire
solution space. C-VNS then performs local search using the Nelder-Mead method,
starting at a random point in the current neighborhood structure. The paper
experimentally evaluates multiple local search methods, along with multiple
ways to reformulate problem 1. The paper concentrates on finding all solutions
of small SNEs (2-20), and the largest SNE described in the paper consisted of
1000 equations.

Tabu Search (TS): Developed initially for solving combinatorial optimization
problems, Tabu Search (TS) is a metaheuristic approach that can be utilized
for global optimization. Starting at an initial feasible solution $x$, TS
follows a pre-defined set of rules to identify a set of ${\mathcal{Y}}$
solutions in the neighborhood of $x$ to explore. Let ${\mathcal{N}}(x)$ denote
the neighborhood of $x$ such that ${\mathcal{Y}}\subset{\mathcal{N}}(x)$.
After evaluating the objective function value at each point
$y\in{\mathcal{Y}}$, the point $y^{*}\in{\mathcal{Y}}$ with the best objective
function value is selected as the new starting point for the next iteration of
TS, even if the initial point $x$ has a better objective function value than
$y^{*}$. By allowing moves to points with worse objective function values, TS
attempts to escape the neighborhood of previously found local minima. While
exploring a discrete or discretized solution space, TS stores a list of
recently visited points called a “tabu list”. To ensure that TS continues to
explore new regions of the solution space, TS prevents future moves to points
in the tabu list. Other memory structures containing information such as the
best points visited so far may also be maintained and utilized in an attempt
to facilitate the performance of TS [134]. TS continues to explore the
solution space until a stopping criteria is met (in the case of SNEs, when the
function value is sufficiently close to zero or upon a maximum number of
iterations).

One variant of TS, named, Directed Tabu Search (DTS), incorporates local
search methods into the TS search procedure. Local search can either be
conducted from the identified $y^{*}\in{\mathcal{Y}}\subset{\mathcal{N}}(x)$
to find a local minimizer $x^{*}$ from which to start the next iteration of
DTS, or local search can be performed immediately from an initial starting
point $x$ to find a local minimizer $x^{*}$ from which to select and evaluate
the points $y\in{\mathcal{Y}}\subset{\mathcal{N}}(x^{*})$ as described in
Algorithm 2.1 of [134].

Numerical results presented in [134] showed that the described DTS method was
tested on SNEs with up to $n=51$ decision variables, and that it was
successful in converging on some but not all of the SNEs within $100\,n^{2}$
function evaluations.

Simulated Annealing (SA): Originally proposed in [135] for solving hard
combinatorial optimization problems, Simulated Annealing (SA) is a well known
metaheuristic that can be utilized to find approximate solutions to global
optimization problems. SA only utilizes evaluations of the objective function
during its search procedure, and can therefore be utilized in cases where
information about the derivatives of the objective function is inaccessible.
Although it is not always guaranteed that SA will converge to a true global
minimum (some papers such as [136] have proved the convergence of SA under
certain specific conditions), SA can be used in cases where the objective
function is non-smooth or even discontinuous over the feasible region [137].

Starting with in initial step vector $v^{0}\in\mathbb{R}^{n}$, an initial
parameter $T^{0}\geqslant 0$ called the “temperature”, and an initial feasible
solution $x^{0}\in\mathcal{D}\subset\mathbb{R}^{n}$, SA attempts to generate a
sequence of successive points $x^{k}\in\mathcal{D}\subset\mathbb{R}^{n}$ for
$k=1,2,\ldots,K$ which approach a global minimum $x^{*}$ of the objective
function. At each iteration of SA, a new candidate solution $x^{\prime}$ is
generated from the previous solution $x^{k}$ by taking random moves along each
coordinate direction in turn. Each new coordinate value $x^{\prime}_{i}$ for
$i=1,2,\ldots,n$ is sampled uniformly from an interval of width $2v^{k}_{i}$
centered around the previous coordinate value $x^{k}_{i}$. If any new
coordinate $x^{\prime}_{i}$ makes $x^{\prime}$ fall outside of the feasible
region $\mathcal{D}$, then a new coordinate value of $x^{\prime}_{i}$ is
randomly generated until a feasible point is found or a maximum number of
attempts is met. When a new feasible point $x^{\prime}$ is found, if
$\varphi(x^{\prime})\leqslant\varphi(x^{k})$, we set $x^{k+1}=x^{\prime}$ so
that $x^{\prime}$ is the starting point for the next iteration of SA.
Alternatively, if $\varphi(x^{\prime})>\varphi(x^{k})$, SA sets
$x^{k+1}=x^{\prime}$ with some probability $P(\text{acceptance})$. For
example, $P(\text{acceptance})$ can be defined as
$P(\text{acceptance})=\exp(-(\varphi(x^{\prime})-\varphi(x^{k}))/T^{k})$ where
$T^{k}$ is the “temperature” parameter. Notice that if $T^{k}=0$, only
improving points will be selected, but as $T^{k}$ increases, the probability
that non-improving moves will be accepted increases. SA accepts non-improving
moves to allow the algorithm to avoid being trapped in local minima. A new
step vector $v^{k+1}\in\mathbb{R}^{n}$ and a new “temperature”
$T^{k+1}\geqslant 0$ can be selected at each iteration as desired. Typically,
the “temperature” $T$ and the elements of the step vector $v$ are gradually
decreased as the number of iterations increases so that the algorithm can take
smaller and smaller improving steps as it (hopefully) approaches a global
minimum.

When tested against the Nelder-Mead simplex method on a variety of global
optimization problems, SA was found to take 500 to 1000 times more function
evaluations than the Nelder-Mead simplex method, but it more reliably found a
global minimum of the objective functions it was tested on than the Nelder-
Mead method [137]. Expanding upon the work presented in [137], article [5]
introduces a metaheuristic which utilizes SA to solve the objective function
described in Eq. (13) which has a penalty designed to repel SA from already
found solutions. Article [5] refers to this penalized objective function as a
“polarization technique”. Article [138] attempts to find multiple solutions to
SNEs by repeatedly solving Eq. (9) with a multistart fuzzy adaptive SA
algorithm. The proposed algorithm successfully found multiple solutions to
SNEs of up to 20 equations and 20 unknowns ($m=n=20$) in less than four
minutes when programmed in C++ and run on a computer with 2 CPU at 2.4 GHz and
512 MB of RAM.

####  4.2.3 Population-based methods

Population-based methods explore the solution space with many agents
simultaneously, typically in a coordinated manner which balances exploration
and exploitation. This section introduces population-based methods which have
been utilized in literature to search for solutions to SNEs via global
optimization.

Evolutionary Algorithms: Evolutionary algorithms are a class of metaheuristic
optimization algorithms, principles of which are inspired by mechanisms of
biological evolution. Often, evolutionary algorithms associate candidate
solutions to biological individuals, which gradually evolve, and only the
fittest survive. For example, article [36] proposes an evolutionary algorithm
for solving SNEs which first ranks all candidate solutions of the current
population by their fitness values, and then performs mutation operations on
the best candidate solutions to generate “offspring” which will serve as the
new set of candidate solutions for the next iteration of the evolutionary
algorithm. Here, fitness is quantified as a mathematical fitness function. The
set of candidate solutions forms a population, that is being iteratively
improved. An attractive property of some evolutionary algorithms is niching
[139], which can help find multiple solution to a SNE by dividing the whole
population into multiple diverse subpopulations [43]. We refer the interested
reader to the book [140] which discusses utilizing evolutionary algorithms for
multimodal optimization.

A large number of research papers are devoted to the application of
evolutionary algorithms for solving SNEs transformed into multiobjective
optimization problems (Eq. (21)). Some of the approaches used for solving
multiobjective optimization problems involve “scalarization”, that is,
conversion to a problem of the form similar to Eq. (8). Evolutionary
algorithms are popular as they may generate multiple solutions at the same
time, which may be used to model the Pareto-front of the optimization problem.
While [18] proposed to use each equation as an objective, thus forming a many-
objective optimization problem, paper [40] proposed to utilize a
multiobjective evolutionary algorithm called MONES to search for solutions to
a SNE reformulated as Eq. (24) where all optimal solutions of the SNE would be
Pareto optimal solutions of Eq. (24).

Genetic Algorithms (GAs): Genetic algorithms (GAs) are one of the most
commonly utilized types of evolutionary algorithms. Memetic algorithms (also
known as genetic local search algorithms) are a class of GAs which utilize a
local search technique to attempt to improve the quality of the approximate
solution(s) found. Article [43] discusses utilizing memetic algorithms for
searching for solutions to SNEs. Paper [89] applied the Random-Key Genetic
Algorithm (RKGA) for finding solutions of a box-constrained SNE, i.e. such
that their solutions $x^{*}\in S$, where
$S=\\{x=(x_{1},x_{2},\ldots,x_{n})\in{\mathbb{R}}^{n}:l_{i}\leqslant
x_{i}\leqslant u_{i}\\}$. The RKGA was originally proposed for solving
discrete optimization problems [141]. At each iteration of the evolutionary
process, RKGA computes the value of the fitness function of its chromosomes;
after which, the most fit chromosomes are placed into an “elite” group, and
remaining chromosomes are considered to be “non-elite”. Then, the next
population is created by combining all the “elite” chromosomes, a small group
of new random chromosomes, and another small group of new chromosomes which
are generated by performing crossover between “elite” and “non-elite”
chromosomes. There are many different variations of the RKGA, some of which
which differ from the original algorithm by allowing chromosomes to have more
than one offspring, and others by modifying the crossover operation so that
“parents” are taken from the entire population. The RKGA also contains an
encoder and a decoder which encode and decode the solution of the optimization
problem to the chromosome representation. In [89], the encoder represents a
solution as a random vector with each element $\chi_{i}$ lying in the interval
[0,1], and the decoder “rescales” this random vector by calculating
$x_{i}=l_{i}+\chi_{i}\,(u_{i}-l_{i})$. In addition to the original RKGA
procedure, article [89] proposes to utilize a local improvement procedure to
attempt to find a better solution in the neighborhood of the current solution.
The gradient-free local search method utilized in [89] is based on the method
introduced in [142]. Article [89] evaluates the method on small-scale problems
from kinematics, chemical engineering, and other fields [51]. Article [34]
utilizes a GA to solve a SNE transformed into a constrained optimization
problem. The constraints are handled by adding a penalty into the objective
function which penalizes solutions for which the constraints are not
satisfied.

GAs have also been utilized within hybrid algorithms for solving SNEs. For
example, article [32] attempted to solve SNEs reformulated as a single
objective optimization problem (Eq. (9)) by using the novel hybrid-GOA-GA
algorithm which combines the grasshopper optimization algorithm (GOA) with a
GA. The hybrid-GOA-GA algorithm was tested on eight benchmark SNEs of up to
ten equations and unknowns, and on the largest SNEs the hybrid-GOA-GA
algorithm converged to a solution in less than 15 seconds on average.

Differential Evolution (DE): Similarly to GAs, algorithms based on
differential evolution (DE) work by generating a population of candidate
solutions which get iteratively improved by combining candidate solutions from
the previous population. New candidate solutions are generated by subsequently
applying mutation and crossover operations to solutions from the previous
generation; then, the solutions are compared with their parents and are
selected and added to the next population if their fitness function value is
better.

Article [143] proposes to hybridize the Grey Wolf Optimization (GWO) algorithm
with a DE algorithm; the algorithm is utilized to search for solutions to a
SNE reformulated as an unconstrained optimization problem. The algorithm was
evaluated on small-scale SNEs and was able to find multiple solutions to SNEs
without using an objective function penalty (Eq. 10) since the population of
candidate solutions can potentially contain multiple solutions to a SNE found
simultaneously in a single run. Article [48] proposes to explore new candidate
solutions using a fuzzy neighborhood technique (based on Euclidean distance).
Paper [46] also reformulates a SNE as a global optimization problem, and
proposes a repulsion-based adaptive DE algorithm which utilizes additive and
multiplicative repulsion techniques (cf equations 9 and 14) along with
adaptive control of DE parameters. Article [49] utilizes niching techniques
with DE in order to search for multiple solutions to a SNE simultaneously in a
single run. It has also been proposed to apply k-means clustering to the set
of candidate solutions, and to treat each cluster as a separate population to
preserve the diversity of the solutions [45]. Article [47] reformulates a SNE
as a multiobjective optimization problem and attempts to solve it using a two-
phase EA. The first phase attempts to maintain the population diversity by
combining a multiobjective optimization technique with a niching technique,
and the second phase attempts to detect promising subregions within the
solution space and explore them by using DE to perform local search. Niching
is also used with DE in [52]. Article [97] applies DE with modified crossover
to a SNE transformed into a tri-objective multiobjective optimization problem.

Particle Swarm Optimization (PSO): Particle Swarm Optimization (PSO)
algorithms generate a set of candidate solutions (“particles”) that are moving
in the search space in accordance to some movement law. Typically, each
particle has a velocity that controls its movement [144]. Particles’ movement
directions are iteratively updated based on the best solution of the
individual particle, and the current best solution of all particles. Article
[120] reformulates a SNE as an optimization problem of the form Eq. (8) and
searches for solutions to it using a hybrid PSO algorithm. The hybrid PSO
algorithm utilized in article [120] is a combination of PSO with the Nelder-
Mead local search method. Here, some of the particles update their positions
according to movement equations, and others update their positions using the
Nelder-Mead method. The paper evaluates the method on small SNEs with less
than seven equations, and looks for only one solution. Paper [86] proposes a
hybrid algorithm combining PSO and Cuckoo search for solving SNEs. The
reformulation utilized in article [86] is the square root of the objective
function shown in Eq. (8) with $\alpha=1$ and $p=2$. Article [145] utilizes
PSO to minimize an objective function of the form shown in Eq. (8) where
$\alpha=1$ and $p=1$. The algorithm utilized in [145] is a modified PSO
algorithm with a new expression for finding the inertia weight. The algorithm
is used to find only one solution, and is tested on small SNEs. Paper [146]
proposes utilizing Quantum behaved Particle Swarm Optimization (QPSO) to solve
a SNE reformulated as Eq. (9). The main difference between the QPSO algorithm
proposed in [147] and the traditional PSO algorithm is that the QPSO algorithm
utilizes principles from quantum mechanics while the traditional PSO algorithm
utilizes principles from Newtonian mechanics. Additionally, the method
presented in article [147] changes the pseudo-random number generator to a
chaotic map. The resulting algorithm is evaluated on small-size problems, and
their experiments show that the proposed method tends to converge to solutions
faster than the ordinary PSO algorithm. Article [148] aims at solving a SNE by
minimizing an objective function of the form shown in 8 with $\alpha=1$ and
$p=2$. The proposed method combines the standard PSO algorithm with Powell’s
conjugate direction (CD) method which can be used to find a local minimum of a
function. The resulting method is used to find one solution and is tested on
SNEs of 5-10 dimensions. PSO was also utilized to search for solutions to SNEs
in [149, 150]. MATLAB code for PSO and a modification of PSO named “Unified
PSO” that aggregates its local and global variants [151] can be found in [152,
144]. PSO has also been hybridized with insights from memetic algorithms as
discussed in article [153].

Spiral Optimization (SPO): Initally introduced in articles [154, 155], spiral
optimization (SPO) methods are multistart methods for global optimization
which are inspired by spiral phenomena in nature. Article [155] initially
proposed a SPO method for solving 2-dimensional continuous optimization
problems, and article [154] expanded upon this work to introduce a SPO
algorithm for $n-$dimensional continuous optimization problems. The
$n-$dimensional SPO algorithm is constructed based on rotation matrices
defined in the $n-$dimensional space. The iterative SPO algorithm explores the
solution space with a population of solutions simultaneously, each of which
explores the solution space by moving in a tightening spiral centered around
the best solution found so far. Upon initialization, the spiral dynamics
optimization algorithm randomly selects a set of starting points within the
domain of interest $\mathcal{D}$, and evaluates the objective function value
at each point. The point with the best objective function value is set as the
initial center point which all other members of the population will spiral
around. The motivation behind the spiral model is that initially when the
spiral is large, the model will explore a diverse set of points across the
solution space. As the number of iterations increases the spiral tightens, and
the search intensifies around the center of the spiral to facilitate
exploitation. One interesting aspect of the SPO algorithm is that it explores
the solution space with a population of candidate solutions which move in a
non-stochastic manner.

Article [87] presents a method which utilizes SPO and clustering to search for
multiple solutions to a SNE. Instead of selecting the set of starting points
for SPO in a purely random fashion, article [87] utilizes the low-discrepancy
Sobol sequence to attempt to select a set of starting points which are more
likely to be uniformly spread across $\mathcal{D}$. By selecting starting
points which are more uniformly spread across $\mathcal{D}$, one can
facilitate the exploration of the model by powerfully leveraging the
advantages presented by SPO’s multistart framework. We refer the interested
reader to Chapter 2 of the book [156] for information on low-discrepancy
sequences and on how to generate random numbers with specified distributions.

To attempt to find multiple roots to a SNE in a single run, article [87]
proposes a method which utilizes a clustering technique to identify different
regions which are likely to contain a solution to a SNE, and simultaneously
utilizes SPO to search for solutions within each cluster. The proposed method
was tested on six SNEs of up to eight equations and eight unknowns, and the
method demonstrated the ability to effectively find multiple solutions to each
SNE in a single run in a matter of seconds.

Other metaphor-based metaheuristics: Cuckoo search models solutions of an
optimization problem as eggs in a nest. Initially, there are multiple nests.
The algorithm performs multiple iterations; at each iteration a randomly
chosen “egg” (a solution) in a random “nest” is replaced by another solution
(modeling a situation where a cuckoo lays her egg in another host bird’s
nest); then the fitness function of the “nests” of solutions is evaluated, and
the best nests are moved to the next generation, while a randomly chosen
fraction of nests of worse quality is dropped (modeling the situation when a
host bird discovers the cuckoo’s egg and leaves the nest). Randomization is
often performed using Levy flights random walk. Cuckoo search was applied to
SNEs in [157, 86, 158].

Monarch butterfly optimization (MBO) algorithms, proposed in [159], are
inspired by the behavior of North American Monarch Butterflies, and their
autumn migration from Northern USA and Southern Canada to Mexico. In MBO, a
population of Monarch butterflies is split into two sub-populations (modeling
populations of butterflies in “Land 1” and “Land 2”). MBO then performs two
operations, one of which is a migration operator which creates new butterflies
which potentially replace butterflies from the previous population, and the
other of which is a butterfly adjusting operator which changes the positions
of the butterflies. These operations are iteratively applied to the two sub-
populations until a stopping criteria is met. MBO is used to search for
solutions to SNEs in [160].

Many other metaphor-based metaheuristics have been used to solve SNEs,
including invasive weed optimization [33, 90], glowworm swarm optimization
[161], hybrid artificial bee colony algorithm [162], imperialist competitive
algorithms [163], modified firefly algorithm [164], harmony search [165, 91],
and a soccer league competition algorithm [166] among others.

Chaotic global optimization methods: Chaotic global optimization methods
typically generate one or more chaotic sequences, and then map these chaotic
sequences to a corresponding sequence of $P$ feasible solutions
$x^{p}=(x_{1}^{p},x_{2}^{p},\ldots,x_{n}^{p})\in\mathcal{D}\subset\mathbb{R}^{n}$
for $p=1,2,\ldots,P$ to the global optimization problem of concern. The set of
generated feasible solutions $x^{p}$ are then used in some way, for example,
they could be utilized as starting points for local search within a multi-
start framework as shown in [167]. The reason why chaotic sequences are
utilized as the basis for the set of starting points $x^{p}$ is because
chaotic sequences have three key traits, namely, (a) pseudo-randomness, (b)
ergodicity, and (c) irregularity [168].

Article [31] presented a two-stage metaheuristic to find solutions to a SNE by
solving Eq. (9) where $\varphi^{0}(x)$ is of the form shown in Eq. (8). The
authors of [31] claim that the first stage of the metaheuristic attempts to
find a hyper-rectangle $S_{2}=I_{S_{2_{1}}}\times
I_{S_{2_{2}}}\times\cdots\times I_{S_{2_{n}}}$ within the feasible region over
which $\varphi^{0}(x)$ is unimodal. Starting with a larger feasible hyper-
rectangle $S_{1}=I_{S_{1_{1}}}\times I_{S_{1_{2}}}\times\cdots\times
I_{S_{1_{n}}}$ such that $S_{2}\subseteq S_{1}$ (meaning that $\varphi^{0}(x)$
is not necessarily unimodal over $S_{1}$), the authors utilize the elements of
a chaotic sequence in a process which iteratively shrinks $S_{1}$ for a set
number of iterations to obtain $S_{2}$. It is important to note that though
they hope to obtain a hyper-rectangle $S_{2}$ over which $\varphi^{0}(x)$ is
unimodal, there is no guarantee that this will occur. After finding $S_{2}$,
the second stage of the metaheuristic presented in [31] utilizes an
$n$-dimensional expansion of golden section search [169] to perform local
search over $S_{2}$. Since $n$-dimensional golden section search is designed
to be utilized on a hyper-rectangle over which the objective function is
unimodal, stage one attempts to find this region. The authors of [31] tested
their algorithm on a variety of SNEs of up to 10 variables and 10 unknowns,
and were able to solve each problem in less than a second. However, they did
not compare their algorithm to any other techniques capable of finding
multiple solutions to a SNE, and there is little evidence to suggest that
stage 1 of their algorithm is capable of consistently finding a hyper-
rectangle $S_{2}$ such that $\varphi^{0}(x)$ is unimodal (as is required to
guarantee that $n$-dimensional golden section search will converge to a global
minimum of $\varphi^{0}(x)$ over $S_{2}$).

###  4.3 Neural network-based optimization

Various neural network-based optimization techniques have been utilized to
search for solutions to SNEs. These approaches construct and train an
artificial neural network (ANN) for each SNE; then, roots of the SNE may be
found at the neurons’ weights. An example of this technique is discussed in
[170], where it is referred to as a Gradient Neural Network and used to search
for the solution to system of linear equations. Similar approaches are
presented in [171] and [172]. In these approaches, solutions to the equations
also correspond to the state vector of the ANN. These neural network based
approaches are different from traditional modern application of ANNs, which
are expected to be able to find solutions to unseen SNEs, after being trained
on many SNEs such as [173]. In [174], Margaris et al. presented a back-
propagation trained four-layered feed-forward ANN framework for solving
algebraic systems of $n$ polynomial equations with $n$ unknown variables. The
four-layered network presented in [174] consists of:

  * (a)

a single neuron input layer whose output synaptic weights are trained to
become the components of one of the SNE’s solutions,

  * (b)

a first hidden layer of $n$ summation neurons which generate each of the
linear terms $\\{x_{1},x_{2},\ldots,x_{n}\\}$ and which have an identity
activation function,

  * (c)

a second hidden layer of product units with custom activation functions which
generate all of the nonlinear higher order terms of the algebraic system by
multiplying the linear quantities received from the first hidden layer (for
example, if a SNE had the nonlinear term $x_{1}^{2}$, a corresponding neuron
will have an activation function that squares the input $x_{1}$),

  * (d)

an output layer of $n$ summation neurons that generate the right hand side of
each nonlinear equation $f_{i}(x)=0$ of the SNE (serving as the estimate of
the training error when compared to each equation’s true right hand side value
of 0),

where the only synaptic weights adjusted during training are those between the
single neuron in the input layer and the neurons of the first hidden layer
which represent each unknown variable of the system. Expanding upon this work,
Awaad et al. utilized a new Manhattan-Adam hybrid updating rule to accelerate
model convergence and overcome the vanishing gradient problem [175]. The main
goal of the paper is to be able to use hardware accelerators that provide
efficient implementation of this update rule. The new updating rule
dynamically adjusts the step size during training, and rotates after a set
number of iterations between using the Manhattan updating rule and the Adam
updating rule to efficiently handle problems with sparse and vanishing
gradients.

Also expanding upon the work of Margaris et al. [174], Goulianas et al.
introduced the General Back-propagation with Adapative Learning Rate (GBALR)
algorithm which allows the activation functions of the second hidden layer to
be any function, including non-algebraic functions [176].

A related technique is called a Zeroing Neural Network (ZNN) [177]. ZNNs aim
at zeroing error functions, the design of which depends on the problem. In the
case of SNEs, similar to the research described above, the goal of ZNNs is to
find a solution to a nonlinear system by applying gradient-based optimization
methods. However, in contrast to previous papers which solve static systems, a
ZNN introduces temporal dependency into its error function; that is,
coefficients of SNEs may be dependent on an additional time parameter $t$.
Also addressing temporal SNEs with ANNs, Zhang utilized a Recurrent Neural
Network (RNN) to find approximate online solutions to a time-variant SNE
emerging in robot kinematics [11].

##  5 Synopsis and concluding remarks

In this article, we have presented part two of a survey on methods for solving
a system of nonlinear equations (SNE). In part two we have presented a
comprehensive survey of techniques which can be utilized to search for
solutions to a SNE via optimization. Since many of the SNEs that arise in real
world applications are considered over a finite bounded domain $\mathcal{D}$,
we first introduced a technique which can be utilized to determine the number
of solutions to a SNE that exist within $\mathcal{D}$. Then, we introduced
various techniques which can be used to transform a SNE into an optimization
problem, and we described a vast array of optimization algorithms which can
then be used to search for solutions to a SNE after transformation. We
discussed both local search algorithms and global search algorithms for
optimization, and we place an emphasis on cutting edge metaheuristics for
global optimization which have been shown to be particularly effective at
finding one or many solutions to SNEs.

Analyzing this literature has led us to conclude that for the hardest classes
of SNEs, many of the most promising techniques that we currently have are
metaheuristics such as those discussed in Section 4.2 of this paper. Although
such techniques are not necessarily guaranteed to converge to a solution in
finite time, many can be applied to ill-formed and challenging classes of SNEs
as they only require evaluations of the objective function (and hence can be
applied when the objective function is not everywhere differentiable).
However, more research needs to be done towards improving these techniques. To
address this need, we are currently exploring how machine learning can be
utilized to enhance the effectiveness of metaheuristics for solving SNEs
through global optimization. Additionally, as described in part one of this
survey, we are also actively developing a new taxonomy of SNEs to facilitate
the identification of new classes of tractable problems as well as the methods
most capable of solving them. Furthermore, we are actively exploring methods
capable of solving systems of nonlinear equations and inequalities.

In part one of this survey, we discussed methods for solving SNEs without
transforming them into optimization problems. In part two of this survey, we
discussed methods for solving SNEs by transforming them into optimization
problems. In part three of this survey, we will present a robust quantitative
comparative analysis of methods capable of searching for solutions to SNEs.

## References

  * [1] Christodoulos A. Floudas.  Global optimization in design and control of chemical process systems.  Journal of Process Control, 10(2):125–134, 2000.
  * [2] Wagner F. Sacco and Nélio Henderson.  Finding all solutions of nonlinear systems using a hybrid metaheuristic with fuzzy clustering means.  Applied Soft Computing, 11(8):5424–5432, 2011.
  * [3] Astrid Holstad.  Numerical solution of nonlinear equations in chemical speciation calculations.  Computational Geosciences, 3(3):229–257, Dec 1999.
  * [4] Hugo Jiménez-Islas, Gloria M. Martínez-González, José L. Navarrete-Bolaños, José E. Botello-Álvarez, and J. Manuel Oliveros-Muñoz.  Nonlinear homotopic continuation methods: A chemical engineering perspective review.  Industrial & Engineering Chemistry Research, 52(42):14729–14742, Oct 2013.
  * [5] Nélio Henderson, Wagner F. Sacco, and Gustavo Mendes Platt.  Finding more than one root of nonlinear equations via a polarization technique: An application to double retrograde vaporization.  Chemical Engineering Research and Design, 88(5):551–561, 2010.
  * [6] Hsiao-Dong Chiang, Tian-Qi Zhao, Jiao-Jiao Deng, and Kaoru Koyanagi.  Homotopy-enhanced power flow methods for general distribution networks with distributed generators.  IEEE Transactions on Power Systems, 29(1):93–100, 2014.
  * [7] Dhagash Mehta, Hung Dinh Nguyen, and Konstantin Turitsyn.  Numerical polynomial homotopy continuation method to locate all the power flow solutions.  IET Generation, Transmission & Distribution, 10(12):2972–2980, 2016.
  * [8] Hsiao-Dong Chiang and Tao Wang.  Novel homotopy theory for nonlinear networks and systems and its applications to electrical grids.  IEEE Transactions on Control of Network Systems, 5(3):1051–1060, 2018.
  * [9] Chao Yang, Jianwen Cao, and Xiao-Chuan Cai.  A fully implicit domain decomposition algorithm for shallow water equations on the cubed-sphere.  SIAM Journal on Scientific Computing, 32(1):418–438, 2010.
  * [10] Matti Schneider, Daniel Wicht, and Thomas Böhlke.  On polarization-based schemes for the fft-based computational homogenization of inelastic materials.  Computational Mechanics, 64(4):1073–1095, Oct 2019.
  * [11] Yunong Zhang.  A set of nonlinear equations and inequalities arising in robotics and its online solution via a primal neural network.  Neurocomputing, 70(1):513–524, 2006.
  * [12] Amir Salimi Lafmejani, Ahmad Kalhor, and Mehdi Masouleh.  A new development of homotopy continuation method, applied in solving nonlinear kinematic system of equations of parallel mechanisms.  2015 3rd RSI International Conference on Robotics and Mechatronics (ICROM), 10 2015.
  * [13] David A Cox, John Little, and Donal O’Shea.  Robotics and automatic geometric theorem proving.  In Ideals, Varieties, and Algorithms, pages 291–343. Springer, 2015\.
  * [14] Ping Ji and Hongtao Wu.  Kinematics analysis of an offset 3-upu translational parallel robotic manipulator.  Robotics Auton. Syst., 42:117–123, 2003.
  * [15] Strachimir Cht. Mavrodiev and Maksym Deliyergiyev.  Modification of the nuclear landscape in the inverse problem framework using the generalized bethe-weizsäcker mass formula.  arXiv e-prints, pages arXiv–1602, 2016.
  * [16] Sani Aji, Poom Kumam, Punnarai Siricharoen, Auwal Bala Abubakar, and Mahmoud Muhammad Yahaya.  A modified conjugate descent projection method for monotone nonlinear equations and image restoration.  IEEE Access, 8:158656–158665, 2020.
  * [17] Jan Verschelde, Pierre Verlinden, and Ronald Cools.  Homotopies exploiting Newton polytopes for solving sparse polynomial systems.  SIAM Journal on Numerical Analysis, 31(3):915–930, 1994.
  * [18] Crina Grosan and Ajith Abraham.  A new approach for solving nonlinear equations systems.  IEEE Transactions on Systems, Man, and Cybernetics - Part A: Systems and Humans, 38(3):698–714, 2008.
  * [19] Ahmad Golbabai, Davood Ahmadian, and Mariyan Milev.  Radial basis functions with application to finance: American put option under jump diffusion.  Mathematical and Computer Modelling, 55(3):1354–1362, 2012.
  * [20] Guang Zhang and Liang Bai.  Existence of solutions for a nonlinear algebraic system.  Discrete dynamics in nature and society, 2009, 2009.
  * [21] Karol Kowalski and Karol Jankowski.  Towards complete solutions to systems of nonlinear equations of many-electron theories.  Phys. Rev. Lett., 81:1195–1198, Aug 1998.
  * [22] Dhagash Mehta.  Numerical polynomial homotopy continuation method and string vacua.  Advances in High Energy Physics, 2011:263937, Oct 2011.
  * [23] Yang Song, Chenlin Meng, Renjie Liao, and Stefano Ermon.  Nonlinear equation solving: A faster alternative to feedforward computation.  arXiv preprint arXiv:2002.03629, 2020.
  * [24] Yunfeng Cai and Ping Li.  Solving the robust matrix completion problem via a system of nonlinear equations.  In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 4162–4172. PMLR, 26–28 Aug 2020.
  * [25] Béla Paláncz, Piroska Zaletnyik, Joseph L Awange, and Erik W Grafarend.  Dixon resultant’s solution of systems of geodetic polynomial equations.  Journal of Geodesy, 82(8):505–511, 2008.
  * [26] Béla Paláncz, Joseph L. Awange, Piroska Zaletnyik, and Robert H. Lewis.  Linear homotopy solution of nonlinear systems of equations in geodesy.  Journal of Geodesy, 84(1):79, Sep 2009.
  * [27] Christian Jansson.  An np-hardness result for nonlinear systems.  Reliable Computing, 4(4):345–350, Nov 1998.
  * [28] Yuri V. Matiyasevich.  Hilbert’s Tenth Problem.  MIT Press, Cambridge, MA, USA, 1993.
  * [29] Mustafa Q Khirallah and MA Hafiz.  Solving system of nonlinear equations using family of Jarratt methods.  International Journal of Differential Equations and Applications, 12(2), 2013.
  * [30] MQ Khirallah and MA Hafiz.  Novel three order methods for solving a system of nonlinear equations.  Bulletin of Mathematical Sciences & Applications, 1(2):01–14, 2012\.
  * [31] J Alikhani Koupaei and Seyed Mohammad Mehdi Hosseini.  A new hybrid algorithm based on chaotic maps for solving systems of nonlinear equations.  Chaos, Solitons & Fractals, 81:233–245, 2015.
  * [32] M. A. El-Shorbagy and Adel M. El-Refaey.  Hybridization of grasshopper optimization algorithm with genetic algorithm for solving system of non-linear equations.  IEEE Access, 8:220944–220961, 2020.
  * [33] Y Ramu Naidu and Akshay Kumar Ojha.  Solving multiobjective optimization problems using hybrid cooperative invasive weed optimization with multiple populations.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 48(6):821–832, 2016.
  * [34] Angel Kuri.  Solution of simultaneous non-linear equations using genetic algorithms.  WSEAS Transactions on Systems, 2, 01 2003.
  * [35] Lin Xiao, Zhijun Zhang, and Shuai Li.  Solving time-varying system of nonlinear equations by finite-time recurrent neural networks with application to motion tracking of robot manipulators.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 49(11):2210–2220, 2018.
  * [36] Huan-Tong Geng, Yi-Jie Sun, Qing-Xi Song, and Ting-Ting Wu.  Research of ranking method in evolution strategy for solving nonlinear system of equations.  In 2009 First International Conference on Information Science and Engineering, pages 348–351, 2009.
  * [37] Jun Pei, Zorica Dražić, Milan Dražić, Nenad Mladenović, and Panos M Pardalos.  Continuous variable neighborhood search (c-vns) for solving systems of nonlinear equations.  INFORMS Journal on Computing, 31(2):235–250, 2019.
  * [38] Abd allah Mousa and Islam Eldesoky.  Genls: Co-evolutionary algorithm for nonlinear system of equations.  Applied Mathematics and Computation, 197:633–642, 04 2008.
  * [39] Wenyin Gong, Yong Wang, Zhihua Cai, and Shengxiang Yang.  A weighted biobjective transformation technique for locating multiple optimal solutions of nonlinear equation systems.  IEEE Transactions on Evolutionary Computation, 21(5):697–713, 2017\.
  * [40] Wu Song, Yong Wang, Han-Xiong Li, and Zixing Cai.  Locating multiple optimal solutions of nonlinear equation systems based on multiobjective optimization.  IEEE Transactions on Evolutionary Computation, 19(3):414–431, 2015\.
  * [41] Sha Qin, Sanyou Zeng, Wei Dong, and Xi Li.  Nonlinear equation systems solved by many-objective hype.  In 2015 IEEE Congress on Evolutionary Computation (CEC), pages 2691–2696, 2015.
  * [42] Jinjin Guo, Binbin Qiu, and Yunong Zhang.  New-type dtz model for solving discrete time-dependent nonlinear equation system with robotic-arm application.  In 2020 10th International Conference on Information Science and Technology (ICIST), pages 153–162, 2020.
  * [43] Zuowen Liao, Wenyin Gong, and Ling Wang.  Memetic niching-based evolutionary algorithms for solving nonlinear equation system.  Expert Systems with Applications, 149:113261, 2020.
  * [44] Zuowen Liao, Wenyin Gong, Xuesong Yan, Ling Wang, and Chengyu Hu.  Solving nonlinear equations system with dynamic repulsion-based evolutionary algorithms.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(4):1590–1601, 2020.
  * [45] Weifeng Gao, Yuting Luo, Jingwei Xu, and Shengqi Zhu.  Evolutionary algorithm with multiobjective optimization technique for solving nonlinear equation systems.  Information Sciences, 541, 07 2020.
  * [46] Wenyin Gong, Yong Wang, Zhihua Cai, and Ling Wang.  Finding multiple roots of nonlinear equation systems via a repulsion-based adaptive differential evolution.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(4):1499–1513, 2020.
  * [47] Weifeng Gao, Genghui Li, Qingfu Zhang, Yuting Luo, and Zhenkun Wang.  Solving nonlinear equation systems by a two-phase evolutionary algorithm.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 51(9):5652–5663, 2021.
  * [48] Wei He, Wenyin Gong, Ling Wang, Xuesong Yan, and Chengyu Hu.  Fuzzy neighborhood-based differential evolution with orientation for nonlinear equation systems.  Knowledge-Based Systems, 182:104796, 2019.
  * [49] Jianye Wu, Wenyin Gong, and Ling Wang.  A clustering-based differential evolution with different crowding factors for nonlinear equations system.  Applied Soft Computing, 98:106733, 2021.
  * [50] Weifeng Gao and Yu Li.  Solving a new test set of nonlinear equation systems by evolutionary algorithm.  IEEE Transactions on Cybernetics, pages 1–10, 2021.
  * [51] Aijuan Song, Guohua Wu, and Witold Pedrycz.  Integrating variable reduction strategy with evolutionary algorithm for solving nonlinear equations systems.  CoRR, abs/2008.04223, 2020.
  * [52] Zuowen Liao, Wenyin Gong, Ling Wang, Xuesong Yan, and Chengyu Hu.  A decomposition-based differential evolution with reinitialization for nonlinear equations systems.  Knowledge-Based Systems, 191:105312, 2020.
  * [53] Shi Cheng, Junfeng Chen, Xiujuan Lei, and Yuhui Shi.  Locating multiple optima via brain storm optimization algorithms.  IEEE Access, 6:17039–17049, 2018.
  * [54] Wenyin Gong, Zuowen Liao, Xianyan Mi, Ling Wang, and Yuanyuan Guo.  Nonlinear equations solving with intelligent optimization algorithms: A survey.  Complex System Modeling and Simulation, 1(1):15–32, 2021.
  * [55] Masoud Ahookhosh, Keyvan Amini, and Somayeh Bahrami.  Two derivative-free projection approaches for systems of large-scale nonlinear monotone equations.  Numerical Algorithms, 64, 10 2013.
  * [56] Stepan Yu. Gatilov.  Properties of nonlinear systems and convergence of the Newton-Raphson method in geometric constraint solving.  volume 32, pages 57––75. NCC, 2011.
  * [57] Susumu Tanabé and Michael N. Vrahatis.  On perturbation of roots of homogeneous algebraic systems.  Mathematics of Computation, 255:1383–1402, 2006.
  * [58] Najib Ullah, Jamilu Sabi’u, and Abdullah Shah.  A derivative-free scaling memoryless broyden– fletcher–goldfarb–shanno method for solving a system of monotone nonlinear equations.  Numerical Linear Algebra with Applications, 28(5):e2374, 2021.
  * [59] José Mario Martínez.  Algorithms for Solving Nonlinear Systems of Equations, pages 81–108.  Springer Netherlands, Dordrecht, 1994.
  * [60] James M. Ortega and Werner C. Rheinboldt.  Iterative Solution of Nonlinear Equations in Several Variables.  Society for Industrial and Applied Mathematics, January 2000.
  * [61] Émile Picard.  Sur le nombre des racines communes à plusieurs équations simultanées.  Journal de Mathématiques Pures et Appliquées, 8($4^{e}$ série):5–24, 1892.
  * [62] Émile Picard.  Traité d’Analyse.  Gauthier-Villars, Paris, 1922.
  * [63] Jane Cronin.  Fixed Points and Topological Degree in Nonlinear Analysis.  American Mathematical Society, Providence, Rhode Island, 1964.
  * [64] Noel G. Lloyd.  Degree Theory.  Oxford University Press, New York, 1978.
  * [65] Donal O’Regan, Yeol Je Cho, and Yu-Qing Chen.  Topological Degree Theory and Applications.  Taylor & Francis Group, Boca Raton, FL, 2006.
  * [66] Enrique Outerelo and Jesús M. Ruiz.  Mapping Degree Theory.  American Mathematical Society, Providence, Rhode Island, 2009.
  * [67] Krzysztof A. Sikorski.  Optimal Solution of Nonlinear Equations.  Oxford University Press, New York, 2001.
  * [68] Jeroen M. Bergamin, Tassos C. Bountis, and Michael N. Vrahatis.  Homoclinic orbits of invertible maps.  Nonlinearity, 15(5):1603–1619, 2002.
  * [69] Ioannis Z. Emiris, Bernard Mourrain, and Michael N. Vrahatis.  Sign methods for counting and computing real roots of algebraic systems.  Technical report, Rapport de Recherche No.3669, INRIA (Institut National de Recherche en Informatique et en Automatique), Sophia Antipolis, France, April 1999.
  * [70] Dimitris J. Kavvadias and Michael N. Vrahatis.  Locating and computing all the simple roots and extrema of a function.  SIAM Journal on Scientific Computing, 17(5):1232–1248, 19965.
  * [71] R. Baker Kearfott.  An efficient degree-computation method for a generalized method of bisection.  Numerische Mathematik, 32:109–127, 1979.
  * [72] Bernard Mourrain, Nicos G. Pavlidis, Dimitris K. Tasoulis, and Michael N. Vrahatis.  Determining the number of real roots of polynomials through neural networks.  Computers and Mathematics with Applications, 51(3-4):527–536, 2006\.
  * [73] Bernard Mourrain, Michael N. Vrahatis, and Yakoubsohn Jean-Claude.  On the complexity of isolating real roots and computing with certainty the topological degree.  Journal of Complexity, 18(2):612–640, 2002.
  * [74] Vassilis P. Plagianakos, Nicos K. Nousis, and Michael N. Vrahatis.  Locating and computing in parallel all the simple roots of special functions using pvm.  Journal of Computational and Applied Mathematics, 133(1-2):545–554, 2001.
  * [75] Chronis Polymilis, Graziano Servizi, Charalampos Skokos, Giorgio Turchetti, and Michael N. Vrahatis.  Topological degree theory and local analysis of area preserving maps.  Chaos, 13(1):94–104, 2003.
  * [76] Frank Stenger.  Computing the topological degree of a mapping in ${R}^{n}$.  Numerische Mathematik, 25:23–38, 1975.
  * [77] Martin Stynes.  An algorithm for numerical calculation of topological degree.  Applicable Analysis, 25:23–38, 1979.
  * [78] Michael N. Vrahatis.  An efficient method for locating and computing periodic orbits of nonlinear mappings.  Journal of Computational Physics, 119:105–119, 1995.
  * [79] Michael N. Vrahatis.  Solving systems of nonlinear equations using the nonzero value of the topological degree.  ACM Transactions on Mathematical Software, 14(4):312–329, Dec. 1988\.
  * [80] Michael N. Vrahatis.  A short proof and a generalization of Miranda’s existence theorem.  Proceedings of the American Mathematical Society, 107(3):701–703, 1989.
  * [81] Michael N. Vrahatis and Kosmas I. Iordanidis.  A rapid generalized method of bisection for solving systems of nonlinear equations.  Numerische Mathematik, 49(2-3):123–138, 1986.
  * [82] Michael N. Vrahatis, Graziano Servizi, Giorgio Turchetti, and Tassos C. Bountis.  A procedure to compute the fixed points and visualize the orbits of a 2D map.  Technical report, CERN SL/93-06 (AP), CERN-SL Division, European Organization for Nuclear Research (CERN), Geneva, Switzerland, February 1993.
  * [83] Peter Franek and Marek Krčál.  Robust satisfiability of systems of equations.  J. ACM, 62(4), sep 2015.
  * [84] Erik Boman, Walter Murray, Gene Golub, and Michael Saunders.  Infeasibility and negative curvature in optimization.  03 1999.
  * [85] John W. Chinneck.  Feasibility and Infeasibility in Optimization.  Springer US, 2008.
  * [86] Abdelmonem M Ibrahim and Mohamed A Tawhid.  A hybridization of cuckoo search and particle swarm optimization for solving nonlinear systems.  Evolutionary Intelligence, 12(4):541–561, 2019.
  * [87] Kuntjoro Sidarto and Adhe Kania.  Finding all solutions of systems of nonlinear equations using spiral dynamics inspired optimization with clustering.  Journal of Advanced Computational Intelligence and Intelligent Informatics, 19:697–707, 09 2015.
  * [88] Michael J Hirsch, Panos M Pardalos, and Mauricio GC Resende.  Solving systems of nonlinear equations with continuous grasp.  Nonlinear Analysis: Real World Applications, 10(4):2000–2006, 2009\.
  * [89] Ricardo M. A. Silva, Mauricio G. C. Resende, and Panos M. Pardalos.  Finding multiple roots of a box-constrained system of nonlinear equations with a biased random-key genetic algorithm.  Journal of Global Optimization, 60(2):289–306, Oct 2014.
  * [90] Ebrahim Pourjafari and Hamed Mojallali.  Solving nonlinear equations systems with a new approach based on invasive weed optimization algorithm and clustering.  Swarm and Evolutionary Computation, 4:33–43, 2012.
  * [91] Gisela C. V. Ramadas, Edite M. G. P. Fernandes, and Ana Maria A. C. Rocha.  Multiple roots of systems of equations by repulsion merit functions.  In Computational Science and Its Applications – ICCSA 2014, pages 126–139. Springer International Publishing, 2014.
  * [92] Gisela C. V. Ramadas, Ana Maria A. C. Rocha, and Edite M. G. P. Fernandes.  Testing nelder-mead based repulsion algorithms for multiple roots of nonlinear systems via a two-level factorial design of experiments.  PLOS ONE, 10(4):1–30, 04 2015.
  * [93] Nélio Henderson, Léa Freitas, and Gustavo Platt.  Prediction of critical points: A new methodology using global optimization.  AIChE Journal, 50:1300 – 1314, 06 2004.
  * [94] Léa Freitas, Gustavo Platt, and Nélio Henderson.  Novel approach for the calculation of critical points in binary mixtures using global optimization.  Fluid Phase Equilibria, 225:29–37, 2004.
  * [95] Abolfazl Pourrajabian Shahrbabak, Ebrahimi Reza, Masoud Mirzaei, and Mehrzad Shams.  Applying genetic algorithms for solving nonlinear algebraic equations.  Applied Mathematics and Computation, 219:11483–11494, 08 2013\.
  * [96] Y. Ramu Naidu and Akshay Kumar Ojha.  Solving multiobjective optimization problems using hybrid cooperative invasive weed optimization with multiple populations.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 48(6):821–832, 2018.
  * [97] Jing-Yu Ji and Man Leung Wong.  An improved dynamic multi-objective optimization approach for nonlinear equation systems.  Information Sciences, 576:204–227, 2021.
  * [98] Jörg Stork, A. E. Eiben, and Thomas Bartz-Beielstein.  A new taxonomy of global optimization algorithms.  Natural Computing, Nov 2020.
  * [99] Rafael Martí, Panos M. Pardalos, and Mauricio G. C. Resende, editors.  Handbook of Heuristics.  Springer International Publishing, 2018.
  * [100] Dimitris Souravlias, Ilias S Kotsireas, Panos M Pardalos, and Konstantinos E Parsopoulos.  Parallel algorithm portfolios with performance forecasting.  Optimization Methods and Software, 34(6):1231–1250, 2019.
  * [101] Dimitris Souravlias, Konstantinos E Parsopoulos, Ilias S Kotsireas, and Panos M Pardalos.  Algorithm Portfolios: Advances, Applications, and Challenges.  Springer, 2021.
  * [102] Jonathan Barzilai and Jonathan M. Borwein.  Two-Point Step Size Gradient Methods.  IMA Journal of Numerical Analysis, 8(1):141–148, 01 1988.
  * [103] Roger Fletcher.  On the barzilai-borwein method.  In Optimization and control with applications, pages 235–256. Springer, 2005.
  * [104] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  CoRR, abs/1609.04747, 2016.
  * [105] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [106] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  J. Mach. Learn. Res., 12(null):2121–2159, jul 2011.
  * [107] Roger Fletcher and Colin M. Reeves.  Function minimization by conjugate gradients.  The Computer Journal, 7(2):149–154, 01 1964.
  * [108] Elijah Polak and Gerard Ribiere.  Note sur la convergence de méthodes de directions conjuguées.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 3(R1):35–43, 1969.
  * [109] Magnus Rudolph Hestenes, Eduard Stiefel, et al.  Methods of conjugate gradients for solving linear systems, volume 49.  NBS Washington, DC, 1952.
  * [110] Yu-Hong Dai and Yaxiang Yuan.  A nonlinear conjugate gradient method with a strong global convergence property.  SIAM Journal on optimization, 10(1):177–182, 1999.
  * [111] Mohammed M. Mahdi, Hasan Hadi Dwail, H. A. Wasi, Karrar Habeeb Hashim, Nabiha kahtan Dreeb, Hussein Ali Hussein, and Mushtak A. K. Shiker.  Solving systems of nonlinear monotone equations by using a new projection approach.  Journal of Physics: Conference Series, 1804(1):012107, feb 2021\.
  * [112] MM Mahdi and Mushtak AK Shiker.  A new class of three-term double projection approach for solving nonlinear monotone equations.  In Journal of Physics: Conference Series, volume 1664, page 012147\. IOP Publishing, 2020.
  * [113] Mohammed Yusuf Waziri, Aliyu Yusuf, and Auwal Bala Abubakar.  Improved conjugate gradient method for nonlinear system of equations.  Computational and Applied Mathematics, 39(4):1–17, 2020.
  * [114] MK Dauda, Mustafa Mamat, MY Waziri, Fadhila Ahmad, and Fatma Susilawati Mohamad.  Inexact CG-method via SR1 update for solving systems of nonlinear equations.  far east journal of mathematical sciences, 100(11):1787, 2016.
  * [115] MK Dauda, AS Magaji, US Shehub, MA Usman, and MY Waziri.  A simple conjugate gradient type method for solving large-scale systems of nonlinear equations.  Malaysian Journal of Computing and Applied Mathematics, 3(2):25–35, 2020.
  * [116] Howard H. Rosenbrock.  An Automatic Method for Finding the Greatest or Least Value of a Function.  The Computer Journal, 3(3):175–184, 01 1960.
  * [117] Robert Hooke and T. A. Jeeves.  “ direct search” solution of numerical and statistical problems.  J. ACM, 8(2):212–229, apr 1961.
  * [118] Richard P. Brent.  An algorithm with guaranteed convergence for finding a zero of a function.  Comput. J., 14:422–425, 1971.
  * [119] John A. Nelder and Roger Mead.  A Simplex Method for Function Minimization.  The Computer Journal, 7(4):308–313, 01 1965.
  * [120] Aijia Ouyang, Yongquan Zhou, and Qifang Luo.  Hybrid particle swarm optimization algorithm for solving systems of nonlinear equations.  In 2009 IEEE International Conference on Granular Computing, pages 460–465. IEEE, 2009.
  * [121] Richard P. Brent.  Algorithms for minimization without derivatives.  Prentice-Hall Series in Automatic Computation. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1973.
  * [122] Andrew R. Conn, Katya Scheinberg, and Luis N. Vicente.  Introduction to derivative-free optimization, volume 8 of MPS/SIAM Series on Optimization.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Programming Society (MPS), Philadelphia, PA, 2009.
  * [123] Reiner Horst, Panos M. Pardalos, and Nguyen V. Thoai.  Introduction to global optimization, volume 48 of Nonconvex Optimization and its Applications.  Kluwer Academic Publishers, Dordrecht, second edition, 2000.
  * [124] Christodoulos A. Floudas, Panos M. Pardalos, Claire S. Adjiman, William R. Esposito, Zeynep H. Gümüş, Stephen T. Harding, John L. Klepeis, Clifford A. Meyer, and Carl A. Schweiger.  Handbook of Test Problems in Local and Global Optimization.  Springer US, 1999.
  * [125] Reiner Horst and Panos M. Pardalos, editors.  Handbook of Global Optimization.  Nonconvex Optimization and Its Applications. Springer US, 1995.
  * [126] Panos M. Pardalos and H. Edwin Romeijn, editors.  Handbook of Global Optimization, Volume 2.  Springer US, 2002.
  * [127] Christodoulos A. Floudas and Panos M. Pardalos.  A collection of test problems for constrained global optimization algorithms.  Lecture Notes in Computer Science. Springer, New York, NY, September 1990\.
  * [128] Eugene L. Lawler and D. E. Wood.  Branch-and-bound methods: A survey.  Operations Research, 14(4):699–719, 1966.
  * [129] M. J. Hirsch, C. N. Meneses, P. M. Pardalos, and M. G. C. Resende.  Global optimization by continuous grasp.  Optimization Letters, 1(2):201–212, Mar 2007.
  * [130] Thomas A. Feo and Mauricio G. C. Resende.  Greedy randomized adaptive search procedures.  Journal of Global Optimization, 6(2):109–133, Mar 1995.
  * [131] R. Baker Kearfott.  Some tests of generalized bisection.  ACM Trans. Math. Softw., 13(3):197–220, September 1987.
  * [132] Keith Meintjes and Alexander P. Morgan.  Chemical equilibrium systems as numerical test problems.  ACM Trans. Math. Softw., 16(2):143–151, June 1990.
  * [133] Nenad Mladenović and Pierre Hansen.  Variable neighborhood search.  Computers & Operations Research, 24(11):1097–1100, 1997.
  * [134] Gisela C.V. Ramadas and Edite M.G.P. Fernandes.  Self-adaptive combination of global tabu search and local search for nonlinear equations.  International Journal of Computer Mathematics, 89(13-14):1847–1864, 2012.
  * [135] Scott Kirkpatrick, Daniel Gelatt, and Mario P. Vecchi.  Optimization by simulated annealing.  Science, 220(4598):671–680, 1983.
  * [136] V. Granville, M. Krivánek, and J. P. Rasson.  Simulated annealing: A proof of convergence.  IEEE Trans. Pattern Anal. Mach. Intell., 16(6):652–656, jun 1994\.
  * [137] A. Corana, M. Marchesi, C. Martini, and S. Ridella.  Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm.  ACM Trans. Math. Softw., 13(3):262–280, sep 1987.
  * [138] Hime A. e Oliveira and Antonio Petraglia.  Solving nonlinear systems of functional equations with fuzzy adaptive simulated annealing.  Applied Soft Computing, 13(11):4349–4357, 2013.
  * [139] Ofer M. Shir.  Niching in Evolutionary Algorithms, pages 1035–1069.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.
  * [140] Mike Preuss.  Multimodal optimization by means of evolutionary algorithms.  Springer, 2015.
  * [141] James C. Bean.  Genetic algorithms and random keys for sequencing and optimization.  ORSA Journal on Computing, 6(2):154–160, 1994.
  * [142] M.J. Hirsch, P.M. Pardalos, and M.G.C. Resende.  Speeding up continuous GRASP.  European Journal of Operational Research, 205(3):507–521, September 2010.
  * [143] Mohamed A Tawhid and Abdelmonem M Ibrahim.  A hybridization of grey wolf optimizer and differential evolution for solving nonlinear systems.  Evolving Systems, 11(1):65–87, 2020.
  * [144] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Particle swarm optimization and intelligence: Advances and applications.  Information Science Publishing (IGI Global), Hershey, PA, USA, 2010.
  * [145] Yugui Li, Yanxu Wei, and Yantao Chu.  Research on solving systems of nonlinear equations based on improved pso.  Mathematical Problems in Engineering, 2015, 2015.
  * [146] Oguz Emrah Turgut, Mert Sinan Turgut, and Mustafa Turhan Coban.  Chaotic quantum behaved particle swarm optimization algorithm for solving nonlinear system of equations.  Computers & Mathematics with Applications, 68(4):508–530, 2014\.
  * [147] Jun Sun, Bin Feng, and Wenbo Xu.  Particle swarm optimization with particles having quantum behavior.  In Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753), volume 1, pages 325–331 Vol.1, 2004.
  * [148] Yuanbin Mo, Hetong Liu, and Qin Wang.  Conjugate direction particle swarm optimization solving systems of nonlinear equations.  Computers & Mathematics with Applications, 57(11-12):1877–1882, 2009.
  * [149] Ivan Amaya, Jorge Cruz, and Rodrigo Correa.  Real roots of nonlinear systems of equations through a metaheuristic algorithm.  Dyna, 78(170):15–23, 2011.
  * [150] Majid Jaberipour, Esmaile Khorram, and Behrooz Karimi.  Particle swarm algorithm for solving systems of nonlinear equations.  Computers & Mathematics with Applications, 62(2):566–576, 2011\.
  * [151] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Parameter selection and adaptation in unified particle swarm optimization.  Mathematical and Computer Modelling, 46(1-2):198–213, 2007.
  * [152] Konstantinos E. Parsopoulos and Michael N. Vrahatis.  Recent approaches to global optimization problems through particle swarm optimization.  Natural Computing, 1(2-3):235–306, 2002.
  * [153] Yiannis G. Petalas, Konstantinos E. Parsopoulos, and Michael N. Vrahatis.  Memetic particle swarm optimization.  Annals of Operations Research, 156(1):99–127, 2007.
  * [154] Kenichi Tamura and Keiichiro Yasuda.  Spiral dynamics inspired optimization.  Journal of Advanced Computational Intelligence and Intelligent Informatics, 15(8):1116–1122, 2011.
  * [155] Kenichi Tamura and Keiichiro Yasuda.  Primary study of spiral dynamics inspired optimization.  IEEJ Transactions on Electrical and Electronic Engineering, 6(S1):S98–S100, 2011.
  * [156] Rüdiger U. Seydel.  Tools for Computational Finance.  Springer London, 2012.
  * [157] Mahdi Abdollahi, Asgarali Bouyer, and Davoud Abdollahi.  Improved cuckoo optimization algorithm for solving systems of nonlinear equations.  The Journal of Supercomputing, 72(3):1246–1269, 2016.
  * [158] Xinming Zhang, Qian Wan, and Youhua Fan.  Applying modified cuckoo search algorithm for solving systems of nonlinear equations.  Neural Computing and Applications, 31(2):553–576, 2019.
  * [159] Gai-Ge Wang, Suash Deb, and Zhihua Cui.  Monarch butterfly optimization.  Neural Computing and Applications, 31(7):1995–2014, Jul 2019.
  * [160] Abdelmonem M Ibrahim and Mohamed A Tawhid.  A hybridization of differential evolution and monarch butterfly optimization for solving systems of nonlinear equations.  Journal of Computational Design and Engineering, 6(3):354–367, 2019\.
  * [161] Yongquan Zhou, Jiakun Liu, and Guangwei Zhao.  Leader glowworm swarm optimization algorithm for solving nonlinear equations systems.  Electrical Review, 88(1):101–106, 2012.
  * [162] Ruimin Jia and Dengxu He.  Hybrid artificial bee colony algorithm for solving nonlinear system of equations.  In 2012 Eighth international conference on computational intelligence and security, pages 56–60. IEEE, 2012.
  * [163] Mahdi Abdollahi, Ayaz Isazadeh, and Davoud Abdollahi.  Imperialist competitive algorithm for solving systems of nonlinear equations.  Computers & Mathematics with Applications, 65(12):1894–1908, 2013\.
  * [164] MKA Ariyaratne, TGI Fernando, and Sunethra Weerakoon.  Solving systems of nonlinear equations using a modified firefly algorithm (modfa).  Swarm and Evolutionary Computation, 48:72–92, 2019.
  * [165] Gisela CV Ramadas and Edite Manuela da GP Fernandes.  Solving systems of nonlinear equations by harmony search.  International Conference on Mathematical Methods in Science and Engineering …, 2013.
  * [166] Naser Moosavian and Babak Kasaee Roodsari.  Soccer league competition algorithm, a new method for solving systems of nonlinear equations.  International Journal of Intelligence Science, 04(01):7–16, 2014\.
  * [167] Ya-Zhong Luo, Guo-Jin Tang, and Li-Ni Zhou.  Hybrid approach for solving systems of nonlinear equations using chaos optimization and quasi-newton method.  Applied Soft Computing, 8(2):1068–1073, 2008.
  * [168] Dixiong Yang, Gang Li, and Gengdong Cheng.  On the efficiency of chaos optimization algorithms for global optimization.  Chaos Solitons and Fractals, 34(4):1366–1375, November 2007.
  * [169] J. Kiefer.  Sequential minimax search for a maximum.  Proceedings of the American Mathematical Society, 4(3):502–506, 1953.
  * [170] Yu-Nong Zhang, Zeng-Hai Chen, and Ke Chen.  Convergence properties analysis of gradient neural network for solving online linear equations.  Acta Automatica Sinica, 35(8):1136–1139, 2009.
  * [171] Guimei Li and Zhezhao Zeng.  A neural-network algorithm for solving nonlinear equation systems.  volume 1, pages 20–23, 12 2008.
  * [172] Karl Mathia.  Solutions of linear equations and a class of nonlinear equations using recurrent neural networks.  Technical report, January 2000.
  * [173] Guillaume Lample and François Charton.  Deep learning for symbolic mathematics.  CoRR, abs/1912.01412, 2019.
  * [174] Athanasios Margaris and Miltiadis Adamopoulos.  Solving nonlinear algebraic systems using artificial neural networks.  CEUR Workshop Proceedings, 284, 01 2007.
  * [175] Tasneem A. Awaad, Abdelrahman M. Elbehery, Amany Abdelhamid, Salma K. Elsokkary, Youssef M. Ali, Khaled Salah, Mohamed Abdel Salam, and M. Watheq El-Kharashi.  An ultrafast neural network-based hardware acceleration for nonlinear systems’ simulators.  Computers & Electrical Engineering, 79:106452, 2019.
  * [176] Konstantinos Goulianas, Athanasios Margaris, Ioannis Refanidis, Konstantinos Diamantaras, and Theofilos Papadimitriou.  A back propagation-type neural network architecture for solving the complete nxn nonlinear algebraic system of equations.  2016\.
  * [177] Long Jin, Shuai Li, Bolin Liao, and Zhijun Zhang.  Zeroing neural networks: A survey.  Neurocomputing, 267:597–604, 2017.