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'|               = 8.88e-14 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.78e-13 ≰ 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:    41
    f(x) calls:    63
    ∇f(x) calls:   63

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:    28
    f(x) calls:    109
    ∇f(x) calls:   109

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/tAcSy/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/tAcSy/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/tAcSy/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")){Expr})() at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Expanders.jl:480
 [6] cd(::getfield(Documenter.Expanders, Symbol("##8#10")){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,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,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] (::getfield(Documenter, Symbol("##2#3")))() at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Documenter.jl:204
 [18] cd(::getfield(Documenter, Symbol("##2#3")), ::String) at ./file.jl:96
 [19] #makedocs#1(::Bool, ::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:doctest,),Tuple{Bool}}}, ::Function) at /home/travis/.julia/packages/Documenter/Qo3Yk/src/Documenter.jl:203
 [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:    34
    f(x) calls:    158
    ∇f(x) calls:   158

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:    215
    ∇f(x) calls:   215

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.