% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Zoom and Scan Algorithm for unconstrained minimum
%  two dimension version
% Side constraints are not checked
%  P. Venkataraman - March 2008
%--------------------------------------------------
%  6.5.1   Scan and Zoom 
% Example 6.2: Zero Order Method
% f = 100(x2 -x1^2)^2 + (1 - x1)^2
%------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% management functions
clear  % clear all variable/information in the workspace - use CAUTION
clc    % position the cursor at the top of the screen
close   %  closes the figure window
format compact  % avoid skipping a line when writing to the command window
warning off  % don't report any warnings like divide by zero etc.

%%%   Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1Low = 0;  x1High = 6;
x2Low = 0;  x2High = 6;

%%%  stopping criteria is based on the tolerance if the
%%% objective function

%%% monitor cpu time
starttime = cputime;

%%% choose current solution at mid point of design region
%%% this will also be center of the design space at current level
x1mid = 0.5*(x1High+x1Low);
x2mid = 0.5*(x2High+x2Low);

%%% the extant of the design region in current scanning level
delx1 = 0.5*(x1High-x1Low);
delx2 = 0.5*(x2High-x2Low);

%%% tolerance on search for variables and function
tolx1 = 1.0e-08;
tolx2 = 1.0e-08;
ftol = 1.0e-8;

%%% define the limits for search in the current level
x1L = x1mid - delx1;   x1U = x1mid + delx1;
x2L = x2mid - delx2;   x2U = x2mid + delx2;

x1m = 0;  % storage for the minimum value of x1
x2m = 0;  % storage for the minimum value of x2
fm = 1.0e10;  % default min - large number

nSL = 50;  % number of scanning levels
nSP = 11; % number of scanning points per level

for ii = 1: nSL
    x1 = linspace(x1L,x1U,nSP);
    x2 = linspace(x2L,x2U,nSP);
    for i = 1:nSP
        for j = 1:nSP
            fij = 100*(x2(j)-x1(i)*x1(i))^2 + (1 - x1(i))^2;
            if fij < fm
                x1m = x1(i);
                x2m = x2(j);
                fm = fij;
            end
        end
    end
    
    %%% print the best values for the current zoom level

    fprintf('\n Example 6.2 ')
    fprintf('\n-------------')
    fprintf('\n Zoom Level (%2i) :\n',ii)
    fprintf('\nx1* = %10.7f',x1m)
    fprintf('\nx2* = %10.7f',x2m)
    fprintf('\n f* = %10.7f\n',fm)
    
    % stopping condition based on the tolerance
    % in the objective function with zoom level
    
    if abs (fij - fm) <= ftol
        fprintf('\nReached tolerance in f %7.3e \n',ftol)
        break;
    end
%  stopping based on tolerance on the variables
%     if abs(x1m - x1mid) <= tolx1
%         fprintf('\n tolerance in x1 <= %7.3e',tolx1) 
%         if abs(x2m - x2mid) <= tolx2
%             fprintf('\n tolerance in x2 <= %7.3e',tolx2)
%             break;
%         end
%     end
%%%  the tolerance on the variables was not effective

    % initialize the values for the next level.
    x1mid = x1m;
    x2mid = x2m;
    delx1 = 0.5*delx1;
    delx2 = 0.5*delx1;

    x1L = x1mid - delx1;   x1U = x1mid + delx1;
    x2L = x2mid - delx2;   x2U = x2mid + delx2;
end

%%%  Hessian Calculation
syms a b
 ff = 100*(b-a^2)^2 +(1-a)^2;

%%% gradient
dfda = diff(ff,a)
dfdb = diff(ff,b)

Hess = [diff(dfda,a) diff(dfda,b);diff(dfdb,a) diff(dfdb,b)]

gradv = [subs(dfda,{a b},{x1m x2m});subs(dfda,{a b},{x1m x2m})]

Hessv = subs(Hess,{a b},{x1m x2m})

determinnt = det(Hessv)  % determinant of Hessv
InvHessv = inv(Hessv)

totaltime = cputime - starttime;
fprintf('\nTotal time (s)= %7.4f \n\n',totaltime)


