Nonlinear Programming Solvers for Unconstrained and Constrained Optimization Problems: a Benchmark Analysis

In this paper we propose a set of guidelines to select a solver for the solution of nonlinear programming problems. With this in mind, we present a comparison of the convergence performances of commonly used solvers for both unconstrained and constrained nonlinear programming problems. The comparison involves accuracy, convergence rate, and convergence speed. Because of its popularity among research teams in academia and industry, MATLAB is used as common implementation platform for the solvers. Our study includes solvers which are either freely available, or require a license, or are fully described in literature. In addition, we diﬀerentiate solvers if they allow the selection of diﬀerent optimal search methods. As result, we examine the performances of 23 algorithms to solve 60 benchmark problems. To enrich our analysis, we will describe how, and to what extent, convergence speed and accuracy can be improved by changing the inner settings of each solver.


Introduction
The current technological era prioritizes, more than ever, high performance and efficiency of complex processes controlled by a set of variables.Examples of these processes are (Lasdon and Warren 1980;Grossmann 1996;Charalambous 1979;Grossmann and Kravanja 1997;Wu and William 1992;Wansuo and Haiying 2010;Rustagi 1994;Ziemba and Vickson 1975): engineering designs, chemical plant reactions, manufacturing processes, grid power management, power generation/conversion process, path planning for autonomous vehicles, climate simulations, etc. Quite often, the search for the best performance, or the highest efficiency, can be transcribed into the form of a Nonlinear Programming (NLP) problem.Namely, the need to minimize (or maximize) a scalar cost function subjected to a set of constraints.In some instances these functions are linear but, in general, one or both of them are characterized by nonlinearities.For simple, one-time use problems, one might successfully use any of the solver available, like fmincon in MATLAB (MathWorks 2020a; MATLAB 2020).Nevertheless, if the NLP derives from some specific applications, like real-time process optimization, then the solver choice begs a more accurate selection.
The first research efforts toward the characterization of optimization solvers started in the 60s'.In (Box 1966) the authors compare eight solvers on twenty benchmark unconstrained NLP problems containing up to 20 variables.Notably, they illustrate techniques to transform particular constrained NLP into equivalent unconstrained problems.The authors of (Levy and Guerra 1976) analyze the convergence properties of two gradient-based solvers applied to 16 test problems.In the last few decades, with the development of new methodologies and optimization applications, more studies aimed to illustrate difference in performance among NLP solvers.Schittkowski et.al.(Schittkowski, Zillober, and Zotemantel 1994) performed the comparison of eleven different mathematical programming codes applied to structural optimization through finite element analysis.George et.al. summarize a qualitative comparison of few optimization methodologies reported by several other sources (George and Raimond 2013).In the research document prepared by Sandia National Laboratory (Gearhart et al. 2013), a study was conducted on four open source Linear Programming (LP) solvers applied to 201 benchmark problems.In (Kronqvist et al. 2018) Kronqvist et.al. carried out a performance comparison of mixed integer NLP solvers limited to convex benchmark problems.Pratiksha Saxena presents a comparison between linear and nonlinear programming techniques on the diet formulation for animals (Saxena 2012).Another work by Hannes Pucher uses the programming language R to analyze multiple nonlinear optimization methods applied to real life problems (Pucher and Stix 2008).State-of-the-art optimization methods were used to compare on their application on L1-regularized classifiers (Yuan et al. 2010).On a similar note, multiple global optimization solvers were compared on a work done by Arnold Neumaier (Neumaier et al. 2005).Authors from (Obayash and Tsukahara 1997;McIlhagga, Husbands, and Ives 1996;Haupt 1995;Hamdy, Nguyen, and Hensen 2016) conducted a performance comparison of optimization techniques for specific applications which includes aerodynamic shape design, integrated manufacturing planning and scheduling, solving electromagnetic problems, and building energy design problems, respectively.Similarly, Frank et. al. conducted a comparison between three optimization methods for solving aerodynamic design problems (Frank and Shubin 1992).In (Karaboga and Basturk 2008), Karaboga et. al. compares the performance of artificial bee colony algorithm with the differential evolution, evolutionary and particle swarm optimization algorithms using multi-dimensional numerical problems.
In this paper we want to provide an explicit comparison of a set of NLP solvers.We include in our comparison popular solvers readily available in MATLAB, a few gradient descent methods that have been extensively used in literature, and a particle swarm optimization.Because of its widespread use among research groups, both in academia and private sector, we have decided to use MATLAB as common implementation platform.For this reason we will focus on all the solvers that are either written on or can be implemented in MATLAB.The NLP problems used in this comparison have been selected amongst the standard benchmark problems (Hedar 2020;Schittkowski 2009;Floudas and Pardalos 1999) with up to thirty variables and a up to nine scalar constraints.The paper is organized as follows.Section 2 describes the statement of unconstrained and contrained NLP problems.In Section 3, we enumerate the NLP solvers included in our analysis and their main features.Subsequently, an overview of the different convergence metrics, and the solvers implementations is carried out in Section 4. The results of the comparison with the benchmark equations is discussed in Section 5. Finally, the main contributions of the paper are outlined in Section 6.

