%% Output Functions for Optimization Toolbox™
%% What Is an Output Function?
% For some problems, you might want output from an optimization algorithm at each iteration. For example, you might want to find the sequence of points that the algorithm computes and plot those points. To do this, create an output function that the optimization function calls at each iteration. See Output Function Syntax for details and syntax.
% This section shows the solver-based approach to output functions. For the problem-based approach, see Output Function for Problem-Based Optimization.
% Generally, the solvers that can employ an output function are the ones that can take nonlinear functions as inputs. You can determine which solvers can have an output function by looking in the Options section of function reference pages.
%% Example: Use an Output Function
% This example shows how to use an output function to monitor the fmincon solution process for solving a constrained nonlinear optimization problem. The output function does the following at the end of each fmincon iteration:
% Plot the current point.
% Store the current point and its corresponding objective function value in a variable called history, and store the current search direction in a variable called searchdir. The search direction is a vector that points in the direction from the current point to the next one.
% Additionally, to make the history available outside of the fmincon function, perform the optimization inside a nested function that calls fmincon and returns the output function variables. For more information about this method of passing information, see Passing Extra Parameters. The runfmincon function at the end of this example contains the nested function call.
%% Objective and Constraint Functions
% The problem is to minimize the function
% f(x) = \exp(x_1)\left (4x_1^2 + 2x_2^2 + 4x_1x_2 + 2x_2 +1\right )
% subject to the nonlinear inequality constraints
% x_1 + x_2 -x_1x_2 \ge 3/2\\
% x_1x_2 \ge -10.
% The objfun function nested in runfmincon implements the objective function. The confun function nested in runfmincon implements the constraint function.
%% Call Solver
% To obtain the solution to the problem and see the history of the fmincon iterations, call the runfmincon function.

[xsol,fval,history,searchdir] = runfmincon;

% The output function creates a plot of the points that fmincon evaluates. Each point is labeled by its iteration number. The optimal point occurs at the eighth iteration. The last two points in the sequence are so close that they overlap.
% The output history is a structure that contains two fields.

disp(history)

% The fval field in history contains the objective function values corresponding to the sequence of points fmincon computes.

disp(history.fval)

% These are the same values displayed in the iterative output in the column with header f(x).
% The x field of history contains the sequence of points that fmincon computes.

disp(history.x)

% The second output argument, searchdir, contains the search directions for fmincon at each iteration. The search direction is a vector pointing from the point computed at the current iteration to the point computed at the next iteration.

disp(searchdir)

%
%% Helper Functions
% The following code creates the runfmincon function, containing the outfun output function, objfun objective function, and confun nonlinear constraint function as nested functions,

function [xsol,fval,history,searchdir] = runfmincon
 
% Set up shared variables with OUTFUN
history.x = [];
history.fval = [];
searchdir = [];
 
% call optimization
x0 = [-1 1];
options = optimoptions(@fmincon,'OutputFcn',@outfun,... 
    'Display','iter','Algorithm','active-set');
[xsol,fval] = fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options);
 
 function stop = outfun(x,optimValues,state)
     stop = false;
 
     switch state
         case 'init'
             hold on
         case 'iter'
         % Concatenate current point and objective function
         % value with history. x must be a row vector.
           history.fval = [history.fval; optimValues.fval];
           history.x = [history.x; x];
         % Concatenate current search direction with 
         % searchdir.
           searchdir = [searchdir;... 
                        optimValues.searchdirection'];
           plot(x(1),x(2),'o');
         % Label points with iteration number and add title.
         % Add .15 to x(1) to separate label from plotted 'o'
           text(x(1)+.15,x(2),... 
                num2str(optimValues.iteration));
           title('Sequence of Points Computed by fmincon');
         case 'done'
             hold off
         otherwise
     end
 end
 
 function f = objfun(x)
     f = exp(x(1))*(4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) +... 
                    2*x(2) + 1);
 end
 
 function [c, ceq] = confun(x)
     % Nonlinear inequality constraints
     c = [1.5 + x(1)*x(2) - x(1) - x(2);
         -x(1)*x(2) - 10];
     % Nonlinear equality constraints
     ceq = [];
 end
end

% Copyright 2020 The MathWorks, Inc.