Feature request: SOS constraints and piecewise linear functions

[FBumann]

Feature Request: Piecewise Linear Constraints (with SOS2 support)

Motivation

Piecewise linear (PWL) constraints are essential for modeling nonlinear relationships in linear/mixed-integer programs — e.g., cost curves, efficiency curves, fuel consumption, or any function that can be approximated by line segments.

The core use case is linking one or more variables through a piecewise linear relationship: given breakpoints that define a nonlinear curve, constrain variables to lie on that curve. For example, linking a generator's power output to its fuel cost, or linking flow to pressure drop.

Currently, pyoframe has no support for piecewise linear constraints. The underlying solver interface (PyOptInterface) already provides model.add_sos_constraint() for SOS2 constraints, so backend support exists — it just needs to be exposed through pyoframe's API.

Background

Piecewise Linear Functions

A piecewise linear function approximates a nonlinear function $f(x)$ using breakpoints $(x_0, x_1, \ldots, x_n)$ and corresponding function values $(f(x_0), f(x_1), \ldots, f(x_n))$. Between breakpoints, the function is linearly interpolated.

The standard convex combination (SOS2) formulation introduces auxiliary $\lambda$ variables:

• $x = \sum_i \lambda_i \cdot x_i$ (input linking)
• $y = \sum_i \lambda_i \cdot f(x_i)$ (output linking)
• $\sum_i \lambda_i = 1$ (convexity)
• $\lambda$ form an SOS2 set (adjacency — at most two adjacent $\lambda_i$ non-zero)

SOS Constraints

• SOS1 (Type 1): At most one variable in the set can be non-zero.
• SOS2 (Type 2): At most two variables can be non-zero, and they must be adjacent in the ordered set. This is the building block for piecewise linear interpolation.

Proposed API

The primary interface should focus on linking variables through piecewise linear relationships:

import pyoframe as pf
import polars as pl

m = pf.Model()

# --- Scalar example ---
m.x = pf.Variable(lb=0, ub=10)

breakpoints_x = [0, 2, 5, 8, 10]
breakpoints_y = [0, 3, 4, 6, 10]

# Link x and y through a piecewise linear curve
m.pwl = pf.PiecewiseLinear(m.x, m.y, breakpoints_x, breakpoints_y)

# --- Indexed example (e.g., different cost curves per generator) ---
m.power = pf.Variable({"gen": ["g1", "g2"]}, lb=0)
m.cost = pf.Variable({"gen": ["g1", "g2"]}, lb=0)

# Breakpoints can be a DataFrame with different curves per index
breakpoints = pl.DataFrame({
    "gen": ["g1", "g1", "g1", "g2", "g2", "g2"],
    "bp":  [   0,   50, 100,    0,   80, 150],
    "power_bp": [0, 50, 100, 0, 80, 150],
    "cost_bp":  [0, 30,  80, 0, 50, 120],
})

m.cost_curve = pf.PiecewiseLinear(
    {"power": m.power, "cost": m.cost},
    breakpoints,
)

m.minimize = m.cost.sum()

Multi-variable linking

