``````

In [1]:

function newtmin( obj, x0; maxIts=100, optTol=1e-6)
# Newton method with Armijo backtracking and some treatment with Hessian, treatment 2
# Minimize a function f using Newton’s method.
# obj:  a function that evaluates the objective value,
# gradient, and Hessian at a point x, i.e.,
# (f, g, H) = obj(x)
# x0: starting point.
# maxIts (optional): maximum number of iterations.
# optTol (optional): optimality tolerance based on
mu=1e-4
epsilon=0.01
x=x0
status = 0
its = 0
(f0,g0,H0)=obj(x0)
(f,g,H)=obj(x)
(V,S)=eig(H)
bar=ones(size(H,1))*epsilon
bar=Diagonal(bar)
bH=S*max(bar,abs(Diagonal(V)))*S'
while status != 1
alpha = 1
xnew = x-alpha*inv(bH)*g
(fnew,gnew,Hnew)=obj(xnew)
sts =-fnew+f-alpha*mu*g'*inv(bH)*g
while sts[1]<0
alpha=alpha/2
xnew = x-alpha*inv(bH)*g
(fnew,gnew,Hnew)=obj(xnew)
sts=-fnew+f-alpha*mu*g'*inv(bH)*g
end
x = x-alpha*inv(bH)*g
(f,g,H)=obj(x)
(V,S)=eig(H)
bar=ones(size(H,1))*epsilon
bar=Diagonal(bar)
bH=S*max(bar,abs(Diagonal(V)))*S'
its = its+1
if norm(g)<= optTol*norm(g0)
status = 1
end
if its>maxIts
status = 1
end
end
return (x, its)
end

``````
``````

Out[1]:

newtmin (generic function with 1 method)

``````
``````

In [2]:

using ForwardDiff
using ReverseDiffSource

``````

# Test problem 6

``````

In [16]:

x0=[-1.2;1]
f(x)=(1-x[1])^2
c(x)=10*(x[2]-x[1]^2)
H(x)=ForwardDiff.hessian(x -> f(x))(x)

``````
``````

Out[16]:

J (generic function with 1 method)

``````
``````

In [17]:

function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=ForwardDiff.hessian(z -> d1(z))(z)
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=ForwardDiff.hessian(z -> d2(z))(z)
result=H(x)+dh1(x)+dh2(x)
return result
end

``````
``````

Out[17]:

nxxLrou (generic function with 1 method)

``````
``````

In [22]:

function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou)   # objective value at x
lagg = nxLrou(x,y,rou)  # gradient at x
lagH = nxxLrou(x,y,rou)   # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
y=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
sts=1
end
x=xnew
end
return x
end

``````
``````

Out[22]:

auglag (generic function with 1 method)

``````
``````

In [27]:

auglag( x0; maxIts=100, optTol=1e-6)

``````
``````

Out[27]:

2-element Array{Float64,1}:
1.0
1.0

``````

# Test problem 7

``````

In [53]:

x0=[2;2]
f(x)=log(1+x[1]^2)-x[2]
c(x)=(1+x[1]^2)^2+x[2]^2-4
g(x)=[2*x[1]/(1+x[1]^2);-1]
H(x)=[(2-2*x[1]^2)/(1+x[1]^2)^2 0;0 0]

``````
``````

Out[53]:

H (generic function with 1 method)

``````
``````

In [54]:

function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(x)=-y*[4+12*x[1]^2 0;0 2]
d2(z)=rou/2*c(z)'*c(z)
dh2(x)=rou*[4+12*x[1]^2 0;0 2]*c(x)+rou*J(x)'*J(x)
result=H(x)+dh1(x)+dh2(x)
return result
end

``````
``````

Out[54]:

nxxLrou (generic function with 1 method)

``````
``````

In [55]:

function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou)   # objective value at x
lagg = nxLrou(x,y,rou)  # gradient at x
lagH = nxxLrou(x,y,rou)   # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
ynew=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
if abs(abs(J(xnew)*g(xnew)/(norm(J(xnew))*norm(g(xnew))))-1)[1]<optTol
sts=1
end
end
x=xnew
y=ynew
end
return x
end
auglag( x0; maxIts=100, optTol=1e-6)

