Alpine Error

Stuart_Rogers · June 21, 2024, 5:33pm

When I use Alpine to solve an optimization problem, I get the error below.

  Starting bound-tightening
ERROR: LoadError: Constraints of type MathOptInterface.ScalarQuadraticFunction{Float64}-in-MathOptInterface.GreaterThan{Float64} are not supported by the solver.

If you expected the solver to support your problem, you may have an error in your formulation. Otherwise, consider using a different solver.

The list of available solvers, along with the problem types they support, is available at https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] _moi_add_constraint(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, f::MathOptInterface.ScalarQuadraticFunction{Float64}, s::MathOptInterface.GreaterThan{Float64})
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/constraints.jl:1004
  [3] add_constraint(model::Model, con::ScalarConstraint{QuadExpr, MathOptInterface.GreaterThan{Float64}}, name::String)
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/constraints.jl:1036
  [4] macro expansion
    @ ~/.julia/packages/JuMP/7rBNn/src/macros/@constraint.jl:173 [inlined]
  [5] amp_post_convhull_constrs(m::Alpine.Optimizer, λ::Dict{Any, Any}, α::Dict{Any, Any}, monomial_idx::Int64, discretization::Dict{Any, Any})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/multilinear.jl:428
  [6] amp_convexify_quadratic_univariate(m::Alpine.Optimizer, k::Vector{Expr}, λ::Dict{Any, Any}, α::Dict{Any, Any}, discretization::Dict{Any, Any})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/multilinear.jl:87
  [7] amp_post_convhull(m::Alpine.Optimizer; kwargs::@Kwargs{use_disc::Dict{Any, Any}})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/multilinear.jl:19
  [8] amp_post_convhull
    @ ~/.julia/packages/Alpine/2DP5q/src/multilinear.jl:1 [inlined]
  [9] #amp_post_convexification#126
    @ ~/.julia/packages/Alpine/2DP5q/src/bounding_model.jl:51 [inlined]
 [10] amp_post_convexification
    @ ~/.julia/packages/Alpine/2DP5q/src/bounding_model.jl:48 [inlined]
 [11] create_obbt_model(m::Alpine.Optimizer, discretization::Dict{Any, Any}, bound::Float64; kwargs::@Kwargs{})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:293
 [12] create_obbt_model(m::Alpine.Optimizer, discretization::Dict{Any, Any}, bound::Float64)
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:286
 [13] optimization_based_bound_tightening(m::Alpine.Optimizer; use_bound::Bool, time_limit::Float64, kwargs::@Kwargs{})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:118
 [14] optimization_based_bound_tightening
    @ ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:51 [inlined]
 [15] bound_tightening_wrapper(m::Alpine.Optimizer; use_bound::Bool, kwargs::@Kwargs{})
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:17
 [16] bound_tightening_wrapper
    @ ~/.julia/packages/Alpine/2DP5q/src/presolve.jl:13 [inlined]
 [17] presolve(m::Alpine.Optimizer)
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/main_algorithm.jl:241
 [18] optimize!(m::Alpine.Optimizer)
    @ Alpine ~/.julia/packages/Alpine/2DP5q/src/main_algorithm.jl:157
 [19] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:380 [inlined]
 [20] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/MathOptInterface.jl:85 [inlined]
 [21] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/cachingoptimizer.jl:316
 [22] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:595
 [23] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:546
 [24] macro expansion
    @ ~/julia_code/Gen7_BPM/create_model.jl:399 [inlined]
 [25] macro expansion
    @ ./timing.jl:395 [inlined]
 [26] top-level scope
    @ ~/julia_code/Gen7_BPM/create_model.jl:397
in expression starting at /home/stuart/julia_code/Gen7_BPM/create_model.jl:395

mtanneau · June 21, 2024, 9:25pm

Hi, could you share a reproducible example please?

The error message suggest you might be using a solver that doesn’t support the type of problems you’re trying to solve.

Stuart_Rogers · June 21, 2024, 10:44pm

The optimization problem has a linear objective with nonconvex quadratic constraints. The unknowns variables are binary or continuous. Alpine is a global MINLP solver, so it should be able to handle this probelm. The example below generates the error using both the legacy and non-legacy JuMP interfaces, which can be toggled with the Boolean use_legacy_model.

