%% Tutorial for the Optimization Toolbox™
% This example shows how to use two nonlinear optimization solvers and how to 
% set options. The nonlinear solvers that we use in this example are |fminunc| 
% and |fmincon|.
% 
% All the principles outlined in this example apply to the other nonlinear 
% solvers, such as |fgoalattain|, |fminimax|, |lsqnonlin|, |lsqcurvefit|, and 
% |fsolve|.
% 
% The example starts with minimizing an objective function, then proceeds 
% to minimize the same function with additional parameters. After that, the example 
% shows how to minimize the objective function when there is a constraint, and 
% finally shows how to get a more efficient and/or accurate solution by providing 
% gradients or Hessian, or by changing some 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 get an idea of 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. For now, this 
% function is simple enough to define as an anonymous function:

fun = @(x) f(x(1),x(2));
%% 
% Take a guess at the solution:

x0 = [-.5; 0];
%% 
% Set optimization options to not use |fminunc|'s default large-scale algorithm, 
% since that algorithm requires the objective function gradient to be provided:

options = optimoptions('fminunc','Algorithm','quasi-newton');
%% 
% View the iterations as the solver calculates:

options.Display = 'iter';
%% 
% Call |fminunc|, an unconstrained nonlinear minimizer:

[x, fval, exitflag, output] = fminunc(fun,x0,options);
%% 
% The solver found a solution at:

uncx = x
%% 
% The function value at the solution is:

uncf = fval
%% 
% We will use the number of function evaluations as a measure of efficiency 
% in this example. The total number of function evaluations is: 

output.funcCount
%% Unconstrained Optimization with Additional Parameters
% We will now pass extra parameters as additional arguments to the objective 
% function. We show two different ways of doing this - using a MATLAB file, or 
% using a nested function.
% 
% Consider the objective function from the previous section:
% 
% $$f(x,y) = x\exp(-(x^2+y^2))+(x^2+y^2)/20.$$
% 
% We parameterize the function with (a,b,c) in the following way:
% 
% $$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. 
% 
% Method 1: MATLAB file Function
% 
% Suppose we have a MATLAB file objective function called |bowlpeakfun| defined 
% as:

mybowlpeakfun=@(x,y)(x-a).*exp(-((x-a).^2)+(y-b).^2)+((x-a).^2+(y-b).^2)./c
%% 
% 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)
%% 
% Method 2: Nested Function
% 
% Consider the following function that implements the objective as a nested 
% function

type nestedbowlpeak
%% 
% In this method, the parameters (a,b,c) are visible to the nested objective 
% function called |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)
%% 
% You can see both methods produced identical answers, so use whichever 
% one you find most convenient.
%% Constrained Optimization Example: Inequalities
% Consider the above 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. Look at 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. We are about to 
% calculate the minimum that was just plotted. Take a guess at the solution:

x0 = [-2 1];
%% 
% Set optimization options: use the interior-point algorithm, and turn on 
% the display of 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, the second for nonlinear equalities. So we write 
% the constraint using the |deal| function to give both outputs:

gfun = @(x) deal(g(x(1),x(2)),[]);
%% 
% Call the nonlinear constrained solver. There are no linear equalities 
% or inequalities or bounds, so pass [ ] for those arguments:

[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],gfun,options);
%% 
% A solution to this problem has been found at:

x
%% 
% The function value at the solution is: 

fval
%% 
% The total number of function evaluations was: 

Fevals = output.funcCount
%% 
% The inequality constraint is satisfied at the solution.

[c, ceq] = gfun(x)
%% 
% Since c(x) is close to 0, the constraint is "active," meaning the constraint 
% affects the solution. Recall the unconstrained solution was found at

uncx
%% 
% and the unconstrained objective function was found to be

uncf
%% 
% The constraint moved the solution, and increased the objective by

fval-uncf
%% Constrained Optimization Example: User-Supplied Gradients
% Optimization problems can be solved more efficiently and accurately if gradients 
% are supplied by the user. This example shows how this may be performed. We again 
% solve 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|, we 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
%% 
% Make a guess at the solution:

x0 = [-2; 1];
%% 
% Set optimization options: we continue to use the same algorithm for comparison 
% purposes. 

options = optimoptions('fmincon','Algorithm','interior-point');
%% 
% We also 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);
%% 
% There should be fewer function counts this time, since |fmincon| does 
% not need to estimate gradients using finite differences. 

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 the current example are similar.
% 
% The solution to this problem has been found at:

xold = x
%% 
% The function value at the solution is: 

minfval = fval
%% 
% The total number of function evaluations was: 

Fgradevals = output.funcCount
%% 
% Compare this to the number of function evaluations without gradients:

Fevals
%% Changing the Default Termination Tolerances
% This time we solve 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,$$
% 
% more accurately by overriding the default termination criteria (options.StepTolerance 
% and options.OptimalityTolerance). We continue to use gradients. The default 
% values for |fmincon|'s interior-point algorithm are options.StepTolerance = 
% 1e-10, options.OptimalityTolerance = 1e-6.
% 
% Override two default termination criteria: termination tolerances on X 
% and fval.

options = optimoptions(options,...
    'StepTolerance',1e-15,...
    'OptimalityTolerance',1e-8);
%% 
% Call the solver:

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
%% 
% We now choose to see more decimals in the solution, in order to see more 
% accurately the difference that the new tolerances make.

format long
%% 
% The optimizer found a solution at:

x
%% 
% Compare this to the previous value:

xold
%% 
% The change is

x - xold
%% 
% The function value at the solution is: 

fval
%% 
% The solution improved by

fval - minfval
%% 
% (this is negative since the new solution is smaller)
% 
% The total number of function evaluations was: 

output.funcCount
%% 
% Compare this to the number of function evaluations when the problem is 
% solved with user-provided gradients but with the default tolerances:

Fgradevals
%% Constrained Optimization Example with User-Supplied Hessian
% If you give not only a gradient, but also a Hessian, solvers are even more 
% accurate and efficient.
% 
% |fmincon|'s interior-point solver takes a Hessian matrix as a separate 
% function (not part of the objective function). The Hessian function H(x,lambda) 
% should evaluate the Hessian of the Lagrangian; see the User's Guide for the 
% definition of this term.
% 
% Solvers calculate the values lambda.ineqnonlin and lambda.eqlin; your Hessian 
% function tells solvers how to use these values.
% 
% In this problem we have but one inequality constraint, so the Hessian is:

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 have been set back to the defaults. There should be fewer 
% function counts this time.

options.Display = 'iter';
%% 
% Call the solver:

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
%% 
% There were fewer, and different iterations this time.
% 
% The solution to this problem has been found at:

x
%% 
% The function value at the solution is: 

fval
%% 
% The total number of function evaluations was: 

output.funcCount
%% 
% Compare this to the number using only gradient evaluations, with the same 
% default tolerances:

Fgradevals
%% 
% _Copyright 1990-2015 The MathWorks, Inc._