%% Tutorial for Optimization Toolbox™
% This tutorial includes multiple examples that show how to use two nonlinear optimization solvers, fminunc and fmincon, and how to set options. The principles outlined in this tutorial apply to the other nonlinear solvers, such as fgoalattain, fminimax, lsqnonlin, lsqcurvefit, and fsolve.
% The tutorial examples cover these tasks:
% Minimizing an objective function
% Minimizing the same function with additional parameters
% Minimizing the objective function with a constraint
% Obtaining a more efficient or accurate solution by providing gradients or a Hessian, or by changing options
%% Unconstrained Optimization Example
% Consider the problem of finding a minimum of the function
% x\exp(-(x^2+y^2))+(x^2+y^2)/20.
% Plot the function to see where it is minimized.

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
fsurf(f,[-2,2],'ShowContours','on')

%
% The plot shows that the minimum is near the point (–1/2,0).
% Usually you define the objective function as a MATLAB® file. In this case, the function is simple enough to define as an anonymous function.

fun = @(x) f(x(1),x(2));

%
% Set an initial point for finding the solution.

x0 = [-.5; 0];

%
% Set optimization options to use the fminunc default 'quasi-newton' algorithm. This step ensures that the tutorial works the same in every MATLAB version.

options = optimoptions('fminunc','Algorithm','quasi-newton');

%
% View the iterations as the solver performs its calculations.

options.Display = 'iter';

%
% Call fminunc, an unconstrained nonlinear minimizer.

[x, fval, exitflag, output] = fminunc(fun,x0,options);

%
% Display the solution found by the solver.

uncx = x

%
% View the function value at the solution.

uncf = fval

%
% The examples use the number of function evaluations as a measure of efficiency. View the total number of function evaluations.

output.funcCount

%
%% Unconstrained Optimization Example with Additional Parameters
% Next, pass extra parameters as additional arguments to the objective function, first by using a MATLAB file, and then by using a nested function.
% Consider the objective function from the previous example.
% f(x,y) = x\exp(-(x^2+y^2))+(x^2+y^2)/20.
% Parameterize the function with (a,b,c) as follows:
% f(x,y,a,b,c) = (x-a)\exp(-((x-a)^2+(y-b)^2))+((x-a)^2+(y-b)^2)/c.
% This function is a shifted and scaled version of the original objective function.
%% MATLAB File Function
% Consider a MATLAB file objective function named bowlpeakfun defined as follows.

type bowlpeakfun

%
% Define the parameters.

a = 2;
b = 3;
c = 10;

%
% Create an anonymous function handle to the MATLAB file.

f = @(x)bowlpeakfun(x,a,b,c)

%
% Call fminunc to find the minimum.

x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');
[x, fval] = fminunc(f,x0,options)

%
%% Nested Function
% Consider the nestedbowlpeak function, which implements the objective as a nested function.

type nestedbowlpeak

% The parameters (a,b,c) are visible to the nested objective function nestedfun. The outer function, nestedbowlpeak, calls fminunc and passes the objective function, nestedfun.
%
% Define the parameters, initial guess, and options:

a = 2;
b = 3;
c = 10;
x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');

%
% Run the optimization:

[x,fval] =  nestedbowlpeak(a,b,c,x0,options)

% Both approaches produce the same answers, so you can use the one you find most convenient.
%% Constrained Optimization Example: Inequalities
% Consider the previous problem with a constraint:
% \mbox{minimize }x\exp(-(x^2+y^2))+(x^2+y^2)/20,
% \mbox{subject to }xy/2 + (x+2)^2 + (y-2)^2/2 \le 2.
% The constraint set is the interior of a tilted ellipse. View the contours of the objective function plotted together with the tilted ellipse.

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
g = @(x,y) x.*y/2+(x+2).^2+(y-2).^2/2-2;
fimplicit(g)
axis([-6 0 -1 7])
hold on
fcontour(f)
plot(-.9727,.4685,'ro');
legend('constraint','f contours','minimum');
hold off

%
% The plot shows that the lowest value of the objective function within the ellipse occurs near the lower-right part of the ellipse. Before calculating the plotted minimum, make a guess at the solution.

x0 = [-2 1];

%
% Set optimization options to use the interior-point algorithm and display the results at each iteration.

options = optimoptions('fmincon','Algorithm','interior-point','Display','iter');

%
% Solvers require that nonlinear constraint functions give two outputs, one for nonlinear inequalities and one for nonlinear equalities. To give both outputs, write the constraint using the deal function.

gfun = @(x) deal(g(x(1),x(2)),[]);

%
% Call the nonlinear constrained solver. The problem has no linear equalities or inequalities or bounds, so pass [ ] for those arguments.