Nonlinear programming problem statements
In general, a constrained NLP problem aims to minimize a nonlinear real scalar objective function, with respect a set of variables, while satisfying a set of nonlinear constraints.If the problem entails the minimization of a function without the presence of constraints the problem is defined as unconstrained (Nocedal and Wright 2006).In the following section, the general form of a nonlinear unconstrained and constrained optimization problems in the minimization format form are thoroughly stated.

Statement
Let x ∈ R n be a real vector with n ≥ 1 components and let f : R n → R be a smooth function.Then, the unconstrained optimization problem is defined as (1)

Optimality conditions
For a one-dimensional function f (x) defined and differentiable over an interval (a, b), the necessary condition for a point x * ∈ (a, b) to be a local maximum or minimum is that f ′ (x * ) = 0.This is also known as Fermat's theorem.The multidimensional extension of this condition states that the gradient must be zero at local optimum point, namely ∇f (x * ) = 0. (2) Eq. 2 is referred as a first order optimality condition, as it expresses in terms of the first order derivatives.

Statement
The constrained optimization problem is formulated as subject to With c(x) a smooth real-valued function on a subset of R n .Notably, c i (x) and c j (x) represent the sets of equality constraints and inequality constraints, respectively.The feasible set is identified as the set of points x that satisfy just the constraints (Eqs.4, 5).It must be pointed out that some of the solvers considered in this study are only able to consider equality constraints.In these instances, we will introduce a set of slack variables s i , and convert Eq. 5 into the following set of equality constraints Such necessary expedient will obviously induce more computational burden on the particular solvers affected by this constraint-type limitation.

Optimality conditions
The measure of first-order optimality for constrained problems derives from the Karush-Kuhn-Tucker (KKT) conditions (Boyd and Vandenberghe 2004).These necessary conditions are defined as follow.Let the objective function f and the constraint functions g i and h j be continuously differentiable functions at x * ∈ R n .If x * is a local optimum and the optimization problem satisfies some regularity conditions (Nocedal and Wright 2006), then there exist the two constants µ i (i = 1, . . ., w) and λ j (j = 1, . . ., ℓ), called KKT multipliers, such that the following four groups of conditions hold: • Primal feasibility: h j (x * ) = 0, for j = 1, . . ., ℓ.

Selection of NLP solvers and algorithms
The selection of the NLP solvers considered in this work is based on the following aspects.First of all, we are only considering algorithms that can be implemented in MATLAB.Secondly, we have included solvers that are either free source or, for commercial software, have a free trial version.The remaining part of this section briefly describes the 23 solvers included in our analysis and the most direct source to each algorithm.

APSO
The Accelerated Particle Swarm Optimization (APSO) is an algorithm developed by Yang at Cambridge University in 2007 and it based on swarm-intelligent search of the optimum (Yang 2014).APSO is an evolution of the standard particle swarm optimization, and developed to accelerate the convergence of the standard version of the algorithm.The standard PSO is characterized by two elements, the swarm, that is the population, and the members of the population, called particles.The search is based on a randomly initialized population that moves in randomly chosen directions.In particular, each particle moves through the searching space, remembers the best earlier positions, velocity, and accelerations of itself and its neighbors.This information is shared among the particles while they dynamically adjust their own position, velocity and acceleration derived from the best position of all particles.The next step starts when all particles have been shifted.Finally, all particles aim to find the global best among all the current best solutions till the objective function no longer improves or after a certain number of iterations (Yang 2014).The standard PSO uses both the current global best and the individual best, whereas the simplified version APSO is able to accelerate the convergence of the algorithm buy using the global best only.Due to the nature of the algorithm, only constrained nonlinear programming problems can be solved.The MATLAB version of the APSO algorithm is provided in (Yang 2014).

BARON
The Branch and Reduced Optimization Navigator (BARON) is a commercial global optimization software that solves both NLP and mixed-integer nonlinear programs (MINLP).BARON uses deterministic global optimization algorithms of the branch and bound search type which, by applying general assumptions, solve the global optimization problem.It comes with embedded linear programming (LP) and NLP solvers, such as CLP/CBC, IPOPT, FilterSD and FilterSQP.By default, BARON selects the NLP solver and may switch between different NLP solvers during the search according to problem characteristics and solver performance.To refer to the default option, the name BARON (auto) is chosen.Unlike many other NLP algorithms, BARON doesn't explicitly require the user to provide an initial guess of the solution but leaves this as an option.If a user doesn't provide the initial guess, then the software shrewdly initializes the variables.In this paper, we use the demo version of the software in conjunction with the MATLAB interface which can be retrieved in (Firm 2021).Must be noted that the free demo version is characterized by some limitations, namely, it can only handle problems with up to ten variables, ten constraints, and it doesn't support trigonometric functions.Details and documentations about BARON software are provided in (Tawarmalani and Sahinidis 2004;Sahinidis n.d.).