A key capability (as in linopy's design) is linking multiple variables simultaneously through shared interpolation weights. This is essential when a single underlying parameter (e.g., operating point) determines multiple outputs (e.g., both cost and emissions).

# Shared lambda weights link power, cost, AND emissions together
m.curves = pf.PiecewiseLinear(
    {"power": m.power, "cost": m.cost, "emissions": m.emissions},
    breakpoints,  # contains columns for each linked variable
)

Formulation Methods

Different formulation methods suit different problem structures (cf. linopy PR #576):

Method	Description	Requirements	Variables
SOS2 (convex combination)	General PWL using SOS2	Any breakpoint order	$\lambda$ continuous + SOS2
Incremental (delta)	Pure LP formulation	Strictly monotonic breakpoints	$\delta$ continuous only
Disjunctive	Disconnected segments	Forbidden zones / gaps	Binary + $\lambda$ per segment

The incremental method is particularly attractive because it requires no SOS2 or binary variables — it's a pure LP formulation that works when breakpoints are monotonic. An method="auto" option could auto-select based on breakpoint structure.

Low-level SOS Support

While the main goal is the piecewise linear API, exposing raw SOS1/SOS2 constraints is also useful for advanced users:

# SOS2: at most two adjacent x[i] are non-zero
m.sos2 = pf.SOS2(m.x, weights=[1, 2, 3, 4, 5])

# SOS1: at most one x[i] is non-zero
m.sos1 = pf.SOS1(m.x)

Design considerations:

• SOS constraints are fundamentally different from linear/quadratic constraints (no LHS/RHS/sense) — they may need a distinct class.
• Solver support: Gurobi and COPT support SOS natively; HiGHS has limited support; Ipopt does not.

Implementation Notes

• PyOptInterface already supports model.add_sos_constraint(variables, SOSType, weights).
• A supports_sos capability flag should be added to the _Solver dataclass in _constants.py.
• The incremental formulation could be implemented without any new solver features (pure LP), making it available on all solvers including HiGHS.
• Dimensioned piecewise linear constraints need careful design — breakpoints may vary per index (e.g., different cost curves per generator).

Reference

• linopy PR #576: Piecewise Linear Constraints — comprehensive implementation with SOS2, incremental, and disjunctive methods
• PyOptInterface SOS constraint docs
• Wikipedia: Special Ordered Set

[FBumann]

@staadecker First of all thank you for creating this package. Looks amazing. Im currently using linopy and im thinking about exploring the up and downsides of uding tables instead of tensors as the fundamental data model.
I will find some limitations with pyoframe compared to linopy for me.
I'll list them here to give you some feedback. Ill probably open PR's for most of them too.

1. Variable and constraint names: I find it restricting to not be able to freely name variables (no special characters, no white space etc.) Add dict-style component access to Model #226
2. ...

[staadecker]

@FBumann thank you for trying out pyoframe! I'm glad you're enjoying it. Don't hesitate to post in the discussion page if you have questions as you get started.

I've moved your comment on special characters to a new issue (#227) so that we can continue the discussion there. If you find other points of feedback, don't hesitate to open a new issue!

As for SOS constraints can you explain your use case. What are you are trying to accomplish and in which context? This will help me assign a priority for this issue.

Looking forward to collaborating!

[FBumann]

@staadecker I like your organized approach on this. Seems like good documenting of decisions.
First of all, im currently exploring the possibilities of porting flixopt from linopy to pyoframe.
Im exploring if the tabular datamodel fits my needs better than xarray (first steps here: https://github.com/FBumann/fluxopt)
While exploring this, i found limitations compared to linopy for my usecase. SO its nothing urgent, take it more as feedback of a potential user, who maybe just didnt adjust his workflow to pyoframe.
About SOS-constraints:
They are the most efficient way to model some milp relations, especially piecewise linear functions.
My use case is about modeling investment costs or efficiency curves as a stepwise function in milp.
Modeling them by hand is possible, but having the underlying modeling framework providing them saves a lot of work for downstream users. Thats why i opened a PR for linopy, and if i would port to pyoframe, id like to add it here too.

See:
PyPSA/linopy#576

[staadecker]

@FBumann thank you for this helpful context. I'm very excited to see how you find Pyoframe compared to Linopy for flixopt.

Regarding SOS constraints I'm open to adding support for such capabilities since it seems to be a relatively often used feature and as you noticed PyOptInterface already supports it. It seems that there are two feature requests here:

1. Ability to create SOS constraints.
2. Some higher-level abstraction for piecewise linear constraints that reformulates the constraints as either a) an LP formulation, b) a SOS2 constraint formulation, c) a binary reformulation.

Here are some messy preliminary thoughts and questions:

SOS constraints

• COPT and Gurobi but not HiGHS supoprt SOS constraints. This is fine enough, although it is good to keep in mind that support for SOS constraints would only be available for academics or paying users.

• We need to think about an appropriate way for creating lists of variables over multiple dimensions. One possible syntax similar to Linopy:

model.my_sos_constr = pf.SOS1(m.VarX, sos_dim="dim1")

• In practice, the first argument should of SOS1 must be of type pf.Expression (since if we want to manipulate a variable it will have to become an expression). So how do we handle expressions with sums? Also, should sos_dim accept multiple dimensions?

• Is it ever desirable to create SOS constraints over different sets of variables (e.g. m.VarX and m.VarY)? Presumably yes. If so, a list of expressions could be accepted (pf.SOS1(m.VarX, m.VarY, sos_dim="dim1")). Alternatively, a (nicer?) syntax could be to let users create SOS constraints from Pyoframe expressions where all the variables that are summed are put into the same SOS constraint. This seems the most flexible but perhaps not intuitive?

model.my_sos_constr = pf.SOS1((m.VarX + m.VarY).sum("dim1"))

• pf.SOS1(expr) or expr.SOS1() ? The later is preferred by users for summation...

• We need to think about an appropriate way to represent the order of the variables for SOS2 constraints in a DataFrame based fashion. Perhaps we need to add a sort function to Pyoframe objects to let users control that order? Alternatively, if we're basing SOS2 constraints from expressions the order could be based on the coefficients (just like PyOptInterface and Gurobi use "weights" as an order). But this feels un-intuitive.

• For all of these, step 1 should be to find and examine several real use cases for SOS constraints and their context. This will help make better design decisions. FlixOpt seems to have a few although it's not clear to me which are SOS2 constraints vs which would be best suited as binary or linear formulations. Linopy also has some examples!

Piecewise abstraction

• If there is a clear way to implement this abstraction even just for the incremental (LP) reformulation that feels like a helpful addition to the library. I remember using an incremental formulation of piecewise constraints when working on energy systems models. Similar questions arise as to what type of datastructure should Pyoframe accept?

• It appears that Gurobi has found that it is often more efficient to use other reformulations instead of SOS2 (and Gurobi will do these automatically if you don't). So is it worth implementing an SOS2 reformulation? I'm sure there are use cases where the answer is yes, but what are those use cases? (Conversely, is it worth implementing a binary reformulation if Gurobi already implements one. Being "smart" about when to reformulate seems like a task for solvers not Pyoframe.)

• More generally, I'd like to understand a) when is the SOS2 reformulation preferred and the LP reformulation cannot be used? (U-shaped convex curves?) b) when is the binary formulation preferred and the SOS2 reformulation not possible? etc.