%% Minimization with Linear Equality Constraints, Trust-Region Reflective Algorithm
% The fmincon trust-region-reflective algorithm can minimize a nonlinear objective function subject to linear equality constraints only (no bounds or any other constraints). For example, minimize
% \[f(x) = \sum\limits_{i = 1}^{n - 1} {\left( {{{\left( {x_i^2} \right)}^{\left( {x_{i + 1}^2 + 1} \right)}} + {{\left( {x_{i + 1}^2} \right)}^{\left( {x_i^2 + 1} \right)}}} \right)} ,\]
% subject to some linear equality constraints. This example takes n = 1000.
%% Create Problem
% The browneq.mat file contains the matrices Aeq and beq, which represent the linear constraints Aeq*x = beq. The Aeq matrix has 100 rows representing 100 linear constraints (so Aeq is a 100-by-1000 matrix). Load the browneq.mat data.

load browneq.mat

% The brownfgh helper function at the end of this example implements the objective function, including its gradient and Hessian.
%% Set Options
% The trust-region-reflective algorithm requires the objective function to include a gradient. The algorithm accepts a Hessian in the objective function. Set the options to include all of the derivative information.

options = optimoptions('fmincon','Algorithm','trust-region-reflective',...
    'SpecifyObjectiveGradient',true,'HessianFcn','objective');

%% Solve Problem
% Set the initial point to –1 for odd indices and +1 for even indices.

n = 1000;
x0 = -ones(n,1);
x0(2:2:n) = 1;

% The problem has no bounds, linear inequality constraints, or nonlinear constraints, so set those parameters to [].

A = [];
b = [];
lb = [];
ub = [];
nonlcon = [];

% Call fmincon to solve the problem.

[x,fval,exitflag,output] = ...
   fmincon(@brownfgh,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

%% Examine Solution and Solution Process
% Examine the exit flag, objective function value, and constraint violation.

disp(exitflag)
disp(fval)
disp(output.constrviolation)

% The exitflag value of 3 indicates that fmincon stops because the change in the objective function value is less than the tolerance FunctionTolerance. The final objective function value is given by fval. Constraints are satisfied, as shown in output.constrviolation, which displays a very small number.
% To calculate the constraint violation yourself, execute the following code.

norm(Aeq*x-beq,Inf)

%
%% Helper Function
% The following code creates the brownfgh helper function.

function [f,g,H] = brownfgh(x)
%BROWNFGH  Nonlinear minimization problem (function, its gradients
% and Hessian).
% Documentation example.         

%   Copyright 1990-2019 The MathWorks, Inc.

% Evaluate the function.
  n = length(x);
  y = zeros(n,1);
  i = 1:(n-1);
  y(i) = (x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1);
  f = sum(y);

% Evaluate the gradient.
  if nargout > 1
     i=1:(n-1);
     g = zeros(n,1);
     g(i) = 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
             2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2);
     g(i+1) = g(i+1)+...
              2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).*log(x(i).^2)+...
              2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2));
  end

% Evaluate the (sparse, symmetric) Hessian matrix
  if nargout > 2
     v = zeros(n,1);
     i = 1:(n-1);
     v(i) = 2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+...
            4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^((x(i+1).^2)-1))+...
            2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2));
     v(i) = v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2);
     v(i+1) = v(i+1)+...
              2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+...
              4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).*((log(x(i).^2)).^2)+...
              2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2));
     v(i+1) = v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i+1).^2).^(x(i).^2-1));
     v0 = v;
     v = zeros(n-1,1);
     v(i) = 4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
            4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2)).*log(x(i).^2);
     v(i) = v(i)+ 4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2);
     v(i) = v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1);
     v1 = v;
     i = [(1:n)';(1:(n-1))'];
     j = [(1:n)';(2:n)'];
     s = [v0;2*v1];
     H = sparse(i,j,s,n,n);
     H = (H+H')/2;
  end
end

% Copyright 2019 The MathWorks, Inc.