[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],gfun,options);

%
% Display the solution found by the solver.

x

%
% View the function value at the solution.

fval

%
% View the total number of function evaluations.

Fevals = output.funcCount

%
% The inequality constraint is satisfied at the solution.

[c, ceq] = gfun(x)

%
% Because c(x) is close to 0, the constraint is active, meaning it affects the solution. Recall the unconstrained solution.

uncx

%
% Recall the unconstrained objective function.

uncf

%
% See how much the constraint moved the solution and increased the objective.

fval-uncf

%
%% Constrained Optimization Example: User-Supplied Gradients
% You can solve optimization problems more efficiently and accurately by supplying gradients. This example, like the previous one, solves the inequality-constrained problem
% \mbox{minimize }x\exp(-(x^2+y^2))+(x^2+y^2)/20,
% \mbox{subject to }xy/2 + (x+2)^2 + (y-2)^2/2 \le 2.
% To provide the gradient of f(x) to fmincon, write the objective function in the form of a MATLAB file.

type onehump

%
% The constraint and its gradient are contained in the MATLAB file tiltellipse.

type tiltellipse

%
% Set an initial point for finding the solution.

x0 = [-2; 1];

%
% Set optimization options to use the same algorithm as in the previous example for comparison purposes.

options = optimoptions('fmincon','Algorithm','interior-point');

%
% Set options to use the gradient information in the objective and constraint functions. Note: these options must be turned on or the gradient information will be ignored.

options = optimoptions(options,...
    'SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true);

%
% Because fmincon does not need to estimate gradients using finite differences, the solver should have fewer function counts. Set options to display the results at each iteration.

options.Display = 'iter';

%
% Call the solver.

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);

%
% fmincon estimated gradients well in the previous example, so the iterations in this example are similar.
% Display the solution found by the solver.

xold = x

%
% View the function value at the solution.

minfval = fval

%
% View the total number of function evaluations.

Fgradevals = output.funcCount

%
% Compare this number to the number of function evaluations without gradients.

Fevals

%
%% Constrained Optimization Example: Changing the Default Termination Tolerances
% This example continues to use gradients and solves the same constrained problem
% \mbox{minimize }x\exp(-(x^2+y^2))+(x^2+y^2)/20,
% \mbox{subject to }xy/2 + (x+2)^2 + (y-2)^2/2 \le 2.
% In this case, you achieve a more accurate solution by overriding the default termination criteria (options.StepTolerance and options.OptimalityTolerance). The default values for the fmincon interior-point algorithm are options.StepTolerance = 1e-10 and options.OptimalityTolerance = 1e-6.
% Override these two default termination criteria.

options = optimoptions(options,...
    'StepTolerance',1e-15,...
    'OptimalityTolerance',1e-8);

%
% Call the solver.

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);

%
% To see the difference made by the new tolerances more accurately, display more decimals in the solution.

format long

%
% Display the solution found by the solver.

x

%
% Compare these values to the values in the previous example.

xold

%
% Determine the change in values.

x - xold

%
% View the function value at the solution.

fval

%
% See how much the solution improved.

fval - minfval

%
% The answer is negative because the new solution is smaller.
% View the total number of function evaluations.

output.funcCount

%
% Compare this number to the number of function evaluations in the example solved with user-provided gradients and the default tolerances.

Fgradevals

%
%% Constrained Optimization Example: User-Supplied Hessian
% If you supply a Hessian in addition to a gradient, solvers are even more accurate and efficient.
% The fmincon interior-point algorithm takes a Hessian matrix as a separate function (not part of the objective function). The Hessian function H(x,lambda) evaluates the Hessian of the Lagrangian; see Hessian for fmincon interior-point algorithm.
% Solvers calculate the values lambda.ineqnonlin and lambda.eqlin; your Hessian function tells solvers how to use these values.
% This example has one inequality constraint, so the Hessian is defined as given in the hessfordemo function.

type hessfordemo

%
% In order to use the Hessian, you need to set options appropriately.

options = optimoptions('fmincon',...
    'Algorithm','interior-point',...
    'SpecifyConstraintGradient',true,...
    'SpecifyObjectiveGradient',true,...
    'HessianFcn',@hessfordemo);

%
% The tolerances are set to their defaults, which should result in fewer function counts. Set options to display the results at each iteration.

options.Display = 'iter';

%
% Call the solver.

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);

%
% The results show fewer and different iterations.
% Display the solution found by the solver.

x

%
% View the function value at the solution.

fval

%
% View the total number of function evaluations.

output.funcCount

%
% Compare this number to the number of function evaluations in the example solved using only gradient evaluations, with the same default tolerances.

Fgradevals

% Copyright 1990-2020 The MathWorks, Inc.