## ConstrainedOptim.jl

Optim.jl style constrained optimization (experimental)
Author JuliaNLSolvers
Popularity
12 Stars
Updated Last
1 Year Ago
Started In
February 2018

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)^2 + 100.0 * (x - x^2)^2

g = -2.0 * (1.0 - x) - 400.0 * (x - x^2) * x
g = 200.0 * (x - x^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x + 1200.0 * x^2
h[1, 2] = -400.0 * x
h[2, 1] = -400.0 * x
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
=============

μ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
``````

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
``````

### 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^2 + x^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 `j`th 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 = x^2 + x^2; c)
function con_jacobian!(J, x)
J[1,1] = 2*x
J[1,2] = 2*x
J
end
function con_h!(h, x, λ)
h[1,1] += λ*2
h[2,2] += λ*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
``````

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*sin(x)-x`.

```function con_c!(c, x)
c = x^2 + x^2     # First constraint
c = x*sin(x)-x # Second constraint
c
end
function con_jacobian!(J, x)
# First constraint
J[1,1] = 2*x
J[1,2] = 2*x
# Second constraint
J[2,1] = x*cos(x)-1.0
J[2,2] = sin(x)
J
end
function con_h!(h, x, λ)
# First constraint
h[1,1] += λ*2
h[2,2] += λ*2
# Second constraint
h[1,1] += λ*x*-sin(x)
h[1,2] += λ*cos(x)
# 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())```