CLP/CBC
The Computational Optimization Infrastructure for Operations Research (COIN-OR) Branch and Cut (CBC) is an open-source mixed integer nonlinear programming solver based on the COIN-OR LP solver (CLP) and the COIN-OR Cut generator library (Cgl).The code has been written primarily by John J. Forrest (COIN-OR 2016).

IPOPT
COIN-OR Interior Point Optimizer (IPOPT) is an open-source solver for large-scale NLP and it has been mainly developed by Andreas Wächter (Wächter and Biegler 2006).IPOPT implements an interior point line search filter method for nonlinear programming models.The problem function are not required to be convex but should be twice continuously differentiable.Mathematical details of the algorithm and documentation can be found in (COIN-OR 2021).

FilterSD
FilterSD is a package of Fortran 77 subroutines for solving nonlinear programming problems and linearly constrained problems in continuous optimization.The NLP solver filterSD aims to find a solution of the NLP problem, where the objective function and the constraint function are continuously differentiable at points that satisfy the bounds on x.The code has been developed to avoid the use of second derivatives, and to prevent storing an approximate reduced Hessian matrix by using a new limited memory spectral gradient approach based on Ritz values.The basic approach is that of Robinson's method, globalised by using a filter and trust region (FilterSD 2020).

FilterSQP
FilterSQP is a Sequential Quadratic Programming solver suitable for solving large, sparse or dense linear, quadratic and nonlinear programming problems.The method implements a trust region algorithm with a filter to promote global convergence.The filter accepts a trial point whenever the objective or the constraint violation is improved compared to all previous iterations.The size of the trust region is reduced if the step is rejected, and increased if it is accepted (Fletcher and Leyffer 1999).

FMINCON
FMINCON is a MATLAB optimization toolbox used to solve constrained nonlinear programming problems.FMINCON provides the user the option to select amongst five different algorithms to solve nonlinear problems: Active-set, Interior-point, Sequential Quadratic Programming, Sequential Quadratic Programming legacy, and Trust region reflective.Four, out of the five, algorithms are implemented in our analysis as one of them, the Trust Region Reflective algorithm, does not accept the type of constraint considered in our benchmark cases.

Active-set
The Active-set, unlike the Interior point (mentioned next), doesn't use a barrier term to ensure that the inequality constraints are met, but it solves the optimal equation by understanding the true active-set.A general active-set algorithm for convex quadratic programming can be found in (Nocedal and Wright 2006).

Interior point
This method, also known as barrier method, is one type of nonlinear problem-solving algorithms that achieves the determination of optimum values by iteratively approaching the optimal solution from the interior of the feasible set (Nocedal and Wright 2006).Since interior point algorithm depends on a feasible set, the following requirements must be met for the method to be used: • the set of feasible interior point should not be empty; • all the iterations should occur in the interior of this feasible set.

Sequential Quadratic Programming
The basic idea behind the Sequential Quadratic Programming (SQP) is to find a minimizer for a subproblem, which is generated as an approximate model of the optimization problem at the current iteration point.This will be then used to define a new iteration point, which is in turn used to define another minimizer, and the process is iterated.SQP is similar to the active set, but some of the differences are listed as follows: • strict feasibility with respect to bounds; • robustness to non-double results; • refactored linear algebra routines; • reformulated feasibility routines.
A general line-search algorithm framework for SQP can be found in (Nocedal and Wright 2006).

Sequential Quadratic Programming legacy
Sequential Quadratic Programming legacy (SQP-legacy) is similar to SQP, with the difference of using a larger memory and, therefore, it is slower to determining the problem solution (Nocedal and Wright 2006).

Trust region reflective
The Trust region reflective algorithm solves a NLP by defining a region that is assumed to represent the objective function as accurately as possible.From the selected trust region, a step is taken and is used as a minimizer.If that specific step doesn't generate a solution, a different region with a reduced size will be selected.Then a new step is executed and considered as new minimizer for the region, and the process is iterated (Nocedal and Wright 2006).Trust region reflective algorithm used by FMINCON requires only bounds or linear equality constraints.Due to this setback, this algorithm isn't included in the analysis.

FMINUNC
FMINUNC is another MATLAB optimization toolbox used to solve unconstrained nonlinear programming problems (MathWorks 2020b).In this case, FMINUNC gives the user the option of choosing between two different algorithms to solve nonlinear minimization problems: Quasi-Newton, and Trust region.

Quasi-Newton
The Quasi-Newton methods build up curvature information at each iteration to formulate a quadratic model problem, with the optimal solution occurring when the sationarity conditions are satisfied.Newton-type methods, as opposed to quasi-Newton methods, calculate the Hessian matrix directly, and proceed in a direction of descent to locate the minimum after a number of iterations, numerically involving a large amount of computation.On the contrary, quasi-Newton methods adopt the observed behavior of the objective function and its gradient to build up curvature information to make an approximation to the Hessian matrix using an appropriate updating technique (MathWorks 2020c).In particular, the quasi-Newton algorithm uses the formula of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) to implement a cubic line search procedure, and for updating the approximation of the Hessian matrix (MathWorks 2020b).

