Apparently, I still cannot get inequality constraints working, but now using Optimization.jl and IPOPT via the OptimizationMOI package.
Problem
The problem is that IPOPT violates my inequality constraint during a step, which then causes an exception in the objective function. The objective function calculates log(variance), and the constraint is there to prevent variance from being negative.
- IPOPT violates this constraint:
-
a - 1/dtmust be negative; - During optimization I get
a - 1/dt == 32.36, which is not negative;
-
-
variancebecomes negative:-2.4941978436429695 < 0, - the
logfunction raises aDomainError.
I’m specifically using IPOPT because it’s an interior point optimizer, so it should never violate constraints, AFAIK. Yet it does.
Code
import Pkg
Pkg.activate(temp=true)
Pkg.add([
Pkg.PackageSpec(name="Optimization", version="3.10.0"),
Pkg.PackageSpec(name="OptimizationMOI", version="0.1.5"),
Pkg.PackageSpec(name="Ipopt", version="1.1.0"),
])
import Optimization
import OptimizationMOI, Ipopt
const AV{T} = AbstractVector{T}
function model_constraints!(out::AV{<:Real}, u::AV{<:Real}, data)
# Model parameters
dt, a, b = u
out[1] = a - 1/dt # Must be < 0
@info "Must be NEGATIVE: $(out[1])"
end
function model_variance(u::AV{T}, data::AV{<:Real}) where T<:Real
# Model parameters
dt, a, b = u
# Compute variance
variance = zeros(T, length(data))
variance[1] = one(T)
for t in 1:(length(data) - 1)
variance[t+1] = (1 - dt * a) * variance[t] + dt * data[t]^2 + dt * b
end
variance
end
function model_loss(u::AV{T}, data::AV{<:Real})::T where T<:Real
variance = model_variance(u, data)
loglik::T = zero(T)
for (r, var) in zip(data, variance)
loglik += -(log(2Ď€) + log(var) + r^2 / var) / 2
end
-loglik / length(data)
end
function model_fit(u0::AV{T}, data::AV{<:Real}) where T<:Real
func = Optimization.OptimizationFunction(
model_loss, Optimization.AutoForwardDiff(),
cons=model_constraints!
)
prob = Optimization.OptimizationProblem(
func, u0, data,
# 0 < dt < 1 && 1 < a < Inf && 0 < b < Inf
lb=T[0.0, 1.0, 0.0], ub=T[1.0, Inf, Inf],
# ^dt ^a ^b ^dt ^a ^b <= model parameters
lcons=T[-Inf], ucons=T[0.0] # a - 1/dt < 0
)
sol = Optimization.solve(prob, Ipopt.Optimizer())
sol.u
end
let
data = [
2.1217711584057386, -0.28350145551002465, 2.3593492969513004, 0.192856733601849, 0.4566485836385113, 1.332717934013979, -1.286716619379847, 0.9868669960185211, 2.2358674776395224, -2.7933975791568098,
1.2555871497124622, 1.276879759908467, -0.8392016987911409, -1.1580875182201849, 0.33201646080578456, -0.17212553408696898, 1.1275285626369556, 0.23041139849229036, 1.648423577528424, 2.384823597473343,
-0.4005518932539747, -1.117737311211693, -0.9490152960583265, -1.1454539355078672, 1.4158585811404159, -0.18926972177257692, -0.2867541528181491, -1.2077459688543788, -0.6397173049620141, 0.66147783407023,
0.049805188778543466, 0.902540117368457, -0.7018417933284938, 0.47342354473843684, 1.2620345361591596, -1.1483844812087018, -0.06487285080802752, 0.39020117013487715, -0.38454491504165356, 1.5125786171885645,
-0.6751768274451174, 0.490916740658628, 0.012872300530924086, 0.46532447715746716, 0.34734421531357157, 0.3830452463549559, -0.8730874028738718, 0.4333151627834603, -0.40396180775692375, 2.0794821773418497,
-0.5392735774960918, 0.6519326323752113, -1.4844713145398716, 0.3688828625691108, 1.010912990717231, 0.5018274939956874, 0.36656889279915833, -0.11403975693239479, -0.6460314660359935, -0.41997005020823147,
0.9652752515820495, -0.37375868692702047, -0.5780729659197872, 2.642742798278919, 0.5076984117208074, -0.4906395089461916, -1.804352047187329, -0.8596663844837792, -0.7510485548262176, -0.07922589350581195,
1.7201304839487317, 0.9024493222130577, -1.8216089665357902, 1.3929269238775426, -0.08410752079538407, 0.6423068180438288, 0.6615201016351212, 0.18546977816594887, -0.717521690742993, -1.0224309324751113,
1.7748350222721971, 0.1929546575877559, -0.1581871639724676, 0.20198379311238596, -0.6919373947349301, -0.9253274269423383, 0.549366272989534, -1.9302106783541606, 0.7197247279281573, -1.220334158468621,
-0.9187468058921053, -2.1452607604834184, -2.1558650694862687, -0.9387913392336701, -0.676637835687265, -0.16621998352492198, 0.5637177022958897, -0.5258315560278541, 0.8413359958184765, -0.9096866525337141
]
# u0 = [0 < dt < 1, 1 < a < 1/dt, 0 < b < Inf]
u0 = [0.3, 2.3333333333333335, 0.33333333333333337]
@assert 0 < u0[1] < 1
@assert 1 < u0[2] < 1 / u0[1]
@assert 0 < u0[3] < Inf
@info "Optimizing..." u0
model_fit(u0, data)
end
I’m setting up the optimization problem like this:
func = Optimization.OptimizationFunction(
model_loss, Optimization.AutoForwardDiff(),
cons=model_constraints!
)
prob = Optimization.OptimizationProblem(
func, u0, data,
# 0 < dt < 1 && 1 < a < Inf && 0 < b < Inf
lb=T[0.0, 1.0, 0.0], ub=T[1.0, Inf, Inf],
# ^dt ^a ^b ^dt ^a ^b <= model parameters
lcons=T[-Inf], ucons=T[0.0] # a - 1/dt < 0
)
I put model_constraints! into OptimizationFunction, as shown in the docs. lcons and ucons in OptimizationProblem are supposed to say that the output of model_constraints! must be between lcons==-Inf and ucons==0.0. In other words, it must be negative, less than ucons==0.0.
However, at some point it becomes positive.
Output
$ julia-1.8 my_model.jl
Activating new project at `/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_fD7r9m`
Updating registry at `~/.julia/registries/General.toml`
Resolving package versions...
Updating `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_fD7r9m/Project.toml`
[b6b21f68] + Ipopt v1.1.0
[7f7a1694] + Optimization v3.10.0
[fd9f6733] + OptimizationMOI v0.1.5
...
┌ Info: Optimizing...
│ u0 =
│ 3-element Vector{Float64}:
│ 0.3
│ 2.3333333333333335
â”” 0.33333333333333337
...
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 3
Number of nonzeros in Lagrangian Hessian.............: 6
[ Info: Must be NEGATIVE: Dual{ForwardDiff.Tag{Optimization.var"#65#82"{Int64}, Float64}}(-1.0,11.111111111111112,1.0,0.0)
[ Info: Must be NEGATIVE: -1.0
[ Info: Must be NEGATIVE: Dual{ForwardDiff.Tag{Optimization.var"#65#82"{Int64}, Float64}}(-1.0,11.111111111111112,1.0,0.0)
Total number of variables............................: 3
variables with only lower bounds: 2
variables with lower and upper bounds: 1
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1
...
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.8621407e+00 0.00e+00 2.10e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
...
9 1.5072838e+00 3.55e-01 3.44e-02 -1.7 5.79e+01 - 6.48e-01 1.00e+00H 1
[ Info: Must be NEGATIVE: -283.8381149547038
[ Info: Must be NEGATIVE: 32.36098192303031
ERROR: LoadError: DomainError with -2.4941978436429695:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math ./math.jl:33
[2] _log(x::Float64, base::Val{:â„Ż}, func::Symbol)
@ Base.Math ./special/log.jl:301
[3] log
@ ./special/log.jl:267 [inlined]
[4] model_loss(u::Vector{Float64}, data::Vector{Float64})
@ Main my_model.jl:37
[5] OptimizationFunction
@ ~/.julia/packages/SciMLBase/QqtZA/src/scimlfunctions.jl:3580 [inlined]
[6] eval_objective(moiproblem::OptimizationMOI.MOIOptimizationProblem{Float64, SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Optimization.var"#57#74"{ForwardDiff.GradientConfig{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, Float64, 3}}}, Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}}, Optimization.var"#60#77"{ForwardDiff.HessianConfig{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, Float64, 3}, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Float64}, Float64, 3}}}, Optimization.var"#56#73"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}}, Optimization.var"#63#80", Optimization.var"#64#81"{SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}}, Optimization.var"#66#83"{ForwardDiff.JacobianConfig{ForwardDiff.Tag{Optimization.var"#65#82"{Int64}, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#65#82"{Int64}, Float64}, Float64, 3}}}}, Optimization.var"#71#88"{Int64, Vector{ForwardDiff.HessianConfig{ForwardDiff.Tag{Optimization.var"#69#86"{Int64}, Float64}, Float64, 3, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#69#86"{Int64}, Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#69#86"{Int64}, Float64}, Float64, 3}, 3}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Optimization.var"#69#86"{Int64}, Float64}, Float64, 3}}}}, Vector{Optimization.var"#69#86"{Int64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, x::Vector{Float64})
@ OptimizationMOI ~/.julia/packages/OptimizationMOI/ZoxrV/src/OptimizationMOI.jl:78
[7] eval_objective(model::Ipopt.Optimizer, x::Vector{Float64})
@ Ipopt ~/.julia/packages/Ipopt/rQctM/src/MOI_wrapper.jl:514
[8] (::Ipopt.var"#eval_f_cb#1"{Ipopt.Optimizer})(x::Vector{Float64})
@ Ipopt ~/.julia/packages/Ipopt/rQctM/src/MOI_wrapper.jl:597
[9] _Eval_F_CB(n::Int32, x_ptr::Ptr{Float64}, x_new::Int32, obj_value::Ptr{Float64}, user_data::Ptr{Nothing})
@ Ipopt ~/.julia/packages/Ipopt/rQctM/src/C_wrapper.jl:38
[10] IpoptSolve(prob::Ipopt.IpoptProblem)
@ Ipopt ~/.julia/packages/Ipopt/rQctM/src/C_wrapper.jl:442
[11] optimize!(model::Ipopt.Optimizer)
@ Ipopt ~/.julia/packages/Ipopt/rQctM/src/MOI_wrapper.jl:727
[12] __solve(prob::SciMLBase.OptimizationProblem{true, SciMLBase.OptimizationFunction{true, Optimization.AutoForwardDiff{nothing}, typeof(model_loss), Nothing, Nothing, Nothing, typeof(model_constraints!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Nothing, Vector{Float64}, Vector{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, opt::Ipopt.Optimizer; maxiters::Nothing, maxtime::Nothing, abstol::Nothing, reltol::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OptimizationMOI ~/.julia/packages/OptimizationMOI/ZoxrV/src/OptimizationMOI.jl:329
[13] __solve
@ ~/.julia/packages/OptimizationMOI/ZoxrV/src/OptimizationMOI.jl:276 [inlined]
[14] #solve#540
@ ~/.julia/packages/SciMLBase/QqtZA/src/solve.jl:84 [inlined]
[15] solve
@ ~/.julia/packages/SciMLBase/QqtZA/src/solve.jl:78 [inlined]
[16] model_fit(u0::Vector{Float64}, data::Vector{Float64})
@ Main my_model.jl:54
[17] top-level scope
@ my_model.jl:77
in expression starting at my_model.jl:58
This bit shows the moment the constraint was violated:
[ Info: Must be NEGATIVE: -283.8381149547038
[ Info: Must be NEGATIVE: 32.36098192303031
ERROR: LoadError: DomainError with -2.4941978436429695
- The first two lines are values of the constraint
a - 1/dt, which must be negative. - However, the second line shows that it became positive!
- The moment the constraint becomes positive, one of the variances becomes negative:
-2.4941978436429695, but it’s not possible to take the log of a non-complex negative value, so I get the error.
Question
What am I doing wrong here? I must be setting up constraints incorrectly, but it seems exactly the same as in Optimization.jl’s documentation.
How to make IPOPT always stay inside the constrained region? I can’t step outside this region because otherwise the variance will become negative.
Package versions
- Julia v1.8.3
- Optimization.jl v3.10.0
- OptimizationMOI 0.1.5
- Ipopt v1.1.0
IPOPT is interior wrt the primal bound constraints (x and the slacks s), as well as the Lagrange multipliers of the bound constraints.