%% Minimization with Bound Constraints and Banded Preconditioner
% This example show how to solve a nonlinear problem with bounds using the fmincon trust-region-reflective algorithm. This algorithm has provisions for added efficiency when the problem is sparse, has an analytic gradient, and has known structure such as its Hessian pattern.
%% Objective Function With Gradient
% For a given n that is a positive multiple of 4, the objective function is
% \[f(x) = 1 + \sum\limits_{i = 1}^n {{{\left| {\left( {3 - 2{x_i}} \right){x_i} - {x_{i - 1}} - {x_{i + 1}} + 1} \right|}^p}}  + \sum\limits_{i = 1}^{n/2} {{{\left| {{x_i} + {x_{i + n/2}}} \right|}^p}} ,\]
% where p = 7/3, x_0 = 0, and x_{n+1} = 0. The tbroyfg function at the end of this example implements the objective function, including its gradient.
% The problem has bounds: -10\le x_i\le 10 for all i. Use n = 800.

n = 800;
lb = -10*ones(n,1);
ub = -lb;

%% Hessian Pattern
% The sparsity pattern of the Hessian matrix has been predetermined and stored in the file tbroyhstr.mat. The sparsity structure for the Hessian of this problem is banded, as you can see in the following spy plot.

load tbroyhstr
spy(Hstr)

% In this plot, the center stripe is itself a five-banded matrix. The following plot shows the matrix more clearly.

spy(Hstr(1:20,1:20))

%% Problem Options
% Set options to use the trust-region-reflective algorithm. This algorithm requires you to set the SpecifyObjectiveGradient option to true.
% Also, use optimoptions to set the HessPattern option to Hstr. When a problem as large as this has obvious sparsity structure, not setting the HessPattern option requires a huge amount of unnecessary memory and computation. This is because fmincon attempts to use finite differencing on a full Hessian matrix of 640,000 nonzero entries.

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

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

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

% There are no linear or nonlinear constraints, so set those parameters to [].

A = [];
b = [];
Aeq = [];
beq = [];
nonlcon = [];

% Call fmincon to solve the problem.

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

%
%% Examine Solution and Solution Process
% Examine the exit flag, objective function value, first-order optimality measure, and number of solver iterations.

disp(exitflag);
disp(fval)
disp(output.firstorderopt)
disp(output.iterations)

% fmincon did not take very many iterations to reach a solution. However, the solution has a relatively high first-order optimality measure, which is the reason that the exit flag is not the more preferable value of 1.
%
%% Improve Solution
% Try using a five-banded preconditioner instead of the default diagonal preconditioner. Using optimoptions, set the PrecondBandWidth option to 2 and solve the problem again. (The bandwidth is the number of upper (or lower) diagonals, not counting the main diagonal.)

options.PrecondBandWidth = 2;
[x2,fval2,exitflag2,output2] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
disp(exitflag2);
disp(fval2)
disp(output2.firstorderopt)
disp(output2.iterations)

% The exit flag and objective function value do not appear to change. However, the number of iterations increased, and the first-order optimality measure decreased considerably. Compute the difference in objective function value.

disp(fval2 - fval)

% The objective function value decreased by a tiny amount. The improvement in the solution is mainly an improvement in the first-order optimality measure, without much improvement in the objective function.
%
%% Helper Function
% This code creates the tbroyfg function.

function [f,grad] = tbroyfg(x,dummy)
%TBROYFG Objective function and its gradients for nonlinear minimization
% with bound constraints and banded preconditioner.
% Documentation example.

%   Copyright 1990-2008 The MathWorks, Inc.

n = length(x);  % n should be a multiple of 4
p = 7/3;
y=zeros(n,1);
i = 2:(n-1);
y(i) = abs((3-2*x(i)).*x(i) - x(i-1) - x(i+1)+1).^p;
y(n) = abs((3-2*x(n)).*x(n) - x(n-1)+1).^p;
y(1) = abs((3-2*x(1)).*x(1) - x(2)+1).^p;
j = 1:(n/2);
z = zeros(length(j),1);
z(j) = abs(x(j) + x(j+n/2)).^p;
f = 1 + sum(y) + sum(z);
%
% Evaluate the gradient.
if nargout > 1
   g = zeros(n,1);
   t = zeros(n,1);
   i = 2:(n-1);
   t(i) = (3-2*x(i)).*x(i) - x(i-1) - x(i+1) + 1;
   g(i) = p*abs(t(i)).^(p-1).*sign(t(i)).*(3-4*x(i));
   g(i-1) = g(i-1) - p*abs(t(i)).^(p-1).*sign(t(i));
   g(i+1) = g(i+1) - p*abs(t(i)).^(p-1).*sign(t(i));
   tt = (3-2*x(n)).*x(n) - x(n-1) + 1;
   g(n) = g(n) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(n));
   g(n-1) = g(n-1) - p*abs(tt).^(p-1).*sign(tt);
   tt = (3-2*x(1)).*x(1)-x(2)+1;
   g(1) = g(1) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(1));
   g(2) = g(2) - p*abs(tt).^(p-1).*sign(tt);
   j = 1:(n/2);
   t(j) = x(j) + x(j+n/2);
   g(j) = g(j) + p*abs(t(j)).^(p-1).*sign(t(j));
   jj = j + (n/2);
   g(jj) = g(jj) + p*abs(t(j)).^(p-1).*sign(t(j));
   grad = g;
end
end

% Copyright 2019 The MathWorks, Inc.
%