Trust region
The trust region algorithm is a subspace trust-region method, based on the interiorreflective Newton method.Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients (PCG) (MathWorks 2020b).In a minimization context, the Hessian matrix can be assumed symmetric, but it is guaranteed to be positive definite only in the neighborhood of a strong minimizer.Algorithm PCG exits when it encounters a direction of negative or zero curvature.The PCG output direction is either a direction of negative curvature or an approximate solution to the Newton system, anyways helping to define the two-dimensional subspace used in the trust-region approach (MathWorks 2020c).

GCMMA
GCMMA, the Globally Convergent Method of Moving Asymptotes, is a modified version of the MMA that evaluates the global optimum value.However, unlike MMA, GMMA consists of the so called inner and outer iterations.The GCMMA follows the same steps as the MMA except for small changes.In GMMA, an approximate subproblem is created for the first outer iteration by replacing the function with convex functions.These subproblems are then solved in order to find the next iteration points otherwise the inner iteration kicks off.For the first inner iteration, a new subproblem will be generated.Then the next iteration points will be calculated by solving these convex subproblems.The algorithm then moves to the next iteration (Svanberg 2002).The GCMMA algorithm is fully described in (Svanberg 2007), and the MATLAB code is freely available at (Svanberg 2020).

KNITRO
ARTELYS KNITRO is a commercially available nonlinear optimization software package developed by Zienna Optimization since 2001 (Knitro 2021b).KNITRO, short for Nonlinear Interior point Trust Region Optimization, is a software package for finding local solutions of both continuous optimization problems, with or without constraints, and discrete optimization problems with integer or binary variables.The KNITRO package provides efficient and robust solution of small or large problems, for both continuous and discrete problems, derivative-free options.It supports the most popular operating systems and several modeling language and programmatic interfaces (Knitro 2021a).Multiple versions of the software are available to download at (Knitro 2021b).In this work, the software free trial license is used, in conjunction with the MATLAB interface.Several algorithms are included in the software, such as Interior point, Activeset, and Sequential Quadratic Programming.The description of these algorithms can be found in Section 3.3.

MIDACO
The Mixed Integer Distributed Ant Colony Optimization (MIDACO) is a global optimization solver that combines an extended evolutionary probabilistic technique, called the Ant Colony Optimization algorithm, with the Oracle Penalty method for constrained handling (MIDACO-Solver, user manual 2021).Ant Colony optimization is modelled using the behavior of ants to find the quickest path between their colony and the source food.Like the majority of evolutionary optimization algorithms, MIDACO considers the objective and constraint functions as black-box functions.MIDACO was created in collaboration of European Space Agency and EADS Astrium to solve constrained mixed-integer non-linear (MINLP) space applications (Schlueter et al. 2013).We use the trial version of MIDACO, in conjunction with the MATLAB interface.The trial version has a limitation, namely, it doesn't support more than four variables per problem.The solver can be downloaded from (MIDACO-Solver, user manual 2021).

MMA
The Method of Moving Asymptotes (MMA) solves nonlinear problem function by generating an approximate subproblem.These convex functions used as subproblems are chosen using gradient information at the current iteration points, and also at parameters that are updated at each iteration stage, called the moving asymptotes.The subproblem is solved at the current iteration point, and the solution is used as the next iteration point.Similarly, a new subproblem is generated at this new iteration point, which again is solved to create the next iteration point (Svanberg 1987).The MMA algorithm is fully described in (Svanberg 2007), and the MATLAB code is freely available at (Svanberg 2020).

MQA
The Modified Quasilinearization Algorithm (MQA) is the modified version of the Standard Quasilinearization Algorithm (SQA) (Eloe and Jonnalagadda 2019; Yeo 1974) described below.These quasilinearization algorithms base their solution search on the linear approximation of the NLP, namely, on the Hessian matrix and gradient of the objective and constraint functions.Ultimately, the goal is the progressive reduction of the performance index.For unconstrained NLP problems, the performance index is defined as Q = f T x f x , where f x is the gradient of the objective function.On the other hand, for constrained NLP problems the performance index is defined as R = P + Q, which comprises both the feasibility index P = h T h, and optimality index Q = F T x F x , with F = f +λ T h, where f is the objective function, h is the constraint function, and λ is the vector of Lagrange multipliers associate to the constraint function.Convergence to the desired solution is achieved when the performance index Q ≤ ε 1 or R ≤ ε 2 , with ε 1 and ε 2 small preselected positive constants, for the unconstrained and constrained case respectively (Miele and Iyer 1971;Miele, Mangiavacchi, and Aggarwal 1974).Unlike SQA, characterized by a unitary step size, MQA reduces progressively the step size 0 < α < 1 to enforce an improvement in optimality.In turn, the main advantage of the MQA, over the SQA, is its descent property: if the stepsize α is sufficiently small, the reduction in the performance index is indeed guaranteed.It must be pointed out that the MQA for NLP problems can only treat equality constraints.Therefore, in our implementation, all the inequality constraints are converted into equality constraints by introducing the slack variables.We have implemented the algorithm on MATLAB in order to solve both unconstrained and constrained NLP problems, based on (Miele and Iyer 1971;Miele, Mangiavacchi, and Aggarwal 1974).

