Note: This functionality is now merged into Optim.jl

ConstrainedOptim

This package adds support for constrained optimization algorithms to the package Optim. We intend to merge the code in ConstrainedOptim with Optim when the interfaces and algorithms in this repository have been tested properly. Feedback is very much appreciated, either via gitter or by creating an issue or PR on github.

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

Multiple nonlinear constraints can be set by considering c(x) as a vector. An equality constraint g(x) = 0 is then defined by setting, say, c(x)_1 = g(x), l_{c,1}= u_{c,1} = 0.

Example usage

We will go through examples of how to use ConstrainedOptim and illustrate how to use the constraints interface with an interior-point Newton optimization algorithm. 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.

This package contains one ConstrainedOptimizer method called IPNewton. To get information on the keywords used to construct method instances, use the Julia REPL help prompt (?).

help?> IPNewton
search: IPNewton

Interior-point Newton
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡

Constructor
=============

IPNewton(; linesearch::Function = ConstrainedOptim.backtrack_constrained_grad,
μ0::Union{Symbol,Number} = :auto,
show_linesearch::Bool = false)

The initial barrier penalty coefficient μ0 can be chosen as a number, or set to :auto to let the
algorithm decide its value, see initialize_μ_λ!.

*Note*: For constrained optimization problems, we recommend
always enabling `allow_f_increases` and `successive_f_tol` in the options passed to `optimize`.
The default is set to `Optim.Options(allow_f_increases = true, successive_f_tol = 2)`.

As of February 2018, the line search algorithm is specialised for constrained interior-point
methods. In future we hope to support more algorithms from LineSearches.jl.

Description
=============

The IPNewton method implements an interior-point primal-dual Newton algorithm for solving nonlinear,
constrained optimization problems. See Nocedal and Wright (Ch. 19, 2006) for a discussion of
interior-point methods for constrained optimization.

References
============

The algorithm was originally written by Tim Holy (@timholy, tim.holy@gmail.com).

•    J Nocedal, SJ Wright (2006), Numerical optimization, second edition. Springer.

•    A Wächter, LT Biegler (2006), On the implementation of an interior-point filter
line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106
(1), 25-57.

Box minimzation

We want to optimize the Rosenbrock function in the box -0.5 ≤ x ≤ 0.5, starting from the point x0=zeros(2). 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 = [1.0, 1.0]
dfc = TwiceDifferentiableConstraints(lx, ux)

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

The output from res is

Results of Optimization Algorithm
 * Algorithm: Interior Point Newton
 * Starting Point: [0.0,0.0]
 * Minimizer: [0.5,0.2500000000000883]
 * Minimum: 2.500000e-01
 * Iterations: 41
 * Convergence: true
   * |x - x'| ≤ 1.0e-32: false
     |x - x'| = 8.88e-14
   * |f(x) - f(x')| ≤ 1.0e-32 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = 1.00e+00
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 63
 * Gradient 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())

The output from res is now

Results of Optimization Algorithm
 * Algorithm: Interior Point Newton
 * Starting Point: [0.0,0.0]
 * Minimizer: [0.9999999998342594,0.9999999996456271]
 * Minimum: 7.987239e-20
 * Iterations: 35
 * Convergence: true
   * |x - x'| ≤ 1.0e-32: false
     |x - x'| = 3.54e-10
   * |f(x) - f(x')| ≤ 1.0e-32 |f(x)|: false
     |f(x) - f(x')| = 3.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: true
     |g(x)| = 8.83e-09
   * Stopped by an increasing objective: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 63
 * Gradient 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())

Generic nonlinear constraints

We now consider the Rosenbrock problem with a constraint on x[1]^2 + x[2]^2.

We pass the information about the constraints to the 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 jth component of c(x), and λ_j the corresponding dual variable in the Lagrangian. Then we want the constraint object to add the values of the 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())

The output from the optimization is

Results of Optimization Algorithm
 * Algorithm: Interior Point Newton
 * Starting Point: [0.0,0.0]
 * Minimizer: [0.45564896414551875,0.2058737998704899]
 * Minimum: 2.966216e-01
 * Iterations: 28
 * Convergence: true
   * |x - x'| ≤ 1.0e-32: true
     |x - x'| = 0.00e+00
   * |f(x) - f(x')| ≤ 1.0e-32 |f(x)|: false
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = 7.71e-01
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 109
 * Gradient 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())

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 hope to fix this in the future.

Multiple constraints

The following example illustrates how to add an additional constraint. In particular, we add a constraint c_2(x) = x[2]*sin(x[1])-x[1].

function con_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 con_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 con_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 x0 = [0.25,0.25].

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())