This example illustrates how `GlobalSearch`

finds
a global minimum efficiently, and how `MultiStart`

finds
many more local minima.

The objective function for this example has many local minima and a unique global minimum. In polar coordinates, the function is

*f*(*r*,*t*)
= *g*(*r*)*h*(*t*),

where

$$\begin{array}{c}g(r)=\left(\mathrm{sin}(r)-\frac{\mathrm{sin}(2r)}{2}+\frac{\mathrm{sin}(3r)}{3}-\frac{\mathrm{sin}(4r)}{4}+4\right)\frac{{r}^{2}}{r+1}\\ h(t)=2+\mathrm{cos}(t)+\frac{\mathrm{cos}\left(2t-\frac{1}{2}\right)}{2}\text{.}\end{array}$$

The global minimum is at *r* = 0, with objective
function 0. The function *g*(*r*)
grows approximately linearly in *r*, with a repeating
sawtooth shape. The function *h*(*t*)
has two local minima, one of which is global.

Write a function file to compute the objective:

function f = sawtoothxy(x,y) [t r] = cart2pol(x,y); % change to polar coordinates h = cos(2*t - 1/2)/2 + cos(t) + 2; g = (sin(r) - sin(2*r)/2 + sin(3*r)/3 - sin(4*r)/4 + 4) ... .*r.^2./(r+1); f = g.*h; end

Create the problem structure. Use the

`'sqp'`

algorithm for`fmincon`

:problem = createOptimProblem('fmincon',... 'objective',@(x)sawtoothxy(x(1),x(2)),... 'x0',[100,-50],'options',... optimoptions(@fmincon,'Algorithm','sqp','Display','off'));

The start point is

`[100,-50]`

instead of`[0,0]`

, so`GlobalSearch`

does not start at the global solution.Validate the problem structure by running

`fmincon`

:[x,fval] = fmincon(problem) x = 45.7332 -107.6469 fval = 555.5422

Create the

`GlobalSearch`

object, and set iterative display:gs = GlobalSearch('Display','iter');

Run the solver:

rng(14,'twister') % for reproducibility [x,fval] = run(gs,problem) Num Pts Best Current Threshold Local Local Analyzed F-count f(x) Penalty Penalty f(x) exitflag Procedure 0 200 555.5 555.5 0 Initial Point 200 1463 1.547e-15 1.547e-15 1 Stage 1 Local 300 1564 1.547e-15 5.858e+04 1.074 Stage 2 Search 400 1664 1.547e-15 1.84e+05 4.16 Stage 2 Search 500 1764 1.547e-15 2.683e+04 11.84 Stage 2 Search 600 1864 1.547e-15 1.122e+04 30.95 Stage 2 Search 700 1964 1.547e-15 1.353e+04 65.25 Stage 2 Search 800 2064 1.547e-15 6.249e+04 163.8 Stage 2 Search 900 2164 1.547e-15 4.119e+04 409.2 Stage 2 Search 950 2356 1.547e-15 477 589.7 387 2 Stage 2 Local 952 2420 1.547e-15 368.4 477 250.7 2 Stage 2 Local 1000 2468 1.547e-15 4.031e+04 530.9 Stage 2 Search GlobalSearch stopped because it analyzed all the trial points. 3 out of 4 local solver runs converged with a positive local solver exit flag. x = 1.0e-07 * 0.0414 0.1298 fval = 1.5467e-15

You can get different results, since

`GlobalSearch`

is stochastic.

The solver found three local minima, and it found the global
minimum near `[0,0]`

.

Write a function file to compute the objective:

function f = sawtoothxy(x,y) [t r] = cart2pol(x,y); % change to polar coordinates h = cos(2*t - 1/2)/2 + cos(t) + 2; g = (sin(r) - sin(2*r)/2 + sin(3*r)/3 - sin(4*r)/4 + 4) ... .*r.^2./(r+1); f = g.*h; end

Create the problem structure. Use the

`fminunc`

solver with the`Algorithm`

option set to`'quasi-newton'`

. The reasons for these choices are:The problem is unconstrained. Therefore,

`fminunc`

is the appropriate solver; see Optimization Decision Table (Optimization Toolbox).The default

`fminunc`

algorithm requires a gradient; see Choosing the Algorithm (Optimization Toolbox). Therefore, set`Algorithm`

to`'quasi-newton'`

.

problem = createOptimProblem('fminunc',... 'objective',@(x)sawtoothxy(x(1),x(2)),... 'x0',[100,-50],'options',... optimoptions(@fminunc,'Algorithm','quasi-newton','Display','off'));

Validate the problem structure by running it:

[x,fval] = fminunc(problem) x = 1.7533 -111.9488 fval = 577.6960

Create a default

`MultiStart`

object:ms = MultiStart;

Run the solver for 50 iterations, recording the local minima:

% rng(1) % uncomment to obtain the same result [x,fval,eflag,output,manymins] = run(ms,problem,50) MultiStart completed some of the runs from the start points. 9 out of 50 local solver runs converged with a positive local solver exit flag. x = -142.4608 406.8030 fval = 1.2516e+03 eflag = 2 output = struct with fields: funcCount: 8586 localSolverTotal: 50 localSolverSuccess: 9 localSolverIncomplete: 41 localSolverNoSolution: 0 message: 'MultiStart completed some of the runs from the start points.↵↵9 out of 50 local solver runs converged with a positive local solver exit flag.' manymins = 1×9 GlobalOptimSolution array with properties: X Fval Exitflag Output X0

You can get different results, since

`MultiStart`

is stochastic.The solver did not find the global minimum near

`[0,0]`

. It found 10 distinct local minima.Plot the function values at the local minima:

histogram([manymins.Fval],10)

Plot the function values at the three best points:

bestf = [manymins.Fval]; histogram(bestf(1:3),10)

`MultiStart`

started `fminunc`

from
start points with components uniformly distributed between –1000
and 1000. `fminunc`

often got stuck in one of the
many local minima. `fminunc`

exceeded its iteration
limit or function evaluation limit 40 times.