PENLAB
PENLAB is a free open source software package implemented in MATLAB for nonlinear optimization, linear and nonlinear semidefinite optimization and any combination of these.It derives from PENNON, the original implementation of the algorithm which is not open source (Fiala, Kočvara, and Stingl 2013).Originally, PENNON was an implementation of the PBM method developed by Ben-Tal and Zibulevsky for problems of structural optimization, that has grown into a stand alone program for solving general problems (Kocvara and Stingl 2003).It is based on a generalized Augmented Lagrangian method pioneered by R. Polyak (Polyak 1992).PENLAB can be freely downloaded from (Kocvara 2017).

SGRA
The Sequential Gradient-Restoration Algorithm (SGRA) is a first order nonlinear programming solver developed by Angelo Miele and his research group in 1969 (COKER 1985;Miele, Huang, and Heideman 1969).It is based on a cyclical scheme whereby, first, the constraints are satisfied to a prescribed accuracy (restoration phase); then, using a first-order gradient method, a step is taken toward the optimal direction to improve the performance index (gradient phase).The performance index is defined as R = P + Q, which includes both the feasibility index P = h T h, and optimality index Q = F T x F x , with F = f + λ T h, where f is the objective function, h is the constraint function, and λ is the vector of Lagrange multipliers associate to the constraint function.Convergence is achieved when the constraint error, and the optimality condition error are P ≤ ε 1 , Q ≤ ε 2 , respectively, with ε 1 , ε 2 small preselected positive constants.It must be pointed out that the SGRA for NLP problems can only treat equality constraints.Therefore, in our implementation, all the inequality constraints are converted into equality constraints by introducing the slack variables.We have programmed the algorithm on MATLAB in order to solve both unconstrained and constrained NLP problems, based on (Miele, Huang, and Heideman 1969).The SGRA version used to solve unconstrained NLP problems differs from the original formulation by the omission of the restoration phase in the iterative process.

SNOPT
The Sparse Nonlinear OPTimizer (SNOPT) is a commercial software package for solving large-scale optimization problems, linear and nonlinear programs.It minimizes a linear or nonlinear function subject to bounds on the variables and sparse linear or non-linear constraints.SNOPT implements a sequential quadratic programming method for solving constrained optimization problems with functions and gradients that are expensive to evaluate, and with smooth nonlinear functions in the objective and constraints (Gill et al. 2001).SNOPT is implemented in Fortran 77 and distributed as source code.In this paper, we use the free trial version of the software in conjunction with the MATLAB interface, that can be retrieved at (Gill et al. 2018).

SOLNP
SOLNP is originally implemented in MATLAB to solve general nonlinear programming problems, characterized by nonlinear smooth functions in the objective and constraints (Ye 1989).Inequality constraints are converted into equality constraints by means of slack variables.The major iteration of SOLNP solves a linearly constrained optimization problem with an augmented Lagrangian objective function.Within the major iteration, as first step it is checked if the solution is feasible for the linear equality constraints of the objective function; if it is not, an interior linear programming procedure is called to find an interior feasible (or near-feasible) solution.Successively, a sequential quadratic programming (QP) solves the linearly constrained problem.If the QP solution is both feasible and optimal, the algorithm stops, otherwise it solves another QP problem as minor iteration.Both major and minor processes repeat until the optimal solution is found or the user-specified maximum number of iterations is reached (Ye 1989).The SOLNP module in MATLAB can be freely downloaded from (Ye 2020).

SQA
The Standard Quasilinearization Algorithm (SQA) is the standard version of the QA, and it uses QA techniques for solving nonlinear problems by generating a sequence of linear problems solutions (Eloe and Jonnalagadda 2019; Yeo 1974).SQA differs from MQA for the value associated to the scaling factor α. As mentioned before, the SQA can only treat equality constraints.Therefore, in our implementation, all the inequality constraints are converted into equality constraints by introducing the slack variables.

Convergence metrics and solvers implementation
In this section we provide the description of the convergence metrics considered in our analysis and narrate the key implementation steps for each solver.

Convergence metrics
The main goal of this paper is to characterize the convergence performance, in terms of speed and accuracy, of the different solvers under analysis.We have selected a number of benchmark NLPs and compared the numerical solutions returned by each solver with the true analytical solution.Moreover, considering that the choice of the initial guesses critically affects the convergence process, we want to also assess the capability to converge to the true optimum, rather than converging into to local minima or not converging at all.With this in mind, we define as converging robustness the quality of a solver to achieve the solution when the search process is initiated from a broad set of initial guesses randomly chosen within the search domain.Finally, to have an accurate assessment of the convergence speed, we will require the solver to repeat the same search several times and average out the total CPU time.As result, given N benchmark test functions, M solvers/algorithms, K randomly generated initial guesses, and Z repeated identical search runs, a total of N × M × K × Z runs have been executed.
The following performance metrics are in order: • Mean error [%]: with f (x) the benchmark test function evaluated at the numerical solution x provided by the solver, f (x * ) the benchmark test function evaluated at the optimal solution f (x * ), E k the error associated to the run from the k-th randomly generated initial guess, Ēn the mean error associated to the n-th benchmark test function, and Ēm the mean error delivered by the m-th solver.The biunivocal choice of the denominator of E k is based on the fact that some benchmark test functions at the optimal solution have zero value; in this case, a value of 0.001 is chosen instead as reference value.
where σ n is the variance correspondent to the n-th benchmark test function, and σm the mean variance delivered by the m-th solver.
• Mean convergence rate [%]: with K conv the number of runs (from a pool of K distinct initial guesses) which successfully reach convergence for the n-th function, γ n the convergence rate for the n-th function, and γm the mean convergence rate delivered by the m-th solver.The convergence rate is computed considering succesfull a run that satisfies the converging threshold conditions E k ≤ E max = 5%, and CP U k ≤ CP U max = 10 s, with CP U k is the CPU time required to the run starting from the k-th initial guess.

