Support constraint programming in JuMP

[dourouc05]

I'm making this new issue to track the progress on CP support in JuMP. (Brief reminder: started with https://github.com/JuliaOpt/JuMP.jl/issues/2014; then, a few changes in https://github.com/JuliaOpt/JuMP.jl/pull/2051 to make more constraints parse.) I'm opening this issue to discuss the progress on this side (maybe this should be marked as 1.0?).

For now, these things have been implemented:

• Several CP sets, an extension to MOI: ConstraintProgrammingExtensions.jl
• New constraints in JuMP to create CP constraints easily: JuCP.jl; this absolutely requires https://github.com/JuliaOpt/JuMP.jl/pull/2051 to be merged
• A first solver is partially wrapped: CPLEXCP.jl

There are many more details about all this in JuCP, including design considerations and next steps. Some of these steps could be made into GSoC projects (dealing with tests, writing new wrappers, mostly), if it's not too late for inclusion.

---

More complete list of things:

• On MOI's side:
  • [ ] Discuss the best way to deal with functions: https://github.com/jump-dev/MathOptInterface.jl/issues/846#issuecomment-816252103
• On JuMP's side:
  • [x] Strict inequality: already possible. https://github.com/dourouc05/JuCP.jl/blob/0c2d0f359e882baf4ff06d404e47d9ee542b05e0/src/sets.jl#L51-L53
  • [x] New syntax (functions as constraints): @constraint(m, constraint(x, y)). PR #2051.
  • [x] New syntax (reification): @constraint(m, x := { constraint(x, y) }), @constraint(m, x := { x < y }). PR #2228 (split from #2051).
  • [ ] New syntax (new expressions in very limited cases): @constraint(m, x == count(y .== 1)). PR #2229 or PR #2241 (split from #2051).
  • [ ] New syntax (new expressions in less limited cases): @constraint(m, x == y[z]) with either y or z variables. No PR for now. Is an overload of getindex enough?
  • [ ] New syntax (Boolean algebra): @constraint(m, x \land y), @constraint(m, x & y). PR #2241.
  • [ ] New syntax (Boolean expressions as constraints): @constraint(m, x) to be understood as @constraint(m, x == true). PR #2051?
  • [ ] Documentation for these new PRs on JuMP's side.
  • [ ] Bug fixes: #2246 (hack in JuCP).
• On CPE's side to extend MOI:
  • [x] Many new sets implemented, review process started: https://github.com/dourouc05/ConstraintProgrammingExtensions.jl
  • [x] Discussion about Element: https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/3
  • [x] Discussion about BinPacking: https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/4
  • [ ] Discussion about constraint mixing: https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/6
  • [ ] FlatZinc:
    • [ ] Output format
      • [x] Standard sets
      • [ ] Nonlinear functions
    • [ ] Input format
      • [x] Standard sets
      • [ ] Nonlinear functions
    • [x] Optimizer interface
  • [ ] XCSP3:
    • [ ] Output format
      • [ ] Standard sets
      • [ ] Nonlinear functions
    • [ ] Input format
      • [ ] Standard sets
      • [ ] Nonlinear functions
    • [ ] Optimizer interface
• On JuCP's side to extend JuMP:
  • [ ] Provide basic, nice syntax for the wrapped constraints:
    • [ ] alldifferent
    • [ ] element
    • [ ] TODO
  • [ ] Provide more natural syntax for constraints that could benefit from it, like element
    • [ ] element
    • [ ] TODO
• On solvers' side:
  • [x] Write generic tests for solvers, specifically for constraint engines (bad support for purely continuous problems).
  • [ ] Wrap CPLEX CP Optimizer in CPLEXCP.jl
    • [x] Base wrapper around Java API
    • [x] CP sets
    • [ ] Support for intervals and interval sequences
  • [ ] Wrap Choco
    • [ ] Low-level API
    • [ ] MOI API
  • [ ] Wrap Chuffed
    • [x] MOI API
  • [ ] Wrap Gecode
    • [ ] Low-level API
    • [ ] MOI API
  • [ ] Wrap JaCoP
    • [ ] Low-level API
    • [ ] MOI API
  • [ ] Wrap LocalSolver (without the usual limitations of CP solver)
    • [ ] Low-level API
    • [ ] MOI API