using JuMP

use_legacy_model = false # Whether to use legacy JuMP nonlinear constraints, expressions, and objectives.

N1 = 16
N = 16
zpad = 2
M = zpad*N

s = rand(M,1)
sqrt_s = sqrt.(s)

Br = rand(M, N1)
Bi = rand(M, N1)

Mt = trunc(Int,floor(M/2))+1
Br = Br[1:Mt,:]
Bi = Bi[1:Mt,:]
s = s[1:Mt]
sqrt_s = sqrt_s[1:Mt]
Me = Mt

using Alpine, Ipopt, HiGHS, Juniper
ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
highs = optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false)
juniper = optimizer_with_attributes(
                  Juniper.Optimizer,
                  MOI.Silent() => true,
                  "mip_solver" => highs,
                  "nl_solver" => ipopt,
                  )
model = Model(
                optimizer_with_attributes(
                Alpine.Optimizer,
                "nlp_solver" => ipopt,
                "mip_solver" => highs,
                "minlp_solver" => juniper,
                ),
                )

@variable(model, y[1:N1], Bin)
if(~use_legacy_model)
        @expression(model, yh, y .- 0.5);
        @variable(model, 1 .>= x_r[1:Me] .>= -1)
        @variable(model, 1 .>= x_i[1:Me] .>= -1)
        @constraint(model, 1 .== x_r.^2+x_i.^2)
        @variable(model, spec_r[1:Me])
        @variable(model, spec_i[1:Me])
        @constraint(model, spec_r .== Br * yh - sqrt_s.*x_r)
        @constraint(model, spec_i .== Bi * yh - sqrt_s.*x_i)
        @expression(model, spec_sq_mag[m in 1:Me], spec_r[m]^2 + spec_i[m]^2)

        #@expression(model, gamma, sum(spec_sq_mag[m] for m in 1:Me))

        @variable(model, gamma >= 0)
        @constraint(model, gamma >= sum(spec_sq_mag[m] for m in 1:Me))

        @objective(model, Min, gamma)
else
        @expression(model, yh[n=1:N1], y[n]-.5)
        @variable(model, 1 .>= x_r[1:Me] .>= -1)
        @variable(model, 1 .>= x_i[1:Me] .>= -1)
        @expression(model, spec_r, Br*yh - sqrt_s.*x_r)
        @expression(model, spec_i, Bi*yh - sqrt_s.*x_i)

        @NLexpression(model, spec_sq_mag[m=1:Me], spec_r[m]^2 + spec_i[m]^2)

        @variable(model, gamma >= 0)
        @NLconstraint(model, gamma >= sum(spec_sq_mag[m] for m in 1:Me))
        @objective(model, Min, gamma)
end

optimize!(model) # Solve the JuMP model.

odow · June 21, 2024, 11:21pm

This just looks like a bug in Alpine, where it assumes that the MIP solver supports quadratic constraints.

Since you’ve asked a number of questions relating to this model recently, let’s step back a bit.

What are you trying to achieve? Is this your full model, or only a subset or it?
Do you need a provably optimal solution to the non-convex problem?
How large will N1 and N be?
Do you have access to Gurobi?

I ran this code:

using JuMP
import Gurobi
begin
    N1 = 16
    N = 16
    zpad = 2
    M = zpad * N
    s = rand(M, 1)
    sqrt_s = sqrt.(s)
    Br = rand(M, N1)
    Bi = rand(M, N1)
    Mt = trunc(Int, floor(M / 2)) + 1
    Br = Br[1:Mt, :]
    Bi = Bi[1:Mt, :]
    s = s[1:Mt]
    sqrt_s = sqrt_s[1:Mt]
    Me = Mt
end;
model = Model(Gurobi.Optimizer);
@variables(model, begin
    y[1:N1], Bin
    -1 <= x_r[1:Me] <= 1
     -1 <= x_i[1:Me] <= 1
    spec_r[1:Me]
    spec_i[1:Me]
    gamma >= 0
end);
@expressions(model, begin
    yh, y .- 0.5
    spec_sq_mag, spec_r.^2 .+ spec_i.^2
end);
@constraints(model, begin
    x_r.^2 .+ x_i.^2 .== 1
    spec_r .== Br * yh - sqrt_s.*x_r
    spec_i .== Bi * yh - sqrt_s.*x_i
    gamma >= sum(spec_sq_mag)
end);
@objective(model, Min, gamma);
optimize!(model)

