Comparison of Six Solvers

Function to Optimize

This example shows how to minimize Rastrigin’s function with six solvers. Each solver has its own characteristics. The characteristics lead to different solutions and run times. The results, examined in Compare Syntax and Solutions, can help you choose an appropriate solver for your own problems.

Rastrigin’s function has many local minima, with a global minimum at (0,0):

$$\text{Ras}(x)=20+x_1^2+x_2^2-10\left(\cos 2\pi x_1+\cos 2\pi x_2\right).$$

Usually you don't know the location of the global minimum of your objective function. To show how the solvers look for a global solution, this example starts all the solvers around the point [20,30], which is far from the global minimum.

The rastriginsfcn.m file implements Rastrigin’s function. This file comes with Global Optimization Toolbox software. This example employs a scaled version of Rastrigin’s function with larger basins of attraction. For information, see Basins of Attraction.

rf2 = @(x)rastriginsfcn(x/10);

 Code for generating the figure

rf3 = @(x,y)reshape(rastriginsfcn([x(:)/10,y(:)/10]),size(x));
fsurf(rf3,[-30 30],'ShowContours','on')
title('rastriginsfcn([x/10,y/10])')
xlabel('x')
ylabel('y')

This example minimizes rf2 using the default settings of fminunc (an Optimization Toolbox™ solver), patternsearch, and GlobalSearch. The example also uses ga and particleswarm with nondefault options to start with an initial population around the point [20,30]. Because surrogateopt requires finite bounds, the example uses surrogateopt with lower bounds of -70 and upper bounds of 130 in each variable.

Six Solution Methods

• fminunc

• patternsearch

• ga

• particleswarm

• surrogateopt

• GlobalSearch

fminunc

To solve the optimization problem using the fminunc Optimization Toolbox solver, enter:

rf2 = @(x)rastriginsfcn(x/10); % objective
x0 = [20,30]; % start point away from the minimum
[xf,ff,flf,of] = fminunc(rf2,x0)

fminunc returns

Local minimum found.

Optimization completed because the size of the gradient is 
less than the default value of the function tolerance.

xf =
   19.8991   29.8486
ff =
   12.9344
flf =
     1
of =

  struct with fields:

       iterations: 3
        funcCount: 15
         stepsize: 1.7776e-06
     lssteplength: 1
    firstorderopt: 5.9907e-09
        algorithm: 'quasi-newton'
          message: 'Local minimum found.…'

• xf is the minimizing point.

• ff is the value of the objective, rf2, at xf.

• flf is the exit flag. An exit flag of 1 indicates xf is a local minimum.

• of is the output structure, which describes the fminunc calculations leading to the solution.

patternsearch

To solve the optimization problem using the patternsearch Global Optimization Toolbox solver, enter:

rf2 = @(x)rastriginsfcn(x/10); % objective
x0 = [20,30]; % start point away from the minimum
[xp,fp,flp,op] = patternsearch(rf2,x0)

patternsearch returns

Optimization terminated: mesh size less than options.MeshTolerance.
xp =
   19.8991   -9.9496
fp =
    4.9748
flp =
     1
op =

  struct with fields:

        function: @(x)rastriginsfcn(x/10)
     problemtype: 'unconstrained'
      pollmethod: 'gpspositivebasis2n'
   maxconstraint: []
    searchmethod: []
      iterations: 48
       funccount: 174
        meshsize: 9.5367e-07
        rngstate: [1x1 struct]
         message: 'Optimization terminated: mesh size less than options.MeshTolerance.'

• xp is the minimizing point.

• fp is the value of the objective, rf2, at xp.

• flp is the exit flag. An exit flag of 1 indicates xp is a local minimum.

• op is the output structure, which describes the patternsearch calculations leading to the solution.

ga

To solve the optimization problem using the ga Global Optimization Toolbox solver, enter:

rng default % for reproducibility
rf2 = @(x)rastriginsfcn(x/10); % objective
x0 = [20,30]; % start point away from the minimum
initpop = 10*randn(20,2) + repmat([10 30],20,1);
opts = optimoptions('ga','InitialPopulationMatrix',initpop);
[xga,fga,flga,oga] = ga(rf2,2,[],[],[],[],[],[],[],opts)

initpop is a 20-by-2 matrix. Each row of initpop has mean [10,30], and each element is normally distributed with standard deviation 10. The rows of initpop form an initial population matrix for the ga solver.

opts is the options that set initpop as the initial population.

The final line calls ga, using the options.

ga uses random numbers, and produces a random result. In this case ga returns:

Optimization terminated: average change in the fitness value
less than options.FunctionTolerance.
xga =
    0.0236   -0.0180
fga =
    0.0017
flga =
     1
oga =

  struct with fields:

    problemtype: 'unconstrained'
       rngstate: [1x1 struct]
    generations: 107
      funccount: 5400
        message: 'Optimization terminated: average change in the fitness value less...'
  maxconstraint: []

• xga is the minimizing point.

• fga is the value of the objective, rf2, at xga.

• flga is the exit flag. An exit flag of 1 indicates xga is a local minimum.

• oga is the output structure, which describes the ga calculations leading to the solution.

particleswarm

Like ga, particleswarm is a population-based algorithm. So for a fair comparison of solvers, initialize the particle swarm to the same population as ga.

rng default % for reproducibility
rf2 = @(x)rastriginsfcn(x/10); % objective
opts = optimoptions('particleswarm','InitialSwarmMatrix',initpop);
[xpso,fpso,flgpso,opso] = particleswarm(rf2,2,[],[],opts)

Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.

