I’m implementing a formula engine in JuMP where users can define expressions like:
z = min(x, y)
and I translate these automatically into JuMP constraints.
I have tried two formulations:
1. SOS1 formulation
function minOf2(model, vars...)
z = @variable(model)
δ = @variable(model, [1:length(vars)])
@constraint(model, δ .>= 0)
for (i, v) in enumerate(vars)
@constraint(model, z <= v)
@constraint(model, z + δ[i] == v)
end
@constraint(model, δ in SOS1())
return z
end
This works for small models, but once the number of min() expressions increases (and therefore SOS1 sets grow), SCIP fails to converge or stalls extremely early.
2. Big-M formulation
I also tried the classic binary + Big-M:
z ≤ v[i]
z ≥ v[i] - M * (1 - b[i])
This works only when M is tight.
In my application, user-provided values don’t always have good bounds, and large M values quickly introduce numerical instability in SCIP.
My question
Given that:
SOS1 scaling is poor in SCIP for many min/max operators
Big-M is unstable when M is not tight
I cannot reliably compute tight bounds because user formulas are arbitrary
What is the recommended way in JuMP/SCIP to model min() reliably in large models?
Are there alternative formulations that scale better?
Any practical advice for building a formula engine on top of JuMP would be greatly appreciated.
First, JuMP already has some high-level interfaces like SOS1 etc. And users aren’t well-motivated to add further wrappers.
Second, your constraint z = min(x, y) is nonconvex intrinsically.
You are mentioning large models, so the answer is simply:
There is no such a universally reliable way. You have to explicitly exploit specific structures of your problem. e.g. to decide tight value of big-M’s.
Users should try their best to build tight models themselves.
If so, the user should stop using a blind big-M method.
My overall opinion is: in real-life large scale MILP,s, we should only formulate our problems using standard MILP language (only Ax <= b constraints, excluding others, e.g. SOS1 constraint is considered vague and blind, z = min(x, y) is considered blind, etc.). And we should try our best to make the polyhedral tight. e.g. x <= 1/2, x Int should be tighten to x <= 0, x Int.
One of the most advisable criteria in mathematical programming is:
There is no foolproof method (e.g. you are asking for a foolproof interface here, which is hopeless).
BTW, Gurobi have some off-the-shelf foolproof constraints APIs
including the nonconvex z = min(x, y) you are mentioning here. But whether it is effective under large-scale problems—you can do some tests and give us feedback. (You can call these “general constraints” APIs from Gurobi.jl)
all variables are continuous. and I am not specifying any objective, In which case the solver just tries to solve the model as just doing a calculation instead optimising (I am also unsure that this might cause issues).
First, if everything is continuous, you could use a local nonlinear solver like Ipopt.jl which supports min(x, y) directly. You can try a global nonlinear solver like Gurobi.jl (but it requires a license).
Second, if you can relax z = min(x, y) into z <= min(x, y), you can add two constraints for z <= x and z <= y. (This isn’t applicable if your problem needs z == min or z >= min)
If you need to use SCIP, and you need the z == min, then the SOS1 formulation is okay. It’s hard to say more without a larger reproducible example of your problem. Your model is non-convex, so it might just be very hard to solve.
The objective is a scalar expression, i.e. \mathbb{R}^1. But there might be more than one distinct z = min(x,y) constraints.
not sure if Ipopt excels at it. I think Ipopt is a nonlinear solver that expect “twice continuously differentiable” ideally. LP is nondifferentiable. Although Ipopt can solve LPs, but the performance might be inadequate compared to dedicated LP solvers, e.g. the barrier algorithm in Gurobi.
According to your description, you’re solving an linear inequality system with some nonconvex z = min(x, y) complicating it. If you input that to Gurobi, I think Gurobi is handling a MILP feasibility system. This is valid usage and has no issues. For an MILP feasibility system, Gurobi will first solve the “root” relaxation problem, this phase is convex and thus efficient. There are two outcomes:
the LP relaxation is infeasible
the LP relaxation is feasible
If it encounters the first case, then the problem has been solved—you get an infeasibility certificate. If it’s the second case, then Gurobi will typically start doing branch-and-bound (if no more cutting planes can be added). This stage is not guaranteed to be efficient so you might be waiting for a prohibitively long time without fetching any valid information.
Which would still be convex. As long as you can maximize the z’s, then z <= min(x, y) will recover z == min(x,y). You need bounds from above though otherwise it just blows up…
If you can find convex formulations it will be multiple orders of magnitude faster than anything based on Ipopt can do, so it’s worth exploring at least I think. Convex.jl can help you validate and find convex forms
Yeah, I will definitely explore more and keep posting here about my findings.
But currently I am thinking of implementing a fallback approach, Where models are first solved using SOS1 constraint, and by the end if the number of SOS1 constraints exceed a threshold or, the model is not converging due to them (I have to parse solver logs for this). I will report the diagnostics to the user’
Model is not converging , probable cause can be unbounded non convex functions.
add realistic bounds to x and y.
Then I will use those bounds to create tight bigM formulations.
