feat: Add sos constraints

[coroa]

Closes #437 (if applicable).

Changes proposed in this Pull Request

Add Special Ordered Sets (SOS) Constraints Support

I added SOS1/2 constraints as an extension of variables. ie. in the example

import linopy

m = linopy.Model()
# Create variables for facility locations
locations = pd.Index([0, 1, 2], name="locations")
build = m.add_variables(coords=[locations], name="build", binary=True)

# Add SOS1: at most one facility can be built
m.add_sos_constraints(build, sos_type=1, sos_dim="locations")

The add_sos_constraints only adds sos_type and sos_dim as attributes to the build variable's dataset.

These attributes have then to be interpreted by the lp file writer or the direct apis. In this draft version this is implemented for the gurobipy direct api as well as for the lp-polars writer.

The coordinates of the chosen dimension are also always used as weights for the SOS. It felt natural, but i am open to other suggestions. Weights are only used by SOS2 constraints to order them, and then only two adjacent variables can be non-zero at the same time. In principle always taking the coordinate order directly (ie use something like np.arange(var.sizes[sos_dim]) as weights) would be the another option.

The documentation that i added was claude written and might still contain mistakes. The code is tested, but the claude tests are to extensive to be helpful. I need to clean them up first.

Core SOS Constraints API

• Model.add_sos_constraints(): New method to add SOS1 and SOS2 constraints to variables
• Multi-dimensional support: SOS constraints work with multi-dimensional variables, creating constraints for each combination of non-SOS dimensions
• Attribute-based tracking: SOS constraints are stored as variable attributes (sos_type, sos_dim)
• Variables.sos property: Easy access to all SOS-constrained variables in a model

SOS Constraint Types

• SOS1 (Type 1): At most one variable in the ordered set can be non-zero
• SOS2 (Type 2): At most two adjacent variables (based on ordering weights) can be non-zero

Solver Support

• File-based solvers: LP file format output with SOS constraint sections
• Gurobi direct API: Native SOS constraint support via gurobipy

📁 Files Added/Modified

Core Implementation

• linopy/model.py: Added add_sos_constraints() method with validation
• linopy/variables.py: Enhanced Variable class with SOS attribute support and Variables.sos property
• linopy/io.py: SOS constraint export for file-based solvers and Gurobi direct API

Documentation & Examples

• doc/sos-constraints.rst: Comprehensive documentation with theory and examples
• doc/index.rst: Updated to include SOS constraints documentation

Testing

to be done

Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• Unit tests for new features were added (if applicable).
• A note for the release notes doc/release_notes.rst of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.

[coroa]

Note that with the implementation here, it is currently unclear how to remove an SOS constraint again or where to access dual values. Is that a thing for SOS constraints?

[FabianHofmann]

Impressive and wonderful! the IO currently only support the polars lp file writing, not the "standard" pandas based, right? for non-breaking compat we would then need to focus on the IO