xpso =

   1.0e-06 *

   -0.8839    0.3073


fpso =

   1.7373e-12


flgpso =

     1


opso =

  struct with fields:

      rngstate: [1x1 struct]
    iterations: 114
     funccount: 2300
       message: 'Optimization ended: relative change in the objective value …'

• xpso is the minimizing point.

• fpso is the value of the objective, rf2, at xpso.

• flgpso is the exit flag. An exit flag of 1 indicates xpso is a local minimum.

• opso is the output structure, which describes the particleswarm calculations leading to the solution.

surrogateopt

surrogateopt does not require a start point, but does require finite bounds. Set bounds of –70 to 130 in each component. To have the same sort of output as the other solvers, disable the default plot function.

rng default % for reproducibility
lb = [-70,-70];
ub = [130,130];
rf2 = @(x)rastriginsfcn(x/10); % objective
opts = optimoptions('surrogateopt','PlotFcn',[]);
[xsur,fsur,flgsur,osur] = surrogateopt(rf2,lb,ub,opts)

Surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

xsur =

   -9.9516   -9.9486


fsur =

    1.9899


flgsur =

     0


osur = 

  struct with fields:

       rngstate: [1×1 struct]
      funccount: 200
    elapsedtime: 3.0454
        message: 'Surrogateopt stopped because it exceeded the function evaluation limit set by ↵'options.MaxFunctionEvaluations'.'

• xsur is the minimizing point.

• fsur is the value of the objective, rf2, at xsur.

• flgsur is the exit flag. An exit flag of 0 indicates that surrogateopt halted because it ran out of function evaluations or time.

• osur is the output structure, which describes the surrogateopt calculations leading to the solution.

GlobalSearch

To solve the optimization problem using the GlobalSearch solver, enter:

rf2 = @(x)rastriginsfcn(x/10); % objective
x0 = [20,30]; % start point away from the minimum
problem = createOptimProblem('fmincon','objective',rf2,...
    'x0',x0);
gs = GlobalSearch;
[xg,fg,flg,og] = run(gs,problem)

problem is an optimization problem structure. problem specifies the fmincon solver, the rf2 objective function, and x0=[20,30]. For more information on using createOptimProblem, see Create Problem Structure.

Note

You must specify fmincon as the solver for GlobalSearch, even for unconstrained problems.

gs is a default GlobalSearch object. The object contains options for solving the problem. Calling run(gs,problem) runs problem from multiple start points. The start points are random, so the following result is also random.

In this case, the run returns:

GlobalSearch stopped because it analyzed all the trial points.

All 4 local solver runs converged with a positive local solver exit flag.

xg =

   1.0e-07 *

   -0.1405   -0.1405


fg =

     0


flg =

     1


og = 

  struct with fields:

                funcCount: 2128
         localSolverTotal: 4
       localSolverSuccess: 4
    localSolverIncomplete: 0
    localSolverNoSolution: 0
                  message: 'GlobalSearch stopped because it analyzed all the trial points.↵↵All 4 local solver runs converged with a positive local solver exit flag.'

• xg is the minimizing point.

• fg is the value of the objective, rf2, at xg.

• flg is the exit flag. An exit flag of 1 indicates all fmincon runs converged properly.

• og is the output structure, which describes the GlobalSearch calculations leading to the solution.

Compare Syntax and Solutions

One solution is better than another if its objective function value is smaller than the other. The following table summarizes the results, accurate to one decimal.

Results	fminunc	patternsearch	ga	particleswarm	surrogateopt	GlobalSearch
solution	[19.9 29.9]	[19.9 -9.9]	[0 0]	[0 0]	[-9.9 -9.9]	[0 0]
objective	12.9	5	0	0	2	0
# Fevals	15	174	5400	2300	200	2178

These results are typical:

• fminunc quickly reaches the local solution within its starting basin, but does not explore outside this basin at all. fminunc has a simple calling syntax.

• patternsearch takes more function evaluations than fminunc, and searches through several basins, arriving at a better solution than fminunc. The patternsearch calling syntax is the same as that of fminunc.

• ga takes many more function evaluations than patternsearch. By chance it arrived at a better solution. In this case, ga found a point near the global optimum. ga is stochastic, so its results change with every run. ga has a simple calling syntax, but there are extra steps to have an initial population near [20,30].

• particleswarm takes fewer function evaluations than ga, but more than patternsearch. In this case, particleswarm found the global optimum. particleswarm is stochastic, so its results change with every run. particleswarm has a simple calling syntax, but there are extra steps to have an initial population near [20,30].

• surrogateopt stops when it reaches a function evaluation limit, which by default is 200 for a two-variable problem. surrogateopt has a simple calling syntax, but requires finite bounds. Although surrogateopt attempts to find a global solution, in this case the returned solution is not the global solution. Each function evaluation in surrogateopt takes a longer time than in most other solvers, because surrogateopt performs many auxiliary computations as part of its algorithm.

• GlobalSearch run takes the same order of magnitude of function evaluations as ga and particleswarm, searches many basins, and arrives at a good solution. In this case, GlobalSearch found the global optimum. Setting up GlobalSearch is more involved than setting up the other solvers. As the example shows, before calling GlobalSearch, you must create both a GlobalSearch object (gs in the example), and a problem structure (problem). Then, you call the run method with gs and problem. For more details on how to run GlobalSearch, see Workflow for GlobalSearch and MultiStart.

Related Topics

• Optimization Problem Setup
• Solver Behavior with a Nonsmooth Problem