``````
``````

Out[55]:

2-element Array{Float64,1}:
-9.51374e-42
1.73206

``````

# Test problem 8

``````

In [6]:

x0=[2;1]
f(x)=-1
c(x)=[x[1]^2+x[2]^2-25; x[1]*x[2]-9]
g(x)=[0;0]
H(x)=[0 0 ;0 0 ]
J(x)=([2*x[1] x[2];2*x[2] x[1]])'
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=-y[1]*[2 0;0 2]-y[2]*[0 1; 1 0]
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=rou*[2 0;0 2]*c(x)[1]+rou*[0 1; 1 0]*c(x)[2]+rou*J(x)'*J(x)
result=H(x)+dh1(x)+dh2(x)
return result
end

``````
``````

Out[6]:

nxxLrou (generic function with 1 method)

``````
``````

In [54]:

function auglag( x0; maxIts=100, optTol=1e-6)
y0=[1;1]
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou)   # objective value at x
lagg = nxLrou(x,y,rou)  # gradient at x
lagH = nxxLrou(x,y,rou)   # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
ynew=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
sts=1
end
x=xnew
y=ynew
end
return x
end
auglag( x0; maxIts=100, optTol=1e-6)

``````
``````

Out[54]:

2-element Array{Float64,1}:
4.6016
1.95585

``````

# Test problem 9

``````

In [119]:

x0=[0;0]
f(x)=sin(pi*x[1]/12)*cos(pi*x[2]/16)
c(x)=4*x[1]-3*x[2]
g(x)=[pi/12*cos(pi*x[1]/12)*cos(pi*x[2]/16); -pi/16*sin(pi*x[1]/12)*sin(pi*x[2]/16)]
H(x)=[-pi^2/(12^2)*sin(pi*x[1]/12)*cos(pi*x[2]/16) -pi^2/(12*16)*cos(pi*x[1]/12)*sin(pi*x[2]/16);-pi^2/(16*12)*cos(pi*x[1]/12)*sin(pi*x[2]/16) -pi^2/(16^2)*sin(pi*x[1]/12)*cos(pi*x[2]/16)]
function Lrou(x,y,rou)
result=f(x)-y'*c(x)+rou/2*c(x)'*c(x)
return result
end
function nxLrou(x,y,rou)
result=g(x)-J(x)'*y+rou*J(x)'*c(x)
result=result[:,1]
#result=convert(Array{Float64,1}, result)
return result
end
function nxxLrou(x,y,rou)
d1(z)=-y'*c(z)
dh1(z)=ForwardDiff.hessian(z -> d1(z))(z)
d2(z)=rou/2*c(z)'*c(z)
dh2(z)=ForwardDiff.hessian(z -> d2(z))(z)
result=H(x)+dh1(x)+dh2(x)
return result
end

``````
``````

Out[119]:

nxxLrou (generic function with 1 method)

``````
``````

In [120]:

function auglag( x0; maxIts=100, optTol=1e-6)
y0=1
c0=max(norm(c(x0)),1)
d0=norm(g(x0)-J(x0)'*y0)
sts=0
y=y0
x=x0
rou=1
while sts==0
function obj(x)
lagf = Lrou(x,y,rou)   # objective value at x
lagg = nxLrou(x,y,rou)  # gradient at x
lagH = nxxLrou(x,y,rou)   # Hessian at x
return (lagf,lagg,lagH)
end
(xnew,its)=newtmin( obj, x; maxIts=10, optTol=1e-6)
if norm(c(xnew))<1/2*norm(c(x))
y=y+rou*c(xnew)
else
rou=rou*10
end
if norm(c(xnew))<optTol*c0
if abs(abs(J(xnew)*g(xnew)/(norm(J(xnew))*norm(g(xnew))))-1)[1]<optTol
sts=1
end
end
x=xnew
end
return x
end

``````
``````

Out[120]:

auglag (generic function with 1 method)

``````
``````

In [121]:

auglag( x0; maxIts=100, optTol=1e-6)

``````
``````

Out[121]:

2-element Array{Float64,1}:
-3.0
-4.0

``````