• export to highspy & gurobipy (and/or raise an error if we don't support for now)
• export with standard lp file (generally also happy to skip, in case we merge in #366 after merging #485 and making a final benchmark)

[fneum]

Great, I will test out now!

Note that with the implementation here, it is currently unclear how to remove an SOS constraint again or where to access dual values. Is that a thing for SOS constraints?

No need to worry about duals since it's MILP class.

Option to remove SOS constraints would be useful.

the IO currently only support the polars lp file writing, not the "standard" pandas based, right?

Would be cool if the standard pandas lp file writing were supported. But probably @coroa just sketched it with polars first.

[coroa]

Indeed, i was only sketching out how it works (it's still a draft).

Currently supported:

• lp-polars
• direct (for gurobipy)

I don't think it's difficult to do the others.

[coroa]

For removing, it would be easy to add:

m.remove_sos_constraints(variable)

[coroa]

Open question:
Should we use the sos_dim coordinates as weights of the SOS constraint or simply do not add weights?

Context: The weight of an SOS constraint is only relevant for ordering the variables in an SOS2 constraint (where 2 adjacent variables can be non-zero).

If we use the coordinate values:

1. one can modify the ordering away from the coordinate order (by using a non-monotonuous coordinate)
2. string coordinates like location names are not allowed (since they do not make good weights).

The current implementation uses coordinate values, but the more i think about it, i think this is more confusing than helpful (i like string coordinates, sometimes) and using the coordinate order directly seems to be intuitive. if you want a different order, it would be fine to reindex the dimension.

Any other thoughts? Am I misunderstanding SOS weights?

[coroa]

2nd open question: should we add a check/raise when the solver does not support sos constraints in the lp file somewhere?

[fneum]

Should we use the sos_dim coordinates as weights of the SOS constraint or simply do not add weights?

It would certainly be more general. If I understand correctly, it's only relevant for SOS2 since SOS1 does not care about order. Mostly, SOS2 would be used for interpolation, in which case the coordinate would very likely be a numerical value. But I think it's better to be explicit.

In Gurobi reference you can find:

An SOS constraint is described using a list of variables and a list of corresponding weights. While the weights have historically had intuitive meanings associated with them, we simply use them to order the list of variables. The weights should be unique.

https://docs.gurobi.com/projects/optimizer/en/current/concepts/modeling/constraints.html

2nd open question: should we add a check/raise when the solver does not support sos constraints in the lp file somewhere?

Yes, definitely an error should be thrown when an attempt is made to solve a model with SOS constraints when either the solver or the IO API (or the combination) does not support it.

Other Feedback

• I had "highs" with "lp-polars" fail when there were no other constraints than SOS.

• A way of accessing and representing SOS constraints would be nice (like m.constraints).

Here are the cases I played with (could be unit tests):

import linopy
import pandas as pd
import numpy as np

m = linopy.Model()
locations = pd.Index([0, 1, 2], name="locations")
build = m.add_variables(coords=[locations], name="build", binary=True)
m.add_sos_constraints(build, sos_type=1, sos_dim="locations")
m.add_objective(build * [1, 2, 3], sense="max")
m.solve(solver_name='gurobi', io_api='lp-polars')
assert np.isclose(build.solution.values, [0, 0, 1]).all()
assert np.isclose(m.objective.value, 3)

m = linopy.Model()
locations = pd.Index([0, 1, 2], name="locations")
build = m.add_variables(coords=[locations], name="build", binary=True)
m.add_sos_constraints(build, sos_type=2, sos_dim="locations")
m.add_objective(build * [1, 2, 3], sense="max")
m.solve(solver_name='gurobi', io_api='direct')
assert np.isclose(build.solution, [0, 1, 1]).all()
assert np.isclose(m.objective.value, 5)

m = linopy.Model()
locations = pd.Index([0, 1, 2], name="locations")
build = m.add_variables(coords=[locations], name="build", binary=True)
m.add_sos_constraints(build, sos_type=2, sos_dim="locations")
m.add_objective(build * [2, 1, 3], sense="max")
m.solve(solver_name='gurobi', io_api='direct')
assert np.isclose(build.solution, [0, 1, 1]).all()
assert np.isclose(m.objective.value, 4)

[coroa]

1. weights: Thanks for this quote of the gurobi docs, which confirms that the weights do not have an inherent meaning any more. For avoiding surprises and for ease of use, i would suggest to use the dimension order directly, since we already tied it now to a dimension (with an order), ie. i'll remove the weights code again.

2. raises/checks: i'll add the remaining io_api and solver implementations, but not for the standard lp writer due to Replace pandas-based LP file writing with polars implementation #496 . Can anyone try to find out which solvers do and don't support it?

3. representation like m.constraints. hmm. yes, that is what i meant what is tricky. in this implementation sos constraints are a variable quality and you get all variables with an sos constraint with: m.variables.sos, which currently looks like:

linopy.model.Variables
----------------------
 * build (locations)

i guess we could improve this representation to:

linopy.model.Variables
----------------------
 * build (locations) - sos1 on locations

[fneum]

representation like m.constraints. hmm. yes, that is what i meant what is tricky. in this implementation sos constraints are a variable quality and you get all variables with an sos constraint with: m.variables.sos, which currently looks like:

I guess this is fine as long as it is possible to distinguish between SOS1 and SOS2 sets.

[coroa]

Rebased on latest master (ie. where only lp-polars remains)

I added the representation that we talked about above. I added the m.remove_sos_constraints(variable).

Solver support seems to be:

1. cbc - probably (it does have an addSOS function in the c interface docs, untested)
2. mosek - explicit no in docs
3. gurobi - yes (tested)
4. cplex - yes (tested)
5. xpress - yes (tested)
6. highs - no (tested)
7. mindopt - no (actually yes, but different format from used lp file format, why??)
8. glpk - no (according to note by Robbie on mailing list)

We unfortunately do not have a good representation of solver features! Do we? Anyone has any prior ideas on how this should look like?