• Mean CPU time [s]:
where CP U z is the mean CPU time per z-th repetition, CP U n is the mean CPU time related to the n-th benchmark test function, and CP U m is the mean CPU time delivered by the m-th solver.

Solvers implementation
In this paper we analyze the convergence performances of the different solvers in terms of robustness, accuracy, and convergence speed.Considering that the user might decide to tune the convergence parameters to favor one of these metrics, we have decided to perform the comparison for three separate implementation scenarios: plug and play (P&P), high accuracy (HA), and quick solution (QS).The plug and play settings, as the name suggests, are the "out-of-the-box" settings of each solver.The high accuracy settings are based on more stringent tolerances and/or on a higher number of maximum iterations with respect to the plug and play settings.This tuning aims to achieve a more precise solution.Finally, the quick solution settings are characterized by more relaxed convergence tolerances, and a lower number of maximum iterations with respect to the plug and play settings.In this scenario the algorithms should reach a less accurate solution but in a shorter time.In general, the objective function, its gradient, the initial conditions, the constraint function (for constrained problems only), and the solver options are elements which are inputted to each solver.The objective function gradient is not necessary for APSO, BARON, MIDACO, and SOLNP, but it is optional for FMINCON/FMINUNC and KNITRO.For GCMMA/MMA, SGRA, and SNOPT, the gradient of both the objective and constraint functions is necessary.MQA/SQA and PENLAB, in addition to these inputs, require the Hessian of the objective function.
In the following subsection, details on each solver, and on the three different solver settings per each solver are described.It must be noted that, in most cases, the settings' names here reported are the same of the solver's options names used in the code implementation.In this way, the reader can have a better understanding of which solver's parameter has been tuned.

APSO
The three settings considered in the analysis are reported in Table 1, where no.particles is the number of particles, no.iterations is the total number of iterations, and γ is a control parameter that multiplies α, one of the two learning parameters or acceleration constants, α and β, the random amplitude of roaming particles and the speed of convergence, respectively.APSO does also require the number of problem variables, no.vars, to be defined but this parameter is, obviously, invariant for the three settings.

BARON
The three settings considered in the analysis are reported in Table 2, with EpsA the absolute termination tolerance, EpsR the relative termination tolerance, and AbsConF easT ol the absolute constraint feasibility tolerance.Due to the limitations of the trial version of the solver, trigonometric functions and problems with more than ten variables are not supported by the solver; for this reason, the following test functions are excluded in the analysis: A.2, A.3, A.4, A.5, A.7, A.11, A.13, A.14, A.16, A.17, A.18, A.22, A.24, A.26 for unconstrained problems, and B.1, B.2, B.5, B.8, B.20 for constrained problems.

FMINCON/FMINUNC
The three settings considered in the analysis are reported in Table 3, with StepT olerance the lower bound on the size of a step, ConstraintT olerance the upper bound on the magnitude of any constraint functions, F unctionT olerance the lower bound on the change in the value of the objective function during a step, and OptimalityT olerance the tolerance for the first-order optimality measure.

GCMMA/MMA
The three settings considered in the analysis are reported in Table 4, where epsimin is a prescribed small positive tolerance that terminates the algorithm, whereas maxoutit is the maximum number of iterations for MMA, and the maximum number of outer iterations for GCMMA.

KNITRO
The three settings considered in the analysis are reported in Table 5, where M axIter is the maximum number of iterations before termination, T olX is a tolerance that terminates the optimization process if the relative change of the solution point estimate is less than that value, T olF un specifies the final relative stopping tolerance for the KKT (optimality) error, and T olCon specifies the final relative stopping tolerance for the feasibility error.

MIDACO
The three settings considered in the analysis are reported in Table 6, where maxeval is the maximum number of function evaluation.It is a distinctive feature of MIDACO that allows the solver to stop exactly after that number of function evaluation.Due to the limitations of the trial version of the solver, test functions with more than four variables are not supported by the solver; for this reason, the following test functions are excluded in the analysis: A.7, A.10, A.11, A.13, A.14, A.15, A.16, A.18, A.22, A.24 for unconstrained problems, and B.1, B.2, B.3, B.4, B.7, B.9, B.10, B.12, B.18, B.19, B.20, B.21, B.22, B.26, B.29, B.30 for constrained problems.