[dourouc05]

New PR #2241 that ought to replace #2229 to add more flexibility in rewriting expressions.

[dourouc05]

I have more general questions that would need a larger audience:

• https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/3
• https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/4

Any comment welcome :)!

[dourouc05]

A version version of the CPLEX CP Optimizer wrapper is being tagged: https://github.com/JuliaRegistries/General/pull/28051. However, it does not meet all guidelines, but that's on purpose (name close to CPLEX.jl, of course; not importable, because the CI environment doesn't have CPLEX; https://github.com/JuliaRegistries/General/pull/28051/checks?check_run_id=1711999279 is unrelated, I think).

[dourouc05]

A few notes on the recent developments. CPLEXCP.jl has a version 0.1.0 tagged; it should really be useful in many cases. Very specific features such as tuning the branching process are not available at all, but the modelling side is quite advanced (except where Concert uses complex-typed variables). There are quite a few remaining questions to answer, have a look at https://github.com/jump-dev/MathOptInterface.jl/issues/1253 (more feedback is welcome there). Similarly, ConstraintProgrammingExtensions.jl has a 0.1.0 tagged, and should be useful to implement other solver wrappers (https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/7).

For longer-term issues, I also added a paragraph to highlight the fact that CP could really benefit from the new NLP features, once they are developed: https://docs.google.com/document/d/1SfwFrIhLyzpRdXyHKMFKv_XKeyJ6mi6kvhnUZ2LEEjM/edit#, under Use Cases (from https://github.com/jump-dev/MathOptInterface.jl/issues/846).

[odow]

To update this issue:

• @dourouc05 wrote a blog post https://jump.dev/blog/constraint-programming-update/
• Development effort should go into https://github.com/dourouc05/ConstraintProgrammingExtensions.jl
• Extending syntax via parse_constraint is now documented: https://jump.dev/JuMP.jl/dev/developers/extensions/#Parse

[odow]

To update this issue:

• MOI v1.4.0 includes a number of new constraint programming sets
• MOI v1.6.0 includes a number of new bridges from constraint programming to MILP formulations

So this mainly needs some documentation to advertise the new features.

[odow]

Adding much better documentation here: https://github.com/jump-dev/JuMP.jl/pull/3202.

I don't know if we'll ever get to a point of supporting a first-class constraint programming interface. Dealing with logical operators in a functional way seems hard.

The best we could probably do is to support them in the nonlinear interface, and then have something like MiniZinc.jl try to parse them, and error if it encounters some nonlinearity that it doesn't support (like sin).

[chriscoey]

That would be awesome! Would we wait for the new nonlinear interface before trying this?

[odow]

Would we wait for the new nonlinear interface before trying this?

I guess, technically, we could attempt this at present, just by parsing the Julia AST.

I think the main issue is that something like @NLconstraint(model, x || y == 1) parses as x || (y == 1), not (x || y) == 1.

[chriscoey]

The way it parses is the way I'd interpret it actually. We can think of constraints as expressions that must evaluate to true, in which case the above says either x (bool variable) is true, or y (int variable) is 1. We would want to support a new Bool type variable (not Bin) I think.

[odow]

The way it parses is the way I'd interpret it actually

The problem currently is that nonlinear constraints must be one- or two-sided (in)equalities. Let me have a play. I think originally I didn't go down this track because we were considering a fully fledged structure, but maybe a mix of function-in-set for things like all-different, and boolean relations gets us most of the way there.

[chriscoey]

Maybe there could be a new type of constraint, a logical constraint, that is similar to nlconstraint but doesn't have this 1/2 sided inequality restriction.

[chriscoey]

