Constraint being violated in a nonlinear model

[barpit20]

using JuMP, Ipopt, Random
n = 1_000
rng = MersenneTwister(1234);
data = randn(rng, n)

mle = Model(with_optimizer(Ipopt.Optimizer, print_level = 0));
@NLparameter(mle, problem_data[i = 1:n] == data[i]);
@variable(mle, μ, start = 0.0);
@variable(mle, σ >= 0.0, start = 1.0);
@NLexpression(mle, likelihood, 
(2 * π * σ^2)^(-n / 2) * exp(-(sum((problem_data[i] - μ)^2 for i in 1:n) / (2 * σ^2)))
);
@NLconstraint(mle, μ == σ^2);
@NLobjective(mle, Max, log(likelihood));

optimize!(mle);

Solutions for the above model violate the constraint μ = σ^2:

julia> @show value(μ);
value(μ) = 0.0

julia> @show value(σ)^2;
value(σ) ^ 2 = 1.0

One looking for statuses, we get:

julia> termination_status(mle)
INVALID_MODEL::TerminationStatusCode = 21

julia> primal_status(mle)
UNKNOWN_RESULT_STATUS::ResultStatusCode = 6

[mlubin]

Thanks for the report. This uncovered a serious bug. (We would have caught this before release with #2033.)

[barpit20]

The above model is slightly different from the one in the examples. If we write it as it is written there, the constraint is satisfied.

using JuMP, Ipopt, Random
n = 1_000
rng = MersenneTwister(1234);
data = randn(rng, n)

model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0))
@variable(model, μ, start = 0.0)
@variable(model, σ >= 0.0, start = 1.0)
@NLconstraint(model, μ == σ^2)
@NLobjective(model, Max, n / 2 * log(1 / (2 * π * σ^2)) -
    sum((data[i] - μ)^2 for i in 1:n) / (2 * σ^2)
)
optimize!(model)

julia> @show value(μ);
value(μ) = 0.6228839414610626

julia> @show value(σ)^2;
value(σ) ^ 2 = 0.6228839414610626

[blegat]

Should be fixed in JuMP master (which requires MOI v0.9.1)

[odow]

This is still not fixed.

(tmp) pkg> st
    Status `/private/tmp/Project.toml`
  [b6b21f68] Ipopt v0.6.0
  [4076af6c] JuMP v0.20.0
  [b8f27783] MathOptInterface v0.9.2

julia> value(μ), value(σ)
(0.0, 1.0)

julia> value(μ) == value(σ)^2
false

[mlubin]

There was a test added in jump-dev/MathOptInterface.jl#850. What happened?

[blegat]

Which value of n and data are you using ? It seems the value of μ and σ you get are simply the starting values but I cannot reproduce it.

[odow]

Oh. The original bug is actually fixed.

I had just copy-pasted from the first example with print_level = 0 and didn't check that it actually solved correctly. It's now running into log(0.0) from the exp(-sum(lots of nonnegative numbers)).