MQA
The three settings considered in the analysis are reported in Table 7, with ε 1 and ε 2 the prescribed small positive tolerances that allow the solver to stop, when the inequality Q ≤ ε 1 or R ≤ ε 2 is met.As mentioned in Section 3.9, MQA for NLP problems can only treat equality constraints, namely all the inequality constraints are converted into equality constraints by introducing the slack variables.In this study, for all the three settings considered in the analysis, a value of 1 is chosen as initial guess for all the slack variables.

Settings P&P HA
QS 1e-4 1e-5 1e-3 4.2.8.PENLAB The three settings considered in the analysis are reported in Table 8, where max_inner_iter is the maximum number of inner iterations, max_outer_iter is the maximum number of outer iterations, mpenalty_min is the lower bound for penalty parameters, inner_stop_limit is the termination tolerance for the inner iterations, outer_stop_limit is the termination tolerance for the outer iterations, kkt_stop_limit is the termination tolerance KKT optimality conditions, and unc_dir_stop_limit is the stopping tolerance for the unconstrained minimization.

Settings P&P
HA QS The three settings considered in the analysis are reported in Table 10, where major_iterations_limit is the limit on the number of major iterations in the SQP method, minor_iterations_limit is the limit on minor iterations in the QP subproblems, major_f easibility_tolerance is the tolerance for feasibility of the nonlinear constraints, major_optimality_tolerance is the tolerance for the dual variables, and minor_f easibility_tolerance is the tolerance for the variables and their bounds.

SOLNP
The three settings considered in the analysis are reported in Table 11, with ρ the penalty parameter in the augmented Lagrangian objective function, maj the maximum number of major iterations, min the maximum number of minor iterations, δ the perturbation parameter for numerical gradient calculation, and ǫ the relative tolerance on optimality and feasibility.During the HA scenario implementation, we learned that different convergence settings are required for unconstrained and constrained problems.This peculiarity might be induced by the stringent tolerances adopted in this scenario.
Table 11.SOLNP settings.Tuning values for the HA scenario are divided for unconstrained (left-side) and constrained (right-side) problems.

Settings P&P
HA QS ρ 1 1 1 maj 10 500|10 10 min 10 500|10 10 δ 1e-5 1e-10|1e-6 1e-3 ǫ 1e-4 1e-12|1e-7 1e-3 4.2.12.SQA The three settings considered in the analysis are reported in Table 12, with ε 1 and ε 2 the prescribed small positive tolerances that allow the solver to stop, when the inequality Q ≤ ε 1 or R ≤ ε 2 is met.As mentioned earlier, SQA can only treat equality constraints.To overcome this limitation, the inequality constraints are converted into equality constraints by introducing slack variables.In this study, for all the three settings considered in the analysis, a value of 1 is chosen for all the slack variables.

Benchmark test functions and results
We present a collection of unconstrained and constrained optimization test problems that are used to validate the performance of the various optimization algorithms presented above for the different implementation scenarios.The comparison results are also discussed in depth in this section.For performance comparison purposes, an equivalent environment and control parameters have been created to run each NLP solver.All outputs tabulated in this paper are calculated using MATLAB software running on a desktop computer with the following specs: Intel(R) Core(TM) i7-6700 CPU 3.40GHz processor, 16.0 GB of RAM, running a 64-bit Windows 10 operating system.To assess the true computational time required by each algorithm to reach convergence, implementation steps that are expected to have an impact on the computer's performance are deactivated during the run for the solution.The internet connection and other unrelated applications are turned off throughout the analysis, ensuring that unnecessary background activities are not accessing computational resources throughout the solvers' performance.A collection of unconstrained and constrained benchmark problems are used to test the solvers based on, but not limited to (Hedar 2020;Schittkowski 2009;Floudas and Pardalos 1999).Specifically, the benchmark problems include combinations of logarithmic, trigonometric, and exponential terms, non-convex and convex functions, a minimum of two to a maximum of thirty variables, and a maximum of nine constraint functions for the constrained optimization problems.For sake of completeness, all the benchmark test functions are listed in Appendix A and B. For each of the test function, dimension, domain and search space, objective function, constraints, and minimum solution are listed.As mentioned in Section 4.2, the comparison between each solver is carried out by considering three different settings: plug and play, high accuracy, and quick solution.In this way we want to assess the robustness, accuracy, and convergence speed of every solver.For each benchmark problem, all solvers use the same set of randomly generated initial guesses.