(let's start by dreaming big, and scale back if there's a clear reason why it's too hard 😁)

[odow]

Hah. This was easier than I thought. I'll make a PR.

julia> using JuMP

julia> import MiniZinc

julia> model = Model(() -> MiniZinc.Optimizer{Int}(MiniZinc.Chuffed()))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: MiniZinc

julia> set_optimizer_attribute(model, "model_filename", "test.mzn")

julia> @variable(model, x, Bin)
x

julia> @variable(model, y, Bin)
y

julia> @constraint(model, [x, y] in MOI.AllDifferent(2))
[x, y] ∈ MathOptInterface.AllDifferent(2)

julia> @NLconstraint(model, (x || y) == true)
(x || y) - 1.0 = 0

julia> @NLconstraint(model, (x && y) == false)
(x && y) - 0.0 = 0

julia> optimize!(model)
Warning: included file "count.mzn" overrides a global constraint file from the standard library. This is deprecated. For a solver-specific redefinition of a global constraint, override "fzn_<global>.mzn" instead.

Warning: included file "global_cardinality_low_up.mzn" overrides a global constraint file from the standard library. This is deprecated. For a solver-specific redefinition of a global constraint, override "fzn_<global>.mzn" instead.


julia> value(x), value(y)
(0.0, 1.0)

julia> print(read("test.mzn", String))
include "alldifferent.mzn";
var bool: x;
var bool: y;
constraint alldifferent([x, y]);
constraint ((x \/ y) - 1) == false;
constraint ((x /\ y) - 0) == false;
solve satisfy;

[odow]

See https://github.com/jump-dev/MiniZinc.jl/pull/21

[odow]

So I'm pretty happy with that PR. It's a little tacky, but we could make it work much nicer if we added something like https://github.com/jump-dev/JuMP.jl/pull/3106.

It's also quite a reasonable approach: use the structured function-in-set approach for combinatorial structure like AllDifferent, and use the nonlinear expression graph for a limited set of nonlinear operators (mainly Boolean operators, plus +, -, etc.) We'll loose the ability to rewrite the boolean expressions via bridges, but we might be able to write a layer that transform a ScalarNonlinearFunction that represents a boolean expression into MILP. And that would pave the way for something like a DCP transformation of ScalarNonlinearFunction into conic. (It'd be opt-in by the user instead of automatic, but that's okay.)

[odow]

Updated documentation is https://jump.dev/JuMP.jl/dev/tutorials/linear/constraint_programming/

[blegat]

We'll loose the ability to rewrite the boolean expressions via bridges, but we might be able to write a layer that transform a ScalarNonlinearFunction that represents a boolean expression into MILP. And that would pave the way for something like a DCP transformation of ScalarNonlinearFunction into conic. (It'd be opt-in by the user instead of automatic, but that's okay.)

Not sure I understand this sentence. Aren't you describing a bridge from ScalarNonliearFunction into structured boolean constraints ? Why does it have to be opt-in and not selected by the bridge hyper-graph ?

[dourouc05]

I think the main issue is that something like @NLconstraint(model, x || y == 1) parses as x || (y == 1), not (x || y) == 1.

Given that, in Julia, || only represents the short-circuiting or, I believe it's totally expected and having another type of variable shouldn't help: it would only confuse typical Julia users. However, Julia offers |, the bitwise or: wouldn't that be better?

Another reason why the current parse seems totally fine: what about using || for disjunctions?

For your original code sample, shouldn't it actually mean @NLconstraint(model, x == 1 || y == 1)?

[odow]

julia> dump(:(x || y == 1))
Expr
  head: Symbol ||
  args: Array{Any}((2,))
    1: Symbol x
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol ==
        2: Symbol y
        3: Int64 1

It lowers as x || (y == 1). (x == 1) || (y == 1) is something different:

julia> dump(:(x == 1 || y == 1))
Expr
  head: Symbol ||
  args: Array{Any}((2,))
    1: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol ==
        2: Symbol x
        3: Int64 1
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol ==
        2: Symbol y
        3: Int64 1

The MiniZinc PR is now working for expressions like x || (b && (y < 5)) (none of the brackets are needed, but useful to show precedence.

julia> dump(:(x || (b && (y < 5))))
Expr
  head: Symbol ||
  args: Array{Any}((2,))
    1: Symbol x
    2: Expr
      head: Symbol &&
      args: Array{Any}((2,))
        1: Symbol b
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol <
            2: Symbol y
            3: Int64 5

julia> dump(:(x || (b && y < 5)))
Expr
  head: Symbol ||
  args: Array{Any}((2,))
    1: Symbol x
    2: Expr
      head: Symbol &&
      args: Array{Any}((2,))
        1: Symbol b
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol <
            2: Symbol y
            3: Int64 5

julia> dump(:(x || b && y < 5))
Expr
  head: Symbol ||
  args: Array{Any}((2,))
    1: Symbol x
    2: Expr
      head: Symbol &&
      args: Array{Any}((2,))
        1: Symbol b
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol <
            2: Symbol y
            3: Int64 5

So I think this is a pretty viable approach for solvers wanting to support boolean expressions.

[pedromxavier]

It seems reasonable to me what @chriscoey said about having a special kind of constraints to handle logical expressions.

Some time ago I tried to implement something like

# Don't ink adjacent vertices with the same color
@logical_constraint(model, [∀i = 1:n, ∀j = 1:n, ∀k = 1:n], neigh[i, j] ⟹ ¬(v[i, k] ∧ v[j, k]))

This would reduce to MOI.CNF ∈ MOI.SAT or any other normal form-in-sat like pair. My MOI/JuMP expertise was insufficient at the time so it went no further than writing the macro to accept some of the syntax and assemble some custom structs. Integrating with MOI and launching as an extension was still somewhat far away.

PS: Thinking about dreaming big, implementing a JuMP Interface for Z3 would be a fantastic milestone in this track. There is an active Julia wrapper for the solver: Z3.jl.

[odow]

See https://github.com/jump-dev/MiniZinc.jl/pull/22.

I think we can extend support for solvers by having them opt-in to supporting certain functions on a symbolic basis only.

[odow]

https://github.com/jump-dev/JuMP.jl/pull/3206 adds syntax for reified constraints.

So I wonder if we're close to closing this issue. The progress in MiniZinc.jl demonstrates that this is now mostly a solver issue, not a missing syntax problem at the JuMP level.

[chriscoey]

@odow echoing @pedromxavier and myself, I still would like to be able to write a logical constraint as x || y instead of (x || y) == true. How difficult would it be to allow this? Either by generalizing @NLconstraint so it doesn't require two-sided bounds, or by adding a new type of @xxxxconstraint (eventually it will all just be @constraint anyway). If the root node isn't an equality/inequality, then the constraint is assumed to be a logical constraint (must evaluate to true).

[odow]

Let's table this question until we have ScalarNonlinearFunction support. For now we won't be modifying the @NLconstraint macro because we don't want to be locked into supporting the syntax going forward.

My worry is also that @constraint(model, x + y) silently assuming that the user meant @constraint(model, x + y == true) is bug-prone.

However, this syntax is currently available:

julia> @constraint(model, x || y)
ERROR: LoadError: At REPL[15]:1: `@constraint(model, x || y)`: Unsupported constraint expression: we don't know how to parse constraints containing expressions of type :||.

If you are writing a JuMP extension, implement `parse_constraint_head(::Function, ::Val{:||}, args...)

So we can write extensions that convert || and && into ScalarNonlinearFunction-in-EqualTo(true).

[chriscoey]

OK we can table until SNF is in.

My worry is also that @constraint(model, x + y) silently assuming that the user meant @constraint(model, x + y == true) is bug-prone.

Do you mean it can be a problem if the user forgets to type out part of their constraint? I don't think this is big concern - there's lots of ways to make silent mistakes while typing models, and users should look at pretty-printouts of models/constraints to check they are defined correctly.

[blegat]

We could have @constraint(model, <arbitrary NL>) create a ScalarNonlinearFunction-in-MOI.EqualTo{true} and in the bridge from ScalarNonlinearFunction-inMOI.EqualTo{Bool}, we can check that the nonlinear function is a boolean expression. Otherwise, we can throw the appropriate error.

[chriscoey]

I've created a MiniZinc.jl PR https://github.com/jump-dev/MiniZinc.jl/pull/23 to update it to use ScalarNonlinearFunction from https://github.com/jump-dev/MathOptInterface.jl/pull/2059

Currently the available functionality is too restricted, so we'll want to add more examples and more supported atoms and function-in-set combinations.

[odow]

We're fairly close with this:

julia> using JuMP, MiniZinc

julia> model = GenericModel{Int}(() -> MiniZinc.Optimizer{Int}("chuffed"))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: MiniZinc

julia> @variable(model, x[1:3], Bin)
3-element Vector{GenericVariableRef{Int64}}:
 x[1]
 x[2]
 x[3]

julia> @constraint(model, (x[1] || x[2] && x[3]) == false)
(x[1] || (x[2] && x[3])) - 0 = 0

julia> @objective(model, Min, sum(x))
x[1] + x[2] + x[3]

julia> optimize!(model)
Warning: included file "count.mzn" overrides a global constraint file from the standard library. This is deprecated. For a solver-specific redefinition of a global constraint, override "fzn_<global>.mzn" instead.


julia> value.(x)
3-element Vector{Int64}:
 0
 0
 0

julia> model = GenericModel{Int}(() -> MiniZinc.Optimizer{Int}("chuffed"))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: MiniZinc

julia> @variable(model, x[1:3], Bin)
3-element Vector{GenericVariableRef{Int64}}:
 x[1]
 x[2]
 x[3]

julia> op_xor = NonlinearOperator(xor, :⊻)
NonlinearOperator(xor, :⊻)

julia> @constraint(model, op_xor(x[1], x[2] && x[3]) == true)
⊻(x[1], x[2] && x[3]) - 1 = 0

julia> @objective(model, Max, sum(x))
x[1] + x[2] + x[3]

julia> optimize!(model)
Warning: included file "count.mzn" overrides a global constraint file from the standard library. This is deprecated. For a solver-specific redefinition of a global constraint, override "fzn_<global>.mzn" instead.


julia> value.(x)
3-element Vector{Int64}:
 0
 1
 1

I still need to debug the case

julia> using JuMP, MiniZinc

julia> model = GenericModel{Bool}(() -> MiniZinc.Optimizer{Bool}("chuffed"))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: MiniZinc

julia> @variable(model, x[1:3], Bin)
3-element Vector{GenericVariableRef{Bool}}:
 x[1]
 x[2]
 x[3]

julia> @constraint(model, (x[1] || x[2] && x[3]) == true)
ERROR: Constraints of type MathOptInterface.ScalarNonlinearFunction-in-MathOptInterface.EqualTo{Int64} 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{MathOptInterface.Utilities.CachingOptimizer{MiniZinc.Optimizer{Bool}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Bool}}}}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Bool}}}, f::MathOptInterface.ScalarNonlinearFunction, s::MathOptInterface.EqualTo{Int64})
   @ JuMP ~/.julia/packages/JuMP/OUdu2/src/constraints.jl:679
 [3] add_constraint(model::GenericModel{Bool}, con::ScalarConstraint{GenericNonlinearExpr{GenericVariableRef{Bool}}, MathOptInterface.EqualTo{Int64}}, name::String)
   @ JuMP ~/.julia/packages/JuMP/OUdu2/src/constraints.jl:706
 [4] macro expansion
   @ ~/.julia/packages/JuMP/OUdu2/src/macros.jl:1343 [inlined]
 [5] top-level scope
   @ REPL[125]:1

[odow]

We also want @constraint(model, x[1] || x[2] == x[3]) and we want it parsed as @constraint(model, (x[1] || x[2]) == x[3]) (if possible)