In [64]:
using ForwardDiff
#Function Definition for Hock-Schittkowski Collection, problems 6 through 9
#Problem 6 ########################################
function F_6(x)
return (1-x[1])^2;
end;
#F_6_x0 = [-1.2,1.0];
F_6_x0 = [0.0,0.0];
function H_6(x)
return 10*(x[2]-x[1]^2);
end;
F_6_grd(x) = ForwardDiff.gradient(x -> F_6(x))(x);
F_6_hes(x) = ForwardDiff.hessian(x -> F_6(x))(x);
H_6_grd(x) = ForwardDiff.gradient(x -> H_6(x))(x);
H_6_hes(x) = ForwardDiff.hessian(x -> H_6(x))(x);
###################################################
#Problem 7 ########################################
function F_7(x)
return log(1.0+x[1]^2)-x[2];
end;
#F_7_x0 = [2.0,2.0];
F_7_x0 = [1.0,sqrt(3)]; #a feasible point
function H_7(x)
return (1+x[1]^2)^2 + x[2]^2 - 4.0;
end;
F_7_grd(x) = ForwardDiff.gradient(x -> F_7(x))(x);
F_7_hes(x) = ForwardDiff.hessian(x -> F_7(x))(x);
H_7_grd(x) = ForwardDiff.gradient(x -> H_7(x))(x);
H_7_hes(x) = ForwardDiff.hessian(x -> H_7(x))(x);
###################################################
#Problem 8 ########################################
function F_8(x)
return -1;
end;
#F_8_x0 = [2.0,1.0];
F_8_x0 =[0.0,0.0];
function H_8_1(x) #first constraint
return x[1]^2 + x[2]^2 - 25.0;
end;
function H_8_2(x) #second constraint
return x[1]*x[2] - 9.0;
end;
F_8_grd(x) = ForwardDiff.gradient(x -> F_8(x))(x);
F_8_hes(x) = ForwardDiff.hessian(x -> F_8(x))(x);
H_8_1_grd(x) = ForwardDiff.gradient(x -> H_8_1(x))(x);
H_8_1_hes(x) = ForwardDiff.hessian(x -> H_8_1(x))(x);
H_8_2_grd(x) = ForwardDiff.gradient(x -> H_8_2(x))(x);
H_8_2_hes(x) = ForwardDiff.hessian(x -> H_8_2(x))(x);
###################################################
#Problem 9 ########################################
function F_9(x)
return sin(pi*x[1]/12.0)*cos(pi*x[2]/12.0);
end;
F_9_x0 = [0.0,0.0]; #a feasible point
function H_9(x)
return 4*x[1] - 3*x[2];
end;
F_9_grd(x) = ForwardDiff.gradient(x -> F_9(x))(x);
F_9_hes(x) = ForwardDiff.hessian(x -> F_9(x))(x);
H_9_grd(x) = ForwardDiff.gradient(x -> H_9(x))(x);
H_9_hes(x) = ForwardDiff.hessian(x -> H_9(x))(x);
###################################################
#Test Problem #####################################
function F_1(x)
return x[1]+x[2];
end;
F_1_x0 = [0.3,0.3]; #a feasible point
function H_1(x)
return x[1]^2 + x[2]^2 - 2;
end;
F_1_grd(x) = ForwardDiff.gradient(x -> F_1(x))(x);
F_1_hes(x) = ForwardDiff.hessian(x -> F_1(x))(x);
H_1_grd(x) = ForwardDiff.gradient(x -> H_1(x))(x);
H_1_hes(x) = ForwardDiff.hessian(x -> H_1(x))(x);
###################################################
Out[64]:
In [2]:
function lagrangian(f, h, x_k0, lambda, rho)
#take in the function definitions and spit
#out the lagrangian for given constraints, lambda and rho
(f_k0, fgrad_k0, fhes_k0) = f(x_k0);
(h_k0, hgrad_k0, hhes_k0) = h(x_k0);
L_k0 = f_k0 - lambda*h_k0 + (rho/2)*(norm(h_k0))^2;
Lgrd_k0 = fgrad_k0 + (-lambda + rho*h_k0)*hgrad_k0;
Lhes_k0 = eye(2);
return (L_k0, Lgrd_k0, Lhes_k0);
end
Out[2]:
In [8]:
function newtmin_BFGS_HW5(f, h, x_k0, lambda, rho, maxIts, optTol)
# obj is a function that evaluates the objective value,
# gradient and hessian a x.
# x0 is the starting point.
# maxIts = maximum interations
# optTol is the optimality tolerance on the gradient
its = 1;
Tol = 0;
err = 0;
fault_count = 0;
delta = 1;
# calculate the gradient, hessian and objective at x
#(f_k0, g_k0, h_k0) = f(x_k0);
(f_k0, g_k0, h_k0) = lagrangian(f, h, x_k0, lambda, rho);
tolFactor = optTol*norm(g_k0);
B_k0 = h_k0;
# Start the Loop
for i=1:maxIts
Tol = norm(g_k0);
if Tol<tolFactor
break;
end
step_k0 = B_k0\g_k0;
alpha = 1;
#(f_k1, g_k1, h_k1) = f(x_k0-alpha*step_k0);
(f_k1, g_k1, h_k1) = lagrangian(f, h, x_k0-alpha*step_k0, lambda, rho);
#Linesearch via Armijo Backtracking
while( f_k1 > f_k0 )
alpha=.5*alpha;
#(f_k1, g_k1, h_k1) = f(x_k0-alpha*step_k0);
(f_k1, g_k1, h_k1) = lagrangian(f, h, x_k0-alpha*step_k0, lambda, rho);
if alpha < ( 1e-2)
break;
end
end
x_k1 = x_k0 - alpha*(step_k0);
#(f_k1, g_k1, h_k1) = f(x_k1);
(f_k1, g_k1, h_k1) = lagrangian(f, h, x_k1, lambda, rho);
y = g_k1 - g_k0;
s = x_k1 - x_k0;
#BFGS Update
dB = (y*y')/dot(s,y) - (B_k0*(s*s')*B_k0)/dot(s,B_k0*s);
B_k0 = B_k0 + dB;
x_k0 = x_k1;
#(f_k0, g_k0, h_k0) = f(x_k0);
(f_k0, g_k0, h_k0) = lagrangian(f, h, x_k0, lambda, rho);
its+=1;
end
x = x_k0;
f = f_k0;
B = B_k0;
#return(f,its,Tol,B,x);
return (x);
end
Out[8]:
In [60]:
function augLagrange(f, h, x_k0, maxIts, optTol)
#implementation of the augmented lagrangian method
lambda_k0 = -.55;
rho_k0 = .001;
beta = 1.5;
gamma = .05;
its = 1;
#initial call to objective and constraint functions
(f_k0, fgrad_k0, fhes_k0) = f(x_k0);
(h_k0, hgrad_k0, hhes_k0) = h(x_k0);
tol = norm(fgrad_k0 - lambda_k0*hgrad_k0)*optTol;
for i=1:maxIts
#check the tolertance condition
Tol_H = norm(h_k0);
Tol_L = norm(fgrad_k0 - lambda_k0*hgrad_k0);
if (Tol_L < tol) && (Tol_H < tol)
break;
end
#find the minimizer for the lambda and rho pair
x_k1 = newtmin_BFGS_HW5(f, h, x_k0, lambda_k0, rho_k0, 200, 1e-6);
(f_k1, fgrad_k1, fhes_k1) = f(x_k1);
(h_k1, hgrad_k1, hhes_k1) = h(x_k1);
if norm(h_k1) <= gamma*norm(h_k0);
rho_k1 = rho_k0;
lambda_k1 = lambda_k0 + rho_k0*h_k1;
else rho_k1 = beta*rho_k0;
lambda_k1 = lambda_k0;
end
#lambda_k1 = lambda_k0 + rho_k1*h_k1;
lambda_k0 = lambda_k1;
rho_k0 = rho_k1;
x_k0 = x_k1;
f_k = f_k1;
h_k = h_k1;
its+=1;
end
lambda = lambda_k0;
rho = rho_k0;
x = x_k0;
f_k = f_k0;
h_k = h_k0;
Tol = norm(fgrad_k0 - lambda_k0*hgrad_k0);
return (x,rho,lambda,f_k,h_k,its,Tol);
end
Out[60]:
In [59]:
(L_k0, Lgrd_k0, Lhes_k0) = lagrangian(x -> (F_1(x), F_1_grd(x), F_1_hes(x)),x -> (H_1(x), H_1_grd(x), H_1_hes(x)), [-1.9,-1.9], 3731, 6471)
Out[59]:
In [68]:
(x,rho,lambda,f_k,h_k,its,Tol)=augLagrange(x -> (F_6(x), F_6_grd(x), F_6_hes(x)),x -> (H_6(x), H_6_grd(x), H_6_hes(x)), F_6_x0, 200, .01)
Out[68]:
In [ ]: