Compare Several Global Solvers, Problem-Based

This example shows how to minimize Rastrigin’s function with several solvers. Each solver has its own characteristics. The characteristics lead to different solutions and run times. The results, summarized in Compare Solvers 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):

ras = @(x, y) 20 + x.^2 + y.^2 - 10*(cos(2*pi*x) + cos(2*pi*y));

Plot the function scaled by 10 in each direction.

rf3 = @(x, y) ras(x/10, y/10);
fsurf(rf3,[-30 30],"ShowContours","on")
title("rastriginsfcn([x/10,y/10])")
xlabel("x")
ylabel("y")

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.

fminunc Solver

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

x = optimvar("x");
y = optimvar("y");
prob = optimproblem("Objective",rf3(x,y));
x0.x = 20;
x0.y = 30;
[solf,fvalf,eflagf,outputf] = solve(prob,x0)

Solving problem using fminunc.

Local minimum found.

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

solf = struct with fields:
    x: 19.8991
    y: 29.8486

fvalf = 
12.9344

eflagf = 
    OptimalSolution

outputf = struct with fields:
             iterations: 3
              funcCount: 5
               stepsize: 1.7773e-06
           lssteplength: 1
          firstorderopt: 2.0461e-09
              algorithm: 'quasi-newton'
                message: 'Local minimum found....'
    objectivederivative: "reverse-AD"
                 solver: 'fminunc'

fminunc solves the problem in very few function evaluations (only five, as shown in the outputf structure), and reaches a local minimum near the start point. The exit flag indicates that the solution is a local minimum.

patternsearch Solver

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

x0.x = 20;
x0.y = 30;
[solp,fvalp,eflagp,outputp] = solve(prob,x0,"Solver","patternsearch")

Solving problem using patternsearch.
patternsearch stopped because the mesh size was less than options.MeshTolerance.

solp = struct with fields:
    x: 19.8991
    y: -9.9496

fvalp = 
4.9748

eflagp = 
    SolverConvergedSuccessfully

outputp = struct with fields:
         function: @(x)fun(x,extraParams)
      problemtype: 'unconstrained'
       pollmethod: 'gpspositivebasis2n'
    maxconstraint: []
     searchmethod: []
       iterations: 48
        funccount: 174
         meshsize: 9.5367e-07
         rngstate: [1x1 struct]
          message: 'patternsearch stopped because the mesh size was less than options.MeshTolerance.'
           solver: 'patternsearch'

Like fminunc, patternsearch reaches a local optimum, as shown in the exit flag exitflagp. The solution is better than the fminunc solution, because it has a lower objective function value. However, patternsearch takes many more function evaluations, as shown in the output structure.

ga Solver

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

rng default % For reproducibility
x0.x = 10*randn(20) + 20;
x0.y = 10*randn(20) + 30; % Random start population near [20,30];
[solg,fvalg,eflagg,outputg] = solve(prob,"Solver","ga")

Solving problem using ga.
ga stopped because it exceeded options.MaxGenerations.

solg = struct with fields:
    x: 0.0064
    y: 7.7057e-04

fvalg = 
8.1608e-05

eflagg = 
    SolverLimitExceeded

outputg = struct with fields:
      problemtype: 'unconstrained'
         rngstate: [1x1 struct]
      generations: 200
        funccount: 9453
          message: 'ga stopped because it exceeded options.MaxGenerations.'
    maxconstraint: []
       hybridflag: []
           solver: 'ga'

ga takes many more function evaluations than the previous solvers, and arrives at a solution near the global minimum. The solver is stochastic, and can reach a suboptimal solution.

particleswarm Solver

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

rng default % For reproducibility
[solpso,fvalpso,eflagpso,outputpso] = solve(prob,"Solver","particleswarm")

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

solpso = struct with fields:
    x: 7.1467e-07
    y: 1.4113e-06

fvalpso = 
4.9631e-12

eflagpso = 
    SolverConvergedSuccessfully

outputpso = struct with fields:
      rngstate: [1x1 struct]
    iterations: 120
     funccount: 2420
       message: 'Optimization ended: relative change in the objective value ...'
    hybridflag: []
        solver: 'particleswarm'

The solver takes fewer function evaluations than ga, and arrives at an even more accurate solution. Again, the solver is stochastic and can fail to reach a global solution.

simulannealbnd Solver

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

