{"id":785,"date":"2023-01-04T06:31:01","date_gmt":"2023-01-04T11:31:01","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=785"},"modified":"2023-01-04T06:31:01","modified_gmt":"2023-01-04T11:31:01","slug":"which-matlab-optimization-functions-can-solve-my-problem","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2023\/01\/04\/which-matlab-optimization-functions-can-solve-my-problem\/","title":{"rendered":"Which MATLAB Optimization functions can solve my problem?"},"content":{"rendered":"
The <\/span>solvers<\/span><\/a> function from <\/span>Optimization toolbox<\/span><\/a> is one of my favourite enhancements of R2022b because it helps improve my knowledge of which algorithms can solve my problems. <\/span><\/div>

Solving Rosenbrock with the default algorithm<\/span><\/h3>
Lets take the classic problem of minimising <\/span>Rosenbrock's function<\/span><\/a> as an example. <\/span><\/div>
$ f(x)=100\\left ( x_2 - x_1^2\\right )^2 + (1-x_1)^2 $<\/span><\/div>
The minimum is at (1,1) which is inside a long, narrow and flat valley. Many Optimization solvers find the valley rather easily but struggle with finding the global mimimum within it.<\/span><\/div>
rosen = @(x,y) 100*(y - x.^2).^2 + (1 - x).^2;<\/span><\/span><\/div><\/div>
fcontour(rosen,[-0.5,2.0],LevelStep=10)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
scatter(1,1,300,<\/span>'r.'<\/span>);<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
title(<\/span>\"Contour plot of Rosenbrock's function. Global minimum shown in red\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Using MATLAB's <\/span>problem based Optimization workflow<\/span><\/a>, we can solve this as follows<\/span><\/div>
% Set up a couple of Optimization variables and apply boundaries<\/span><\/span><\/div><\/div>
x = optimvar(<\/span>\"x\"<\/span>,LowerBound=-3,UpperBound=3);<\/span><\/span><\/div><\/div>
y = optimvar(<\/span>\"y\"<\/span>,LowerBound = 0,UpperBound=9);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
%Define the objective function<\/span><\/span><\/div><\/div>
obj = 100*(y - x^2)^2 + (1 - x)^2;<\/span><\/span><\/div><\/div>
%Set up the Optimization problem<\/span><\/span><\/div><\/div>
prob = optimproblem(Objective=obj);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
%choose a starting guess<\/span><\/span><\/div><\/div>
x0.x = -2.1;<\/span><\/span><\/div><\/div>
x0.y = 2.2;<\/span><\/span><\/div><\/div>
% solve the problem<\/span><\/span><\/div><\/div>
[sol,fval,exitflag,solution] = solve(prob,x0)<\/span><\/span><\/div>
Solving problem using lsqnonlin.\r\n\r\nLocal minimum possible.\r\n\r\nlsqnonlin stopped because the final change in the sum of squares relative to \r\nits initial value is less than the value of the function tolerance.\r\n\r\n<stopping criteria details><\/div><\/div>
sol = struct with fields:<\/span><\/span><\/div>
x: 1.0000\r\n y: 1.0000\r\n<\/div><\/div><\/div>
fval = 4.9503e-16<\/div>
exitflag = <\/span><\/div>
FunctionChangeBelowTolerance\r\n<\/div><\/div><\/div>
solution = struct with fields:<\/span><\/span><\/div>
firstorderopt: 1.2147e-06\r\n iterations: 26\r\n funcCount: 27\r\n cgiterations: 0\r\n algorithm: 'trust-region-reflective'\r\n stepsize: 3.4494e-05\r\n message: 'Local minimum possible.\u21b5\u21b5lsqnonlin stopped because the final change in the sum of squares relative to \u21b5its initial value is less than the value of the function tolerance.\u21b5\u21b5<stopping criteria details>\u21b5\u21b5Optimization stopped because the relative sum of squares (r) is changing\u21b5by less than options.FunctionTolerance = 1.000000e-06.'\r\n objectivederivative: \"forward-AD\"\r\n solver: 'lsqnonlin'\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
MATLAB analysed the problem and decided to use <\/span>lsqnonlin<\/span><\/a> as the solver. This turned out rather well for us as it went on to find the minimum without any difficulty. There were only 27 function evaluations.<\/span><\/div>

What other algorithms could solve this problem?<\/span><\/h3>
You might ask yourself, however, out of all of the Optimization solvers in both Optimization Toolbox and Global Optimization Toolbox, which of them are capable of solving this problem as well? This is what the <\/span>solvers<\/span><\/a> function does for us. Give it an Optimization problem and it will return the default solver that MATLAB would choose for that problem.<\/span><\/div>
[autosolver,allsolvers] = solvers(prob);<\/span><\/span><\/div><\/div>
fprintf(<\/span>\"The default solver for this problem is %s\\n\"<\/span>,autosolver)<\/span><\/span><\/div>
The default solver for this problem is lsqnonlin<\/div><\/div><\/div><\/div><\/div>
It also returns all other solvers that could tackle the input problem.<\/span><\/div>
allsolvers'<\/span><\/span><\/div>
ans = <\/span>10\u00d71 string<\/span><\/div><\/div>
\"lsqnonlin\"
\"lsqcurvefit\"
\"fmincon\"
\"ga\"
\"patternsea\u2026
\"surrogateopt\"
\"particlesw\u2026
\"simulannea\u2026
\"gamultiobj\"
\"paretosearch\"
<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Let's explicity instruct MATLAB to use one of these:<\/span> <\/span>fmincon<\/span><\/a><\/div>
[sol,fval,exitflag,solution] = solve(prob,x0,Solver=<\/span>\"fmincon\"<\/span>)<\/span><\/span><\/div>
Solving problem using fmincon.\r\n\r\nLocal minimum found that satisfies the constraints.\r\n\r\nOptimization completed because the objective function is non-decreasing in \r\nfeasible directions, to within the value of the optimality tolerance,\r\nand constraints are satisfied to within the value of the constraint tolerance.\r\n\r\n<stopping criteria details><\/div><\/div>
sol = struct with fields:<\/span><\/span><\/div>
x: 1.0000\r\n y: 1.0000\r\n<\/div><\/div><\/div>
fval = 1.9355e-13<\/div>
exitflag = <\/span><\/div>
OptimalSolution\r\n<\/div><\/div><\/div>
solution = struct with fields:<\/span><\/span><\/div>
iterations: 46\r\n funcCount: 71\r\n constrviolation: 0\r\n stepsize: 7.3444e-05\r\n algorithm: 'interior-point'\r\n firstorderopt: 9.0385e-07\r\n cgiterations: 14\r\n message: 'Local minimum found that satisfies the constraints.\u21b5\u21b5Optimization completed because the objective function is non-decreasing in \u21b5feasible directions, to within the value of the optimality tolerance,\u21b5and constraints are satisfied to within the value of the constraint tolerance.\u21b5\u21b5<stopping criteria details>\u21b5\u21b5Optimization completed: The relative first-order optimality measure, 9.038468e-07,\u21b5is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint\u21b5violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'\r\n bestfeasible: [1\u00d71 struct]\r\n objectivederivative: \"forward-AD\"\r\n constraintderivative: \"auto\"\r\n solver: 'fmincon'\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
So <\/span>fmincon<\/span> can and did solve the problem but it had to work much harder with 71 function evaluations. Just because the <\/span>solvers<\/span> function determines that an algorithm is <\/span>capable<\/span> of tackling a problem doesn't necessarily mean its a good fit. An example of a really bad fit for this problem is <\/span>Simulated Annealing<\/span><\/a>. <\/span><\/div>
[sol,fval,exitflag,solution]= solve(prob,x0,Solver=<\/span>\"simulannealbnd\"<\/span>)<\/span><\/span><\/div>
Solving problem using simulannealbnd.\r\nOptimization terminated: change in best function value less than options.FunctionTolerance.<\/div><\/div>
sol = struct with fields:<\/span><\/span><\/div>
x: 2.8612\r\n y: 8.1900\r\n<\/div><\/div><\/div>
fval = 3.4653<\/div>
exitflag = <\/span><\/div>
SolverConvergedSuccessfully\r\n<\/div><\/div><\/div>
solution = struct with fields:<\/span><\/span><\/div>
iterations: 2247\r\n funccount: 2262\r\n message: 'Optimization terminated: change in best function value less than options.FunctionTolerance.'\r\n rngstate: [1\u00d71 struct]\r\n problemtype: 'boundconstraints'\r\n temperature: [2\u00d71 double]\r\n totaltime: 0.2231\r\n solver: 'simulannealbnd'\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
Yikes! That's rather far from the actual solution and we've already burned though <\/span>over a thousand<\/span> function evaluations. Perhaps simulated annealing isn't such a good idea for this function after all. <\/span><\/div>
Let's try every solver that's valid for this problem using default settings for each one<\/span><\/div>
variableNames = [<\/span>\"x solution\"<\/span>,<\/span>\"y solution\"<\/span>,<\/span>\"objective value\"<\/span>,<\/span>\"function evaluations\"<\/span>];<\/span><\/span><\/div><\/div>
results = array2table(NaN(numel(allsolvers), numel(variableNames)));<\/span><\/span><\/div><\/div>
results.Properties.VariableNames = variableNames;<\/span><\/span><\/div><\/div>
results.Properties.RowNames = allsolvers;<\/span><\/span><\/div><\/div>
 <\/div><\/div>
for <\/span>solverNum = 1:numel(allsolvers)<\/span><\/span><\/div><\/div>
solver = allsolvers(solverNum);<\/span><\/span><\/div><\/div>
<\/span>% I want each solver to be quiet while it works. <\/span><\/span><\/div><\/div>
options=optimoptions(solver,Display=<\/span>\"None\"<\/span>);<\/span><\/span><\/div><\/div>
<\/span>if <\/span>solver==<\/span>\"surrogateopt\"<\/span><\/span><\/div><\/div>
options=optimoptions(options,PlotFcn=[]);<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
[sol,fval,exitflag,solution] = solve(prob,x0,Solver=solver,Options=options);<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
<\/span>%Some solvers have solution.funcCount, others have solution.funccount<\/span><\/span><\/div><\/div>
<\/span>if <\/span>isfield(solution,<\/span>'funcCount'<\/span>)<\/span><\/span><\/div><\/div>
funcount = solution.funcCount;<\/span><\/span><\/div><\/div>
<\/span>else<\/span><\/span><\/div><\/div>
funcount = solution.funccount;<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
results.(<\/span>\"x solution\"<\/span>)(solverNum) = sol.x;<\/span><\/span><\/div><\/div>
results.(<\/span>\"y solution\"<\/span>)(solverNum) = sol.y;<\/span><\/span><\/div><\/div>
results.(<\/span>\"objective value\"<\/span>)(solverNum) = fval;<\/span><\/span><\/div><\/div>
results.(<\/span>\"function evaluations\"<\/span>)(solverNum) = funcount;<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
results<\/span><\/span><\/div>
results = <\/span>10\u00d74 table <\/span><\/div>
 <\/span><\/th>x solution<\/span><\/th>y solution<\/span><\/th>objective value<\/span><\/th>function evaluations<\/span><\/th><\/tr><\/thead>
1 lsqnonlin<\/span><\/th>1<\/span><\/td>1<\/span><\/td>0<\/span><\/td>27<\/span><\/td><\/tr>
2 lsqcurvefit<\/span><\/th>1<\/span><\/td>1<\/span><\/td>0<\/span><\/td>27<\/span><\/td><\/tr>
3 fmincon<\/span><\/th>1<\/span><\/td>1<\/span><\/td>0<\/span><\/td>71<\/span><\/td><\/tr>
4 ga<\/span><\/th>0.9443<\/span><\/td>0.8904<\/span><\/td>0.0033<\/span><\/td>9453<\/span><\/td><\/tr>
5 patternsearch<\/span><\/th>1.8658<\/span><\/td>3.4852<\/span><\/td>0.7511<\/span><\/td>734<\/span><\/td><\/tr>
6 surrogateopt<\/span><\/th>0.9885<\/span><\/td>0.9768<\/span><\/td>0.0001<\/span><\/td>200<\/span><\/td><\/tr>
7 particleswarm<\/span><\/th>0.9996<\/span><\/td>0.9993<\/span><\/td>0<\/span><\/td>5200<\/span><\/td><\/tr>
8 simulannealbnd<\/span><\/th>2.3423<\/span><\/td>5.4883<\/span><\/td>1.8021<\/span><\/td>2151<\/span><\/td><\/tr>
9 gamultiobj<\/span><\/th>1.0121<\/span><\/td>1.0244<\/span><\/td>0.0001<\/span><\/td>5100<\/span><\/td><\/tr>
10 paretosearch<\/span><\/th>0.9497<\/span><\/td>0.9088<\/span><\/td>0.0074<\/span><\/td>9000<\/span><\/td><\/tr><\/tbody><\/table><\/div><\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
It seems that MATLAB's default choice for this problem, <\/span>lsqnonlin<\/span>, was a good one. Indeed it was the best one. That's not always the case though.<\/span><\/div>

Which algorithm is best for Rastigrin's function?<\/span><\/h3>
Consider <\/span>Rastigrin's function<\/span><\/a>. This is another classic test function in numerical optimization with many local minima and a global minimum at (0,0):<\/span><\/div>
ras = @(x, y) 20 + x.^2 + y.^2 - 10*(cos(2*pi*x) + cos(2*pi*y));<\/span><\/span><\/div><\/div><\/div>
Plotting it shows just how many local minima there are<\/span><\/div>
fsurf(ras,[-3 3],<\/span>\"ShowContours\"<\/span>,<\/span>\"on\"<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>\"rastriginsfcn(x,y)\"<\/span>)<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"x\"<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>\"y\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
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 [2,3], which is far from the global minimum. It also uses finite bounds of -70 to 130 for each variable<\/span><\/div>
x = optimvar(<\/span>\"x\"<\/span>,<\/span>\"LowerBound\"<\/span>,-70,<\/span>\"UpperBound\"<\/span>,130);<\/span><\/span><\/div><\/div>
y = optimvar(<\/span>\"y\"<\/span>,<\/span>\"LowerBound\"<\/span>,-70,<\/span>\"UpperBound\"<\/span>,130);<\/span><\/span><\/div><\/div>
prob = optimproblem(<\/span>\"Objective\"<\/span>,ras(x,y));<\/span><\/span><\/div><\/div>
x0.x = 2;<\/span><\/span><\/div><\/div>
x0.y = 3;<\/span><\/span><\/div><\/div><\/div>
This time, MATLAB chooses <\/span>fmincon<\/span><\/a> as the default solver which gets stuck in one of the local minima <\/span><\/div>
solve(prob,x0)<\/span><\/span><\/div>
Solving problem using fmincon.\r\n\r\nLocal minimum found that satisfies the constraints.\r\n\r\nOptimization completed because the objective function is non-decreasing in \r\nfeasible directions, to within the value of the optimality tolerance,\r\nand constraints are satisfied to within the value of the constraint tolerance.\r\n\r\n<stopping criteria details><\/div><\/div>
ans = struct with fields:<\/span><\/span><\/div>
x: -1.9899\r\n y: -2.9849\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
As before, we use <\/span>solvers<\/span><\/a> to get the list of everything that could tackle this problem<\/span><\/div>
[autosolver,allsolvers] = solvers(prob)<\/span><\/span><\/div>
autosolver = <\/span>\"fmincon\"<\/div><\/div><\/div>
allsolvers = <\/span>1\u00d78 string<\/span><\/div><\/div>
\"fmincon\" \"ga\" \"patternse\u2026 \"surrogate\u2026 \"particles\u2026 \"simulannea\u2026 \"gamultiobj\" \"paretosea\u2026
<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Sweeping across all valid solvers again and we see that, this time, global Optimization functions such as <\/span>patternsearch<\/span><\/a>, <\/span>particle swarm<\/span><\/a> and <\/span>simulated annealing<\/span><\/a> do much better.<\/span><\/div>
variableNames = [<\/span>\"x\"<\/span>,<\/span>\"y\"<\/span>,<\/span>\"objective value\"<\/span>,<\/span>\"function evaluations\"<\/span>];<\/span><\/span><\/div><\/div>
results = array2table(NaN(numel(allsolvers), numel(variableNames)));<\/span><\/span><\/div><\/div>
results.Properties.VariableNames = variableNames;<\/span><\/span><\/div><\/div>
results.Properties.RowNames = allsolvers;<\/span><\/span><\/div><\/div>
 <\/div><\/div>
for <\/span>solverNum = 1:numel(allsolvers)<\/span><\/span><\/div><\/div>
solver = allsolvers(solverNum);<\/span><\/span><\/div><\/div>
<\/span>% I want each solver to be quiet while it works. <\/span><\/span><\/div><\/div>
options=optimoptions(solver,Display=<\/span>\"None\"<\/span>);<\/span><\/span><\/div><\/div>
<\/span>if <\/span>solver==<\/span>\"surrogateopt\"<\/span><\/span><\/div><\/div>
options=optimoptions(options,PlotFcn=[]);<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
[sol,fval,exitflag,solution] = solve(prob,x0,Solver=solver,Options=options);<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
<\/span>%Some solvers have solution.funcCount, others have solution.funccount<\/span><\/span><\/div><\/div>
<\/span>if <\/span>isfield(solution,<\/span>'funcCount'<\/span>)<\/span><\/span><\/div><\/div>
funcount = solution.funcCount;<\/span><\/span><\/div><\/div>
<\/span>else<\/span><\/span><\/div><\/div>
funcount = solution.funccount;<\/span><\/span><\/div><\/div>
<\/span>end<\/span><\/span><\/div><\/div>
results.(<\/span>\"x\"<\/span>)(solverNum) = sol.x;<\/span><\/span><\/div><\/div>
results.(<\/span>\"y\"<\/span>)(solverNum) = sol.y;<\/span><\/span><\/div><\/div>
results.(<\/span>\"objective value\"<\/span>)(solverNum) = fval;<\/span><\/span><\/div><\/div>
results.(<\/span>\"function evaluations\"<\/span>)(solverNum) = funcount;<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
results<\/span><\/span><\/div>
results = <\/span>8\u00d74 table <\/span><\/div>
 <\/span><\/th>x<\/span><\/th>y<\/span><\/th>objective value<\/span><\/th>function evaluations<\/span><\/th><\/tr><\/thead>
1 fmincon<\/span><\/th>-1.9899<\/span><\/td>-2.9849<\/span><\/td>12.9344<\/span><\/td>39<\/span><\/td><\/tr>
2 ga<\/span><\/th>1.9899<\/span><\/td>2.9849<\/span><\/td>12.9344<\/span><\/td>3343<\/span><\/td><\/tr>
3 patternsearch<\/span><\/th>0<\/span><\/td>0<\/span><\/td>0<\/span><\/td>99<\/span><\/td><\/tr>
4 surrogateopt<\/span><\/th>-0.9982<\/span><\/td>0.0094<\/span><\/td>1.0146<\/span><\/td>200<\/span><\/td><\/tr>
5 particleswarm<\/span><\/th>0.9950<\/span><\/td>0<\/span><\/td>0.9950<\/span><\/td>1500<\/span><\/td><\/tr>
6 simulannealbnd<\/span><\/th>-0<\/span><\/td>0<\/span><\/td>0<\/span><\/td>2864<\/span><\/td><\/tr>
7 gamultiobj<\/span><\/th>-0<\/span><\/td>0<\/span><\/td>0<\/span><\/td>5100<\/span><\/td><\/tr>
8 paretosearch<\/span><\/th>-0.9950<\/span><\/td>0<\/span><\/td>0.9950<\/span><\/td>7153<\/span><\/td><\/tr><\/tbody><\/table><\/div><\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
It's worth bearing in mind that these 'solver sweeps' only consider the default settings for each solver. Each solver has a range of options that can improve performance and solution quality that haven't been considered here but we have made a good start in exploring the algorithms we can bring to bear on our problem.<\/span><\/div>
To learn more about the Optimization algorithms available in MATLAB, take a look at the documentation. <\/span><\/div>