MathOptInterface.OTHER_ERROR when trying to use ISRES of NLopt through JuMP

Specific Domains Optimization (Mathematical)
,
1

I am trying to minimize a nonlinear function with nonlinear inequality constraints with NLopt and JuMP.
In my test code below, I am minimizing a function with a known global minima.
Local optimizers such as LD_MMA fails to find this global minima, so I am trying to use global optimizers of NLopt that allow nonlinear inequality constraintes.

However, when I check my termination status, it says “termination_status(model) = MathOptInterface.OTHER_ERROR”. I am not sure which part of my code to check for this error.
What could be the cause?

I am using JuMP since in the future I plan to use other solvers such as KNITRO as well, but should I rather use the NLopt syntax?

Below is my code:

# THIS IS A CODE TO SOLVE FOR THE TOYMODEL
# THE EQUILIBRIUM IS CHARACTERIZED BY A NONLINEAR SYSTEM OF ODEs OF INCREASING FUCTIONS B(x) and S(y)
# THE GOAL IS TO APPROXIMATE B(x) and S(y) WITH POLYNOMIALS
# FIND THE POLYNOMIAL COEFFICIENTS THAT MINIMIZE THE LEAST SQUARES OF THE EQUILIBRIUM EQUATIONS

# load packages
using Roots, NLopt, JuMP

# model primitives and other parameters
k = .5 # equal split
d = 1 # degree of polynomial
nparam = 2*d+2 # number of parameters to estimate
m = 10 # number of grids
m -= 1
vGrid = range(0,1,m) # discretize values
c1 = 0 # lower bound for B'() and S'()
c2 = 2 # lower and upper bounds for offers
c3 = 1 # lower and upper bounds for the parameters to be estimated

# objective function to be minimized
function obj(α::T...)  where {T<:Real}
    # split parameters
    αb = α[1:d+1] # coefficients for B(x)
    αs = α[d+2:end] # coefficients for S(y)

    # define B(x), B'(x), S(y), and S'(y)
    B(v) = sum([αb[i] * v .^ (i-1) for i in 1:d+1])
    B1(v) = sum([αb[i] * (i-1) * v ^ (i-2) for i in 2:d+1])
    S(v) = sum([αs[i] * v .^ (i-1) for i in 1:d+1])
    S1(v) = sum([αs[i] * (i-1) * v ^ (i-2) for i in 2:d+1])

    # the equilibrium is characterized by the following first order conditions
    #FOCb(y) = B(k * y * S1(y) + S(y)) - S(y)
    #FOCs(x) = S(- (1-k) * (1-x) * B1(x) + B(x)) - B(x)
    function FOCb(y)
        sy = S(y)
        binv = find_zero(q -> B(q) - sy, (-c2, c2))
        return k * y * S1(y) + sy - binv
    end
    function FOCs(x)
        bx = B(x)
        sinv = find_zero(q -> S(q) - bx, (-c2, c2))
        return (1-k) * (1-x) * B1(x) - B(x) + sinv
    end

    # evaluate the FOCs at each grid point and return the sum of squares
    Eb = [FOCb(y) for y in vGrid]
    Es = [FOCs(x) for x in vGrid]
    E = [Eb; Es]
    return E' * E
end

# this is the actual global minimum
αa = [1/12, 2/3, 1/4, 2/3]
obj(αa...)

# do optimization
model = Model(NLopt.Optimizer)
set_optimizer_attribute(model, "algorithm", :GN_ISRES)
@variable(model, -c3 <= α[1:nparam] <= c3)
@NLconstraint(model, [j = 1:m], sum(α[i] * (i-1) * vGrid[j] ^ (i-2) for i in 2:d+1) >= c1) # B should be increasing
@NLconstraint(model, [j = 1:m], sum(α[d+1+i] * (i-1) * vGrid[j] ^ (i-2) for i in 2:d+1) >= c1) # S should be increasing
register(model, :obj, nparam, obj, autodiff=true)
@NLobjective(model, Min, obj(α...))
println("")
println("Initial values:")
for i in 1:nparam
    set_start_value(α[i], αa[i]+rand()*.1)
    println(start_value(α[i]))
end
JuMP.optimize!(model)
println("")
@show termination_status(model)
@show objective_value(model)
println("")
println("Solution:")
sol = [value(α[i]) for i in 1:nparam]

My output:

Initial values:
0.11233072522513032
0.7631843020124309
0.3331559403539963
0.7161240026812674

termination_status(model) = MathOptInterface.OTHER_ERROR
objective_value(model) = 0.19116585196576466

Solution:
4-element Vector{Float64}:
 0.11233072522513032
 0.7631843020124309
 0.3331559403539963
 0.7161240026812674
2

You have multiple issues:

  • range(0,1,m) should be range(0,1; length = m) (how did this work otherwise?)
  • Sometimes your objective function errors because the root doesn’t exist. If I run with Ipopt, I get
    ERROR: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

Here’s what I would do:

using JuMP
import Ipopt
import Roots

function main()
    k, d, c1, c2, c3, m = 0.5, 1, 0, 2, 1, 10
    nparam = 2 * d + 2
    m -= 1
    vGrid = range(0, 1; length = m)
    function obj(α::T...)  where {T<:Real}
        αb, αs = α[1:d+1], α[d+2:end]
        B(v) = sum(αb[i] * v^(i-1) for i in 1:d+1)
        B1(v) = sum(αb[i] * (i-1) * v^(i-2) for i in 2:d+1)
        S(v) = sum(αs[i] * v^(i-1) for i in 1:d+1)
        S1(v) = sum(αs[i] * (i-1) * v^(i-2) for i in 2:d+1)
        function FOCb(y)
            sy = S(y)
            binv = Roots.fzero(q -> B(q) - sy, zero(T))
            return k * y * S1(y) + sy - binv
        end
        function FOCs(x)
            bx = B(x)
            sinv = Roots.fzero(q -> S(q) - bx, zero(T))
            return (1-k) * (1-x) * B1(x) - B(x) + sinv
        end
        return sum(FOCb(x)^2 + FOCs(x)^2 for x in vGrid)
    end
    αa = [1/12, 2/3, 1/4, 2/3]
    model = Model(Ipopt.Optimizer)
    @variable(model, -c3 <= α[i=1:nparam] <= c3, start = αa[i]+ 0.1 * rand())
    @constraints(model, begin
        [j = 1:m], sum(α[i] * (i-1) * vGrid[j]^(i-2) for i in 2:d+1) >= c1
        [j = 1:m], sum(α[d+1+i] * (i-1) * vGrid[j]^(i-2) for i in 2:d+1) >= c1
    end)
    register(model, :obj, nparam, obj; autodiff = true)
    @NLobjective(model, Min, obj(α...))
    optimize!(model)
    print(solution_summary(model))
    return value.(α)
end
main()
1 Like