% Chapter 6: Numerical Techniques
%            For Unconstrained Optimization
%--------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Applied Optimization with Matlab Programming
% Dr. P.Venkataraman
% Second Edition,  John Wiley
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------------------------------
% Example 6.1: Graphical Solution
%  6.1.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x1 x2 	% symbolic variables
%f = 3*(sin(0.5 +0.25*x1*x2))*cos(x1)
f = 3+(x1-1.5*x2)^2+(x2-2)^2;

% computing the derivative
dfdx1 = diff(f,x1)
dfdx2 = diff(f,x2)

% computing the Hessian
Hess = [diff(dfdx1,x1) diff(dfdx1,x2);diff(dfdx2,x1) diff(dfdx2,x2)]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Analytic solution

sol = solve(dfdx1,dfdx2)
x1s = double(sol.x1)
x2s = double(sol.x2)

%% default values for minimum
x1sol = 0;
x2sol = 0;
fsol = 1.0e08;
nosol = 0;
fprintf('\nAnalytical solution')
fprintf('\n-----------------------')
for i = 1:length(x1s)
    fv = subs(f,{x1,x2},{x1s(i),x2s(i)});
    Hessv = subs(Hess,{x1,x2},{x1s(i),x2s(i)});
    eigv = eig(Hessv);
    fprintf('\nSolution  : %3i',i)
    fprintf('\n----------------')
    fprintf('\n x1 = %7.4f    x2 = %7.4f',x1s(i),x2s(i))
    fprintf('\n f = %7.4f',fv)
    fprintf('\n eigenvalues :  %7.4f , %7.4f',eigv(1),eigv(2))
    if (x1s(i) >= 0) & (x1s(i) <= 5)
        if (x2s(i) >= 0) & (x2s(i) <= 8)
            if (eigv(1) >=0) & (eigv(2) >= 0)
                fprintf('\n This is a possible minimum')
                if fv < fsol
                    fsol = fv;
                    x1sol = x1s(i);
                    x2sol = x2s(i);
                    nosol = nosol+1;
                end
            else

            end
        end

    end
   if (i == length(x1s)) & (nosol == 0)
            fprintf('\n\n There are no analytical solutions\n\n')
        end
end


%%%%%%%   end analytical solution %%%%%%%%%%%%%%%%%%%

%%%%%%5%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Numerical solution using fsolve
%%%  Create the function m-file Example6_1.m

x0 = [2 2];
xns = fsolve(@FOC_Example6_1a,x0);
fnv = subs(f,{x1,x2},{xns(1),xns(2)});
Hessvn = subs(Hess,{x1,x2},{xns(1), xns(2)});
eign = eig(Hessvn);

%%% print information
fprintf('\nNumerical solution')
fprintf('\n-----------------------')

fprintf('\n Initial Guess:  x1 = %7.4f    x2 = %7.4f',x0(1),x0(2))
fprintf('\n Final Solution: x1 = %7.4f    x2 = %7.4f',xns(1),xns(2))
fprintf('\n f = %7.4f',fnv)
fprintf('\n eigenvalues :  %7.4f , %7.4f\n\n',eign(1),eign(2))

%%%%   end Numerical Solution %%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Graphical Solution
x1g = linspace(0,5,101);
x2g = linspace(0,5,101);
[X1 X2] = meshgrid(x1g,x2g);
%ff = 3*(sin(0.5+0.25*X2.*X1)).*cos(X1);
ff = 3+(X1-1.5*X2).^2 + (X2-2).^2;
%ff = 4*X.*(180-X)/(40500 - X.*(180-X))

minf = min(min(ff));
maxf = max(max(ff));
delf = maxf - minf;
cvlow = minf + 0.0000005*delf;
cvhigh = maxf - 0.05*delf;
cv = linspace(cvlow,cvhigh,40);
[c1 h1] = contourf(x1g,x2g,ff,cv);
clabel(c1);
colorbar
grid

xlabel('x_1')
ylabel('x_2')
title('Example 6.1 - Graphical Solution')
