Main Content

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`

SolverTo 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: '...'
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`

SolverTo 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. Optimization terminated: mesh size 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: 'Optimization terminated: mesh size 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`

SolverTo 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. Optimization terminated: maximum number of generations exceeded.

`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: 'Optimization terminated: maximum number of generations exceeded.'
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`

SolverTo 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`

SolverTo 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. Optimization terminated: change in best function value less than options.FunctionTolerance.

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

fvalsim = 1.8311e-05

eflagsim = SolverConvergedSuccessfully

`outputsim = `*struct with fields:*
iterations: 1967
funccount: 1986
message: 'Optimization terminated: change in best function value less than options.FunctionTolerance.'
rngstate: [1x1 struct]
problemtype: 'unconstrained'
temperature: [2x1 double]
totaltime: 1.9293
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: -0.0033
y: 4.7219e-04

fvalsur = 2.2456e-05

eflagsur = SolverLimitExceeded

`outputsur = `*struct with fields:*
elapsedtime: 6.8243
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.

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.0017923 1.8311e-05 1986 surrogateopt -0.0033311 0.00047219 2.2456e-05 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 take 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.

`solve`

| `patternsearch`

| `ga`

| `particleswarm`

| `simulannealbnd`

| `surrogateopt`