rng default % For reproducibility
x0.x = 20;
x0.y = 30;
[solsim,fvalsim,eflagsim,outputsim] = solve(prob,x0,"Solver","simulannealbnd")

Solving problem using simulannealbnd.
simulannealbnd stopped because the change in best function value is less than options.FunctionTolerance.

solsim = struct with fields:
    x: 0.0025
    y: 0.0018

fvalsim = 
1.8386e-05

eflagsim = 
    SolverConvergedSuccessfully

outputsim = struct with fields:
     iterations: 1967
      funccount: 1986
        message: 'simulannealbnd stopped because the change in best function value is less than options.FunctionTolerance.'
       rngstate: [1x1 struct]
    problemtype: 'unconstrained'
    temperature: [2x1 double]
      totaltime: 1.1845
         solver: 'simulannealbnd'

The solver takes about the same number of function evaluations as particleswarm, and reaches a good solution. This solver, too, is stochastic.

surrogateopt Solver

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
x = optimvar("x","LowerBound",-70,"UpperBound",130);
y = optimvar("y","LowerBound",-70,"UpperBound",130);
prob = optimproblem("Objective",rf3(x,y));
options = optimoptions("surrogateopt","PlotFcn",[]);
[solsur,fvalsur,eflagsur,outputsur] = solve(prob,...
    "Solver","surrogateopt",...
    "Options",options)

Solving problem using surrogateopt.
surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

solsur = struct with fields:
    x: 9.9494
    y: -9.9502

fvalsur = 
1.9899

eflagsur = 
    SolverLimitExceeded

outputsur = struct with fields:
        elapsedtime: 3.7364
          funccount: 200
    constrviolation: 0
               ineq: [1x1 struct]
           rngstate: [1x1 struct]
            message: 'surrogateopt stopped because it exceeded the function evaluation limit set by ...'
             solver: 'surrogateopt'

The solver takes relatively few function evaluations to reach a solution near the global optimum. However, each function evaluation takes much more time than those of the other solvers.

Compare Solvers and Solutions

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

sols = [solf.x solf.y;
    solp.x solp.y;
    solg.x solg.y;
    solpso.x solpso.y;
    solsim.x solsim.y;
    solsur.x solsur.y];
fvals = [fvalf;
    fvalp;
    fvalg;
    fvalpso;
    fvalsim;
    fvalsur];
fevals = [outputf.funcCount;
    outputp.funccount;
    outputg.funccount;
    outputpso.funccount;
    outputsim.funccount;
    outputsur.funccount];
stats = table(sols,fvals,fevals);
stats.Properties.RowNames = ["fminunc" "patternsearch" "ga" "particleswarm" "simulannealbnd" "surrogateopt"];
stats.Properties.VariableNames = ["Solution" "Objective" "# Fevals"];
disp(stats)

                              Solution            Objective     # Fevals
                      ________________________    __________    ________

    fminunc               19.899        29.849        12.934         5  
    patternsearch         19.899       -9.9496        4.9748       174  
    ga                 0.0063672    0.00077057    8.1608e-05      9453  
    particleswarm     7.1467e-07    1.4113e-06    4.9631e-12      2420  
    simulannealbnd      0.002453     0.0018028    1.8386e-05      1986  
    surrogateopt          9.9494       -9.9502        1.9899       200  

These results are typical:

• fminunc quickly reaches the local solution within its starting basin, but does not explore outside this basin at all. Because the objective function has analytic derivatives, fminunc uses automatic differentiation and takes very few function evaluations to reach an accurate local minimum.

• patternsearch takes more function evaluations than fminunc, and searches through several basins, arriving at a better solution than fminunc.

• ga takes many more function evaluations than patternsearch. By chance it arrives at a better solution. In this case, ga finds a point near the global optimum. ga is stochastic, so its results change with every run. ga requires 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 finds a point with lower objective function value than patternsearch or ga. Because particleswarm is stochastic, its results change with every run. particleswarm requires extra steps to have an initial population near [20,30].

• simulannealbnd takes about the same number of function evaluations as particleswarm. In this case, simulannealbnd finds a good solution, but not as good as particleswarm. The solver is stochastic and can arrive at a suboptimal solution.

• surrogateopt stops when it reaches a function evaluation limit, which by default is 200 for a two-variable problem. surrogateopt requires finite bounds. surrogateopt attempts to find a global solution, and in this case succeeds. 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.