Notebook

Nonlinear constrained optimization¶

The nonlinear constrained optimization interface in Optim assumes that the user can write the optimization problem in the following way.

$\min_{x\in\mathbb{R}^n} f(x) \quad \text{such that}\\ l_x \leq \phantom{c(}x\phantom{)} \leq u_x \\ l_c \leq c(x) \leq u_c.$

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.

In [1]:

using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide

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

In [2]:

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 minimization¶

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).

In [3]:

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())

Out[3]:

 * Status: success

 * Candidate solution
    Final objective value:     2.500000e-01

 * Found with
    Algorithm:     Interior Point Newton

 * 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:   1  (vs limit Inf)
    Iterations:    43
    f(x) calls:    68
    ∇f(x) calls:   68

Like the rest of Optim, you can also use autodiff=:forward and just pass in fun.

If we only want to set lower bounds, use ux = fill(Inf, 2)

In [4]:

ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

Out[4]:

 * Status: success (objective increased between iterations)

 * Candidate solution
    Final objective value:     7.987239e-20

 * Found with
    Algorithm:     Interior Point Newton

 * 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

In [5]:

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())

Out[5]:

 * Status: success

 * Candidate solution
    Final objective value:     5.998937e-19

 * Found with
    Algorithm:     Interior Point Newton

 * 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

$c(x)_1 = x_1^2 + x_2^2.$

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.

In [6]:

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$ .

In [7]:

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())

Out[7]:

 * Status: success

 * Candidate solution
    Final objective value:     2.966216e-01

 * Found with
    Algorithm:     Interior Point Newton

 * 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.

In [8]:

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 ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:136

Stacktrace:
  [1] initial_state(method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol}, options::Optim.Options{Float64, Nothing}, d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var"##231".con_c!), typeof(Main.var"##231".con_jacobian!), typeof(Main.var"##231".con_h!), Float64}, initial_x::Vector{Float64})
    @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:137
  [2] optimize
    @ ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:289 [inlined]
  [3] optimize(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var"##231".con_c!), typeof(Main.var"##231".con_jacobian!), typeof(Main.var"##231".con_h!), Float64}, initial_x::Vector{Float64}, method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol})
    @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:289
  [4] top-level scope
    @ ~/work/Optim.jl/Optim.jl/docs/src/examples/generated/ipnewton_basics.ipynb:3
  [5] eval
    @ ./boot.jl:430 [inlined]
  [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2734
  [7] #47
    @ ~/.julia/packages/Literate/TujCy/src/Literate.jl:932 [inlined]
  [8] task_local_storage(body::Literate.var"#47#49"{String, Bool, Module, String}, key::Symbol, val::String)
    @ Base ./task.jl:315
  [9] #46
    @ ~/.julia/packages/Literate/TujCy/src/Literate.jl:930 [inlined]
 [10] (::IOCapture.var"#5#9"{Core.TypeofBottom, Literate.var"#46#48"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.PipeEndpoint, Base.PipeEndpoint})()
    @ IOCapture ~/.julia/packages/IOCapture/Y5rEA/src/IOCapture.jl:170
 [11] with_logstate(f::IOCapture.var"#5#9"{Core.TypeofBottom, Literate.var"#46#48"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.PipeEndpoint, Base.PipeEndpoint}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging ./logging/logging.jl:522
 [12] with_logger(f::Function, logger::Base.CoreLogging.ConsoleLogger)
    @ Base.CoreLogging ./logging/logging.jl:632
 [13] capture(f::Literate.var"#46#48"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Any})
    @ IOCapture ~/.julia/packages/IOCapture/Y5rEA/src/IOCapture.jl:167
 [14] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool, continue_on_error::Bool)
    @ Literate ~/.julia/packages/Literate/TujCy/src/Literate.jl:928
 [15] execute_block
    @ ~/.julia/packages/Literate/TujCy/src/Literate.jl:912 [inlined]
 [16] execute_notebook(nb::Dict{Any, Any}; inputfile::String, fake_source::String, softscope::Bool, continue_on_error::Bool)
    @ Literate ~/.julia/packages/Literate/TujCy/src/Literate.jl:828
 [17] (::Literate.var"#40#42"{Dict{String, Any}})()
    @ Literate ~/.julia/packages/Literate/TujCy/src/Literate.jl:804
 [18] cd(f::Literate.var"#40#42"{Dict{String, Any}}, dir::String)
    @ Base.Filesystem ./file.jl:112
 [19] jupyter_notebook(chunks::Vector{Literate.Chunk}, config::Dict{String, Any})
    @ Literate ~/.julia/packages/Literate/TujCy/src/Literate.jl:803
 [20] notebook(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{execute::Bool})
    @ Literate ~/.julia/packages/Literate/TujCy/src/Literate.jl:740
 [21] top-level scope
    @ ~/work/Optim.jl/Optim.jl/docs/generate.jl:22
 [22] include(fname::String)
    @ Main ./sysimg.jl:38
 [23] top-level scope
    @ ~/work/Optim.jl/Optim.jl/docs/make.jl:14
 [24] include(mod::Module, _path::String)
    @ Base ./Base.jl:557
 [25] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:323
 [26] _start()
    @ Base ./client.jl:531

Out[8]:

 * Status: success

 * Candidate solution
    Final objective value:     2.966216e-01

 * Found with
    Algorithm:     Interior Point Newton

 * 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

$c(x)_2 = x_2\sin(x_1)-x_1$

In [9]:

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)$ .

In [10]:

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())

Out[10]:

 * Status: success

 * Candidate solution
    Final objective value:     1.000000e+00

 * Found with
    Algorithm:     Interior Point Newton

 * 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

This notebook was generated using Literate.jl.