Results for unconstrained optimization problems
A collection of 30 unconstrained optimization test problems is used to validate the performance of the optimization algorithms.The benchmark test functions for unconstrained global optimization are listed in Appendix A. For the purpose of this analysis, given N = 30 benchmark test functions, M = 16 solvers and algorithms, K = 50 randomly generated initial guesses, and Z = 3 iterations, a set of N × M × K × Z runs are executed.Tables 13, 14, and 15 report the results for the plug and play (P&P), high accuracy (HA), and quick solution (QS) settings, respectively.From the analysis of the results for the P&P settings, Table 13, we observe that BARON (auto) and BARON (ipopt) are able to reach the minimum mean error and variance, the highest convergence rate, but they are not the fastest ones to reach the solution.Moreover, BARON (sd), BARON (sqp), SNOPT, and PENLAB are able to obtain good results in terms of mean error and variance.Overall, PENLAB is also able to reach a convergence rate similar to BARON (auto) and BARON (ipopt), with the advantage of being 10 times faster than them.The worst results in terms of accuracy and convergence rate are obtained by SOLNP and SGRA.For the HA settings, Table 14, we can observe similar trends.In general, as expected, all the solvers manage to achieve a more accurate solution as they reduce the average error, increase their convergence rate, and increase the average convergence time.MIDACO is now able to reach the highest convergence rate together with all the versions of BARON.Overall PENLAB is the solver which delivers a good trade-off in performance.With respect the P&P settings, SOLNP significantly improves its convergence rate, whereas SGRA just slightly increase its performances.It is interesting to observe that, KNITRO (interior-point) and KNITRO (sqp), aside improving their convergence rate, increase their mean error and variance increase.Despite our effort, we are not sure how to explain this unexpected behaviour.Regarding the QS settings, Table 15, generally all the solvers reduce their convergence time and also decrease their convergence rate except for BARON (auto), BARON (ipopt), and BARON (sqp) which remain unaltered.SQA, FMINUNC, SOLNP, and FMINUNC are amongst the fastest to reach the solution but their convergence rate is quite low.In addition, conversely to all the other solvers that experience a smaller CPU time, BARON is not always able to achieve a faster CPU time with respect to the P&P settings.The same happens to the SGRA, probably due to its intrinsic iterative nature.

Conclusions
In this paper we provide an explicit comparison of a set of NLP solvers.The comparison includes popular solvers which are readily available in MATLAB, a few gradient descent methods that have been extensively used in literature, and a particle swarm optimization.Because of its widespread use among research groups, both in academia and private sector, we have used MATLAB as common implementation platform.Constrained and unconstrained NLP problems have been selected amongst the standard benchmark problems with up to thirty variables and a up to nine scalar constraints.
Results for the unconstrained problems show that BARON is the algorithm that deliver the best convergence rate and accuracy but it is the slowest.PENLAB is the algorithm that has the best trade off between accuracy, convergence rate, and speed.For the constrained NLP problems, again, BARON is the solver which delivers excellent accuracy and convergence rate but is amongst the slower.FMINCON, KNITRO, SNOPT, and MIDACO are the one that are able to deliver a fair compromise of accuracy, convergence rate, and speed.

Table 14 .
All unconstrained problems, high accuracy (HA) settings.Solvers ranked w.r.t.mean error.

Table 15 .
All unconstrained problems, quick solution (QS) settings.Solvers ranked w.r.t.mean CPU time.A collection of 30 unconstrained optimization test problems is used to validate the performance of the optimization algorithms.The benchmark test functions for constrained global optimization are listed in Appendix B. For the purpose of the analysis, given N = 30 benchmark test functions, M = 21 solvers and algorithms, K = 50 randomly generated initial guesses, and Z = 3 iterations, a set of N × M × K × Z runs are executed.Tables16, 17, and 18 report the results for the P&P, HA, and QS settings, respectively.From the analysis of the results for the P&P settings, Table16, we observe that all the versions of BARON are able to reach almost the highest accuracy and the best convergence rate but they are not the fastest to reach the solution.MIDACO is able to achieve the second best convergence rate, with an average CPU time that is more than 50% faster than BARON.PENLAB obtains the best mean error and variance but this performance is tempered by a low convergence rate, together with the SGRA, MQA, and SQA which are also quite slow to reach solution.FMINCON (interior-point), KNITRO (interior-point), and SNOPT reach a convergence rate lower than BARON and MIDACO, but they are significantly faster.Regarding the HA settings, Table17, similar consideration can be made for BARON, but in this case the CPU time is considerably increasing.MIDACO shows an improvement in the convergence rate, reaching values very similar to BARON.PENLAB still obtains the best mean error and variance, but it has one of the lowest convergence rate, together with the SGRA.In general, most of the solvers increase their convergence rate, and decrease their mean error, except for GCMMA and PENLAB.Regarding the QS settings, Table 18, generally all the solvers decrease their convergence rate except for BARON and PENLAB.Same considerations about BARON and PENLAB can be done as in the two previous scenarios.MIDACO reports a significant decrease in the convergence rate.The different versions of BARON have similar CPU time with respect to the P&P settings.FMINCON (interior-point), KNITRO (interior-point), and SNOPT reach a convergence rate lower than BARON, but they are significantly faster.The worst results in terms of convergence rate and CPU time are obtained by MQA and SQA.

Table 17 .
All constrained problems, high accuracy (HA) settings.Solvers ranked w.r.t.mean error.

Table 18 .
All constrained problems, quick solution (QS) settings.Solvers ranked w.r.t.mean CPU time.