# Modeling min() in JuMP at scale: SOS1 stops converging, Big-M is numerically unstable
**URL:** https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251
**Category:** Optimization (Mathematical)
**Tags:** jump, gurobi
**Created:** [December 1, 2025, 11:31am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251 "2025-12-01T11:31:17Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@Bhanuprakash\_Bahoju](https://discourse.julialang.org/u/Bhanuprakash_Bahoju)
#### Post date: [December 1, 2025, 11:31am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/1 "2025-12-01T11:31:17Z")
I’m implementing a formula engine in JuMP where users can define expressions like:
```julia-auto
z = min(x, y)
```
and I translate these automatically into JuMP constraints.
I have tried two formulations:
* * *
## **1. SOS1 formulation**
```julia-auto
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:
```julia-auto
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.
---
### Author:  [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [December 1, 2025, 1:59pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/2 "2025-12-01T13:59:53Z")
Welcome to Discourse,
but I think we lack motivation to implement this:
- 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.
> [@Bhanuprakash\_Bahoju](#):
>
> What is the recommended way in JuMP/SCIP to model `min()` reliably in large models?
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.
> [@Bhanuprakash\_Bahoju](#):
>
> Are there alternative formulations that scale better?
Users should try their best to build tight models themselves.
> [@Bhanuprakash\_Bahoju](#):
>
> In my application, user-provided values don’t always have good bounds
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
> **[Constraints - Gurobi Optimizer Reference Manual](https://docs.gurobi.com/projects/optimizer/en/current/concepts/modeling/constraints.html#simple-constraints)**
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)
---
### Author:  [@langestefan](https://discourse.julialang.org/u/langestefan)
#### Post date: [December 1, 2025, 4:10pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/3 "2025-12-01T16:10:04Z")
Can you say something about the variables involved? For example are z, x, y integer or continuous variables? Are you minimizing or maximizing z?
SCIP supports [indicator constraints](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints), that should perform much better than SOS1 or Big-M.
---
### Author:  [@Bhanuprakash\_Bahoju](https://discourse.julialang.org/u/Bhanuprakash_Bahoju)
#### Post date: [December 1, 2025, 5:44pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/4 "2025-12-01T17:44:44Z")
Thank you for the reply.
I learnt some new insights from it.
Will try out the Gurobi’s reference and update you.
---
### Author:  [@Bhanuprakash\_Bahoju](https://discourse.julialang.org/u/Bhanuprakash_Bahoju)
#### Post date: [December 1, 2025, 5:55pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/5 "2025-12-01T17:55:01Z")
> [@langestefan](#):
>
> For example are z, x, y integer or continuous variables?
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).
> [@langestefan](#):
>
> SCIP supports [indicator constraints] that should perform much better than SOS1 or Big-M.
doesn’t indicator constraints use
```julia-auto
MathOptInterface.Bridges.Constraint.IndicatorToMILPBridge
```
underneath?, correct me if I am wrong.
but this throws FiniteDomainError if bounds are not specified for the variables.
---
### Author:  [@langestefan](https://discourse.julialang.org/u/langestefan)
#### Post date: [December 1, 2025, 6:38pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/6 "2025-12-01T18:38:48Z")
> [@Bhanuprakash\_Bahoju](#):
>
> 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).
Why dont you just maximize `min(x,y)` then? Thats convex 🙂
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [December 1, 2025, 8:55pm UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/7 "2025-12-01T20:55:17Z")
Hi @bhanuprakash_bahoju, welcome to the forum 😄
There are a few options to consider.
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.
---
### Author:  [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [December 2, 2025, 12:59am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/8 "2025-12-02T00:59:13Z")
> [@langestefan](#):
>
> Why dont you just maximize `min(x,y)` then? Thats convex 🙂
The objective is a scalar expression, i.e. \mathbb{R}^1. But there might be more than one distinct `z = min(x,y)` constraints.
> [@odow](#):
>
> if everything is continuous, you could use a local nonlinear solver like [Ipopt.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Ipopt) which supports `min(x, y)` directly.
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.
> [@Bhanuprakash\_Bahoju](#):
>
> 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).
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.
---
### Author:  [@langestefan](https://discourse.julialang.org/u/langestefan)
#### Post date: [December 2, 2025, 9:02am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/9 "2025-12-02T09:02:38Z")
> [@WalterMadelim](#):
>
> The objective is a scalar expression, i.e. \mathbb{R}^1. But there might be more than one distinct `z = min(x,y)` constraints.
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…
---
### Author:  [@Bhanuprakash\_Bahoju](https://discourse.julialang.org/u/Bhanuprakash_Bahoju)
#### Post date: [December 2, 2025, 9:36am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/10 "2025-12-02T09:36:10Z")
> [@langestefan](#):
>
> As long as you can maximize the `z`’s, then `z <= min(x, y)` will recover `z == min(x,y)`
That’s a very cool idea, But there can be other nonconvex functions like z == max(x,y) in a model.
---
### Author:  [@langestefan](https://discourse.julialang.org/u/langestefan)
#### Post date: [December 2, 2025, 9:48am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/11 "2025-12-02T09:48:49Z")
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
---
### Author:  [@Bhanuprakash\_Bahoju](https://discourse.julialang.org/u/Bhanuprakash_Bahoju)
#### Post date: [December 2, 2025, 11:13am UTC](https://discourse.julialang.org/t/modeling-min-in-jump-at-scale-sos1-stops-converging-big-m-is-numerically-unstable/134251/12 "2025-12-02T11:13:52Z")
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’
```julia-auto
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.