Nonlinear constrained optimization
Tip
This example is also available as a Jupyter notebook: ipnewton_basics.ipynb
The nonlinear constrained optimization interface in Optim assumes that the user can write the optimization problem in the following way.
For equality constraints on $x_j$ or $c(x)_j$ you set those particular entries of bounds to be equal, $l_j=u_j$. Likewise, setting $l_j=-\infty$ or $u_j=\infty$ means that the constraint is unbounded from below or above respectively.
Constrained optimization with IPNewton
We will go through examples on how to use the constraints interface with the interior-point Newton optimization algorithm IPNewton.
Throughout these examples we work with the standard Rosenbrock function. The objective and its derivatives are given by
fun(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end
function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;
Optimization interface
To solve a constrained optimization problem we call the optimize method
optimize(d::AbstractObjective, constraints::AbstractConstraints, initial_x::Tx, method::ConstrainedOptimizer, options::Options)
We can create instances of AbstractObjective and AbstractConstraints using the types TwiceDifferentiable and TwiceDifferentiableConstraints from the package NLSolversBase.jl.
Box minimzation
We want to optimize the Rosenbrock function in the box $-0.5 \leq x \leq 0.5$, starting from the point $x_0=(0,0)$. Box constraints are defined using, for example, TwiceDifferentiableConstraints(lx, ux).
x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)
res = optimize(df, dfc, x0, IPNewton())
* Status: success
* Candidate solution
Minimizer: [5.00e-01, 2.50e-01]
Minimum: 2.500000e-01
* Found with
Algorithm: Interior Point Newton
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 4.39e-10 ≰ 0.0e+00
|x - x'|/|x'| = 8.79e-10 ≰ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 1.00e+00 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 43
f(x) calls: 68
∇f(x) calls: 68
If we only want to set lower bounds, use ux = fill(Inf, 2)
ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 7.987239e-20
* Found with
Algorithm: Interior Point Newton
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 3.54e-10 ≰ 0.0e+00
|x - x'|/|x'| = 3.54e-10 ≰ 0.0e+00
|f(x) - f(x')| = 2.40e-19 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00
|g(x)| = 8.83e-09 ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 35
f(x) calls: 63
∇f(x) calls: 63
Defining "unconstrained" problems
An unconstrained problem can be defined either by passing Inf bounds or empty arrays. Note that we must pass the correct type information to the empty lx and ux
lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 5.998937e-19
* Found with
Algorithm: Interior Point Newton
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 1.50e-09 ≰ 0.0e+00
|x - x'|/|x'| = 1.50e-09 ≰ 0.0e+00
|f(x) - f(x')| = 1.80e-18 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00
|g(x)| = 7.92e-09 ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 34
f(x) calls: 63
∇f(x) calls: 63
Generic nonlinear constraints
We now consider the Rosenbrock problem with a constraint on
We pass the information about the constraints to optimize by defining a vector function c(x) and its Jacobian J(x).
The Hessian information is treated differently, by considering the Lagrangian of the corresponding slack-variable transformed optimization problem. This is similar to how the CUTEst library works. Let $H_j(x)$ represent the Hessian of the $j$th component $c(x)_j$ of the generic constraints. and $\lambda_j$ the corresponding dual variable in the Lagrangian. Then we want the constraint object to add the values of $H_j(x)$ to the Hessian of the objective, weighted by $\lambda_j$.
The Julian form for the supplied function $c(x)$ and the derivative information is then added in the following way.
con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
J
end
function con_h!(h, x, λ)
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
end;
Note that con_h! adds the λ-weighted Hessian value of each element of c(x) to the Hessian of fun.
We can then optimize the Rosenbrock function inside the ball of radius $0.5$.
lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
* Status: success
* Candidate solution
Minimizer: [4.56e-01, 2.06e-01]
Minimum: 2.966216e-01
* Found with
Algorithm: Interior Point Newton
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 7.71e-01 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 32
f(x) calls: 111
∇f(x) calls: 111
We can add a lower bound on the constraint, and thus optimize the objective on the annulus with inner and outer radii $0.1$ and $0.5$ respectively.
lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
┌ Warning: Initial guess is not an interior point
└ @ Optim ~/.julia/packages/Optim/SFpsz/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:111
Stacktrace:
[1] initial_state(::IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}, ::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{typeof(Main.ex-ipnewton_basics.con_c!),typeof(Main.ex-ipnewton_basics.con_jacobian!),typeof(Main.ex-ipnewton_basics.con_h!),Float64}, ::Array{Float64,1}) at /home/travis/.julia/packages/Optim/SFpsz/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:112
[2] optimize(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::TwiceDifferentiableConstraints{typeof(Main.ex-ipnewton_basics.con_c!),typeof(Main.ex-ipnewton_basics.con_jacobian!),typeof(Main.ex-ipnewton_basics.con_h!),Float64}, ::Array{Float64,1}, ::IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}) at /home/travis/.julia/packages/Optim/SFpsz/src/multivariate/solvers/constrained/ipnewton/interior.jl:196 (repeats 2 times)
[3] top-level scope at none:0
[4] eval at ./boot.jl:319 [inlined]
[5] (::getfield(Documenter.Expanders, Symbol("##8#10")){Module,Expr})() at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Expanders.jl:480
[6] cd(::getfield(Documenter.Expanders, Symbol("##8#10")){Module,Expr}, ::String) at ./file.jl:96
[7] #7 at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Expanders.jl:479 [inlined]
[8] (::getfield(Documenter.Utilities, Symbol("##18#19")){getfield(Documenter.Expanders, Symbol("##7#9")){Documenter.Documents.Page,Module,Expr},Base.PipeEndpoint,Base.PipeEndpoint,Pipe,Array{UInt8,1}})() at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Utilities/Utilities.jl:591
[9] with_logstate(::getfield(Documenter.Utilities, Symbol("##18#19")){getfield(Documenter.Expanders, Symbol("##7#9")){Documenter.Documents.Page,Module,Expr},Base.PipeEndpoint,Base.PipeEndpoint,Pipe,Array{UInt8,1}}, ::Base.CoreLogging.LogState) at ./logging.jl:395
[10] with_logger(::Function, ::Logging.ConsoleLogger) at ./logging.jl:491
[11] withoutput at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Utilities/Utilities.jl:589 [inlined]
[12] runner(::Type{Documenter.Expanders.ExampleBlocks}, ::Markdown.Code, ::Documenter.Documents.Page, ::Documenter.Documents.Document) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Expanders.jl:478
[13] dispatch(::Type{Documenter.Expanders.ExpanderPipeline}, ::Markdown.Code, ::Vararg{Any,N} where N) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Selectors.jl:168
[14] expand(::Documenter.Documents.Document) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Expanders.jl:31
[15] runner(::Type{Documenter.Builder.ExpandTemplates}, ::Documenter.Documents.Document) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Builder.jl:178
[16] dispatch(::Type{Documenter.Builder.DocumentPipeline}, ::Documenter.Documents.Document) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Selectors.jl:168
[17] #2 at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Documenter.jl:204 [inlined]
[18] cd(::getfield(Documenter, Symbol("##2#3")){Documenter.Documents.Document}, ::String) at ./file.jl:96
[19] #makedocs#1 at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Documenter.jl:203 [inlined]
[20] (::getfield(Documenter, Symbol("#kw##makedocs")))(::NamedTuple{(:doctest,),Tuple{Bool}}, ::typeof(makedocs)) at ./none:0
[21] top-level scope at none:0
[22] include at ./boot.jl:317 [inlined]
[23] include_relative(::Module, ::String) at ./loading.jl:1044
[24] include(::Module, ::String) at ./sysimg.jl:29
[25] exec_options(::Base.JLOptions) at ./client.jl:266
[26] _start() at ./client.jl:425
* Status: success
* Candidate solution
Minimizer: [4.56e-01, 2.06e-01]
Minimum: 2.966216e-01
* Found with
Algorithm: Interior Point Newton
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 7.71e-01 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 36
f(x) calls: 162
∇f(x) calls: 162
Note that the algorithm warns that the Initial guess is not an interior point. IPNewton can often handle this, however, if the initial guess is such that c(x) = u_c, then the algorithm currently fails. We may fix this in the future.
Multiple constraints
The following example illustrates how to add an additional constraint. In particular, we add a constraint function
function con2_c!(c, x)
c[1] = x[1]^2 + x[2]^2 ## First constraint
c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
c
end
function con2_jacobian!(J, x)
# First constraint
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
# Second constraint
J[2,1] = x[2]*cos(x[1])-1.0
J[2,2] = sin(x[1])
J
end
function con2_h!(h, x, λ)
# First constraint
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
# Second constraint
h[1,1] += λ[2]*x[2]*-sin(x[1])
h[1,2] += λ[2]*cos(x[1])
# Symmetrize h
h[2,1] = h[1,2]
h
end;
We generate the constraint objects and call IPNewton with initial guess $x_0 = (0.25,0.25)$.
x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
* Status: success
* Candidate solution
Minimizer: [-1.60e-19, -1.95e-18]
Minimum: 1.000000e+00
* Found with
Algorithm: Interior Point Newton
Initial Point: [2.50e-01, 2.50e-01]
* Convergence measures
|x - x'| = 6.90e-10 ≰ 0.0e+00
|x - x'|/|x'| = 3.55e+08 ≰ 0.0e+00
|f(x) - f(x')| = 1.38e-09 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 1.38e-09 ≰ 0.0e+00
|g(x)| = 2.00e+00 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 29
f(x) calls: 219
∇f(x) calls: 219
Plain Program
Below follows a version of the program without any comments. The file is also available here: ipnewton_basics.jl
using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide
fun(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end
function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;
x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)
res = optimize(df, dfc, x0, IPNewton())
ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
J
end
function con_h!(h, x, λ)
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
end;
lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
function con2_c!(c, x)
c[1] = x[1]^2 + x[2]^2 ## First constraint
c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
c
end
function con2_jacobian!(J, x)
# First constraint
J[1,1] = 2*x[1]
J[1,2] = 2*x[2]
# Second constraint
J[2,1] = x[2]*cos(x[1])-1.0
J[2,2] = sin(x[1])
J
end
function con2_h!(h, x, λ)
# First constraint
h[1,1] += λ[1]*2
h[2,2] += λ[1]*2
# Second constraint
h[1,1] += λ[2]*x[2]*-sin(x[1])
h[1,2] += λ[2]*cos(x[1])
# Symmetrize h
h[2,1] = h[1,2]
h
end;
x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl
This page was generated using Literate.jl.