and got

julia> optimize!(model)
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[x86] - Darwin 23.5.0 23F79)

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 34 rows, 85 columns and 612 nonzeros
Model fingerprint: 0x56cf91ae
Model has 18 quadratic constraints
Variable types: 69 continuous, 16 integer (16 binary)
Coefficient statistics:
  Matrix range     [1e-03, 1e+00]
  QMatrix range    [1e+00, 1e+00]
  QLMatrix range   [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [3e+00, 5e+00]
  QRHS range       [1e+00, 1e+00]
Presolve time: 0.00s
Presolved: 154 rows, 153 columns, 851 nonzeros
Presolved model has 34 quadratic constraint(s)
Presolved model has 34 bilinear constraint(s)

Solving non-convex MIQCP

Variable types: 137 continuous, 16 integer (16 binary)

Root relaxation: objective 0.000000e+00, 17 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    0.00000    0   17          -    0.00000      -     -    0s
     0     0    0.00000    0   19          -    0.00000      -     -    0s
... lines omitted ...
H1533161 302241                       0.7527551    0.72356  3.88%   5.6  278s
 1541606 303778    0.72870   66   25    0.75276    0.72363  3.87%   5.5  280s
 1574142 307722    0.75089   75   22    0.75276    0.72393  3.83%   5.5  285s
 1606070 312226    0.72916   61   31    0.75276    0.72419  3.80%   5.5  290s
 1637880 316423    0.74987   67   26    0.75276    0.72446  3.76%   5.5  295s
 1669783 320571 infeasible   70         0.75276    0.72470  3.73%   5.5  300s

I stopped it there. So even Gurobi finds this problem very difficult to solve.

Stuart_Rogers · June 22, 2024, 3:13am

What are you trying to achieve? Is this your full model, or only a subset or it?
This is the full model, except that s, Br, and Bi are constructed elsewhere rather than being random.
Do you need a provably optimal solution to the non-convex problem?
No.
How large will N1 and N be?
N1 = N = 256.
Do you have access to Gurobi?
No. Gurobi may be struggling for this particular instance because s, Br, and Bi are constructed randomly.

Stuart_Rogers · June 22, 2024, 5:18pm

The error went away when I switched Alpine’s MIP solver from HiGHS to COPT, which supports quadratic constraints.

odow · June 22, 2024, 10:34pm

On installation issues

Okay. So you’re having a bunch of problems that all stem from the fact you have installed many packages in the same environment, and they have slight incompatibilities.

You can’t upgrade Bonmin_jll because it depends on Ipopt_jll, and SCIP also depends on Ipopt_jll, but it requires a slightly older version, etc.

I would encourage you to use separate Julia environments to install only the packages you actually need.

For example:

] activate /tmp/ipopt_env
] add JuMP Ipopt

] activate /tmp/scip_env
] add JuMP SCIP

] activate /tmp/bonmin_env
] add JuMP AmplNLWriter Bonmin_jll

On SCIP and Alpine

SCIP and Alpine are not developed or maintained by the core JuMP developers. Alpine has a number of known bugs: Issues · lanl-ansi/Alpine.jl · GitHub, and so does SCIP: Issues · scipopt/SCIP.jl · GitHub.

If they work for your problem, great. If they don’t, there are other options.

The most important point is that SCIP and Alpine are global MINLP solvers. They will try to find the provably optimal solution. If you don’t need that, then a local MINLP solver like Juniper is a good choice.

It seems that you have access to COPT, since that can solve MIQP problems, have you tried just using COPT?

Complex number support

I also realized that you’re really trying to model this:

using JuMP
import Juniper
import Ipopt
import HiGHS
begin
    N1, N, zpad = 16, 16, 2
    Me = trunc(Int, floor(zpad * N / 2)) + 1
    s = rand(Me)
    B = rand(ComplexF64, Me, N1)
end;
model = Model(
    optimizer_with_attributes(
        Juniper.Optimizer,
        "mip_solver" =>
            optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true),
        "nl_solver" =>
            optimizer_with_attributes(Ipopt.Optimizer, MOI.Silent() => true),
    ),
)
@variable(model, y[1:N1], Bin)
@variable(model, -1-1im <= x[1:Me] <= 1+1im, set = ComplexPlane())
@variable(model, spec[1:Me] in ComplexPlane())
@constraint(model, abs2.(x) .== 1)
@constraint(model, spec .== B * (y .- 0.5) - sqrt.(s) .* x)
@objective(model, Min, sum(abs2, spec));
optimize!(model)

Stuart_Rogers · June 23, 2024, 1:34am

COPT generates an error when I try to solve the model constructed with the non-legacy interface. COPT solves convex MIQCP, but not non-convex MIQCP: GitHub - COPT-Public/COPT-Release: COPT (Cardinal Optimizer) release notes and links. This problem is nonconvex.

Thanks for rewriting my model using complex numbers. The matrix B is actually a weighted M-point DFT matrix, where each row of the DFT matrix is weighted by a vector w. Since the input sequence has length N1 <= M, it is zero-padded to length M prior to the weighted DFT matrix, which is why B is comprised of N1 columns rather than M columns. Zero-padding the input sequence has the effect of interpolating the spectrum. Since the input sequence is real-valued and due to conjugate symmetry of the DFT operator, only Me spectral elements need to be computed, which is why B is comprised of Me rows rather than M rows. x is a vector of unit-magnitude complex numbers and sqrt.(s) defines a magnitude profile.

$ julia create_model.jl
Cardinal Optimizer v7.1.3. Build date Apr 29 2024
Copyright Cardinal Operations 2024. All Rights Reserved
Model fingerprint: 8e817c8

Using Cardinal Optimizer v7.1.3 on Linux
Hardware has 24 cores and 48 threads. Using instruction set X86_NATIVE (1)
Minimizing a MIQCP problem

The original problem has:
    34 rows, 85 columns and 612 non-zero elements
    16 binaries
    35 quadratic constraints

Quadratic constraint 2 is not convex
ERROR: LoadError: API call failed: non-convex problem
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] _check_ret
    @ ~/.julia/packages/COPT/mKfAT/src/MOI/MOI_wrapper.jl:166 [inlined]
  [3] _check_ret
    @ ~/.julia/packages/COPT/mKfAT/src/MOI/MOI_wrapper.jl:455 [inlined]
  [4] _check_ret_optimize
    @ ~/.julia/packages/COPT/mKfAT/src/MOI/MOI_wrapper.jl:180 [inlined]
  [5] _optimize!(model::COPT.Optimizer)
    @ COPT ~/.julia/packages/COPT/mKfAT/src/MOI/MOI_wrapper.jl:2857
  [6] optimize!(model::COPT.Optimizer)
    @ COPT ~/.julia/packages/COPT/mKfAT/src/MOI/MOI_wrapper.jl:2895
  [7] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:380 [inlined]
  [8] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/MathOptInterface.jl:85 [inlined]
  [9] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{COPT.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/cachingoptimizer.jl:316
 [10] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:595
 [11] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:546

odow · June 23, 2024, 2:22am

Ah i didnt realize it didnt support nonconvex.

In that case, just use Juniper.

Stuart_Rogers · June 23, 2024, 2:25am

Can expressions support complex numbers? Is the model below incorrect?

@variable(model, y[1:N1], Bin)
@variable(model, -1-1im <= x[1:Me] <= 1+1im, set = ComplexPlane())
@constraint(model, abs2.(x) .== 1)
@expression(model, spec, B * (y .- 0.5) - sqrt.(s) .* x)
@objective(model, Min, sum(abs2, spec))
optimize!(model)

odow · June 23, 2024, 2:27am

It looks correct, but as with everything, you should double check with appropriate tests on small instances, etc.

Youll probably find that this formulation does not improve the performance. Its usually better to have more linear equality constraints and a sparser quadratic objective.

Stuart_Rogers · June 24, 2024, 4:48pm

SCIP.jl is able to solve both complex number models (your version with more variables and constraints and my version with more expressions), which both have quadratic constraints and objectives due to abs2. I had thought SCIP.jl could only solve models created with the legacy interface, which would require @NLconstraint and an epigraph formulation of the objective to handle abs2.

odow · June 24, 2024, 7:04pm

JuMP has specialised support for quadratic so this doesn’t use any of the non-linear stuff.