Find Global or Multiple Local Minima

Function to Optimize

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

g(r)=(sin(r)sin(2r)2+sin(3r)3sin(4r)4+4)r2r+1h(t)=2+cos(t)+cos(2t12)2.

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.

 Code for Generating the Figure

Single Global Minimum Via GlobalSearch

  1. 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
  2. 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.

  3. Add an artificial bound, and validate the problem structure by running fmincon:

    problem.lb = -Inf([2,1]);
    [x fval] = fmincon(problem)
    
    x =
    
       45.6965 -107.6645
    
    
    fval =
    
      555.6941
  4. Create the GlobalSearch object, and set iterative display:

    gs = GlobalSearch('Display','iter');
  5. 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.7                                   555.7            0    Initial Point
         200     1479   1.547e-15                               1.547e-15            1    Stage 1 Local
         300     1580   1.547e-15     5.858e+04        1.074                              Stage 2 Search
         400     1680   1.547e-15      1.84e+05         4.16                              Stage 2 Search
         500     1780   1.547e-15     2.683e+04        11.84                              Stage 2 Search
         600     1880   1.547e-15     1.122e+04        30.95                              Stage 2 Search
         700     1980   1.547e-15     1.353e+04        65.25                              Stage 2 Search
         800     2080   1.547e-15     6.249e+04        163.8                              Stage 2 Search
         900     2180   1.547e-15     4.119e+04        409.2                              Stage 2 Search
         950     2372   1.547e-15           477        589.7          387            2    Stage 2 Local
         952     2436   1.547e-15         368.4          477        250.7            2    Stage 2 Local
        1000     2484   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].

Multiple Local Minima Via MultiStart

  1. 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
  2. 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 in the Optimization Toolbox™ documentation.

    • The default fminunc algorithm requires a gradient; see Choosing the Algorithm in the Optimization Toolbox documentation. 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'));
  3. Validate the problem structure by running it:

    [x fval] = fminunc(problem)
    
    x =
        8.4420 -110.2602
    
    fval =
      435.2573
  4. Create a default MultiStart object:

    ms = MultiStart;
  5. Run the solver for 50 iterations, recording the local minima:

    [x fval eflag output manymins] = run(ms,problem,50)
    
    MultiStart completed some of the runs from the start points.
    
    16 out of 50 local solver runs converged with a positive 
    local solver exit flag.
    
    x =
     -379.3434  559.6154
    
    fval =
      1.7590e+003
    
    eflag =
         2
    
    output = 
                         funcCount: 7803
                  localSolverTotal: 50
                localSolverSuccess: 16
             localSolverIncomplete: 34
             localSolverNoSolution: 0
                           message: 'MultiStart completed some of the runs from ...'
    manymins = 
      1x16 GlobalOptimSolution
    
      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 16 distinct local minima.

  6. Plot the function values at the local minima:

    hist([manymins.Fval])

    Plot the function values at the three best points:

    bestf = [manymins.Fval];
    hist(bestf(1:3))

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 34 times.

Was this topic helpful?