Generalized Disjunctive Programming (GDP)

Generalized Disjunctive Programming (GDP)#

Many engineering design problems involve discrete choices that activate or deactivate groups of constraints. For example:

  • A chemical plant that can operate in high-throughput or high-purity mode

  • A supply chain where exactly one of several routes must be selected

  • A scheduling problem where tasks run on one of several machines

Generalized Disjunctive Programming (GDP) [Raman and Grossmann, 1994] provides a natural way to model these problems using disjunctions (“either A or B”) and logical propositions (“if X then Y”), rather than manually encoding big-M constraints.

This tutorial covers:

  1. Basic disjunctions with either_or()

  2. Indicator constraints with if_then()

  3. Logical propositions: implies(), iff(), at_least(), at_most(), exactly()

  4. A full superstructure design example

  5. Big-M vs Hull reformulation tradeoffs

import os

os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"

import discopt.modeling as dm
import numpy as np

print("discopt loaded successfully")

discopt loaded successfully

1. Disjunctions with either_or()#

The simplest GDP construct is a disjunction: exactly one group of constraints must hold.

Consider a variable \(x \in [0, 10]\) that must lie in one of two ranges:

\[\begin{bmatrix} x \leq 3 \end{bmatrix} \;\veebar\; \begin{bmatrix} x \geq 7 \end{bmatrix}\]

We want to minimize \((x - 5)^2\), so the optimal solution chooses the range closest to 5.

m = dm.Model("basic_disjunction")
x = m.continuous("x", lb=0, ub=10)

m.minimize((x - 5) ** 2)

# Either x <= 3 OR x >= 7
m.either_or(
    [
        [x <= 3],  # disjunct A
        [x >= 7],  # disjunct B
    ],
    name="excluded_middle",
)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: 4.0)")
print(f"x = {result.x['x']:.6f}  (expected: 3.0 or 7.0)")
print(f"Nodes:     {result.node_count}")

Status:    optimal
Objective: 4.000000  (expected: 4.0)
x = 3.000000  (expected: 3.0 or 7.0)
Nodes:     3

The solver picks \(x = 3\) (or \(x = 7\)) since both give objective value 4.0. Under the hood, either_or() introduces a binary variable for each disjunct and reformulates the constraints using big-M (by default).

Multi-constraint disjuncts#

Each disjunct can contain multiple constraints. Here we model a process unit that operates in one of two modes, each with different flow and quality bounds:

m = dm.Model("two_mode_unit")
flow = m.continuous("flow", lb=0, ub=100)
quality = m.continuous("quality", lb=0, ub=1)

# Maximize revenue = flow * quality
m.maximize(flow * quality)

# Mode A: high flow, low quality
# Mode B: low flow, high quality
m.either_or(
    [
        [flow >= 60, flow <= 100, quality >= 0.2, quality <= 0.4],  # mode A
        [flow >= 10, flow <= 30, quality >= 0.8, quality <= 1.0],  # mode B
    ],
    name="operating_mode",
)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Revenue:   {result.objective:.4f}")
print(f"Flow:      {result.x['flow']:.4f}")
print(f"Quality:   {result.x['quality']:.4f}")
print(f"Mode selected: {'A (high flow)' if result.x['flow'] > 50 else 'B (high quality)'}")

Status:    optimal
Revenue:   30.0000
Flow:      30.0000
Quality:   1.0000
Mode selected: B (high quality)

2. Indicator Constraints with if_then()#

Indicator constraints activate a set of constraints only when a binary variable equals 1:

\[y = 1 \;\Rightarrow\; \text{constraints hold}\]

This is useful for modeling “build or don’t build” decisions. If we build a unit (\(y = 1\)), it must operate within certain limits. If we don’t build it (\(y = 0\)), its flow is zero.

m = dm.Model("indicator_example")

# Two optional processing units
build = m.binary("build", shape=(2,))
flow = m.continuous("flow", shape=(2,), lb=0, ub=50)

# Revenue minus fixed costs
revenue_per_unit = np.array([3.0, 5.0])
fixed_cost = np.array([20.0, 40.0])

m.maximize(
    dm.sum(lambda i: revenue_per_unit[i] * flow[i] - fixed_cost[i] * build[i], over=range(2))
)

# If built, flow must be between 5 and 50
for i in range(2):
    m.if_then(build[i], [flow[i] >= 5, flow[i] <= 50], name=f"unit{i}_active")
    m.subject_to(flow[i] <= 50 * build[i], name=f"unit{i}_bigM")  # flow=0 if not built

# Total flow capacity
m.subject_to(dm.sum(flow) <= 60, name="total_capacity")

result = m.solve()
print(f"Status:    {result.status}")
print(f"Profit:    {result.objective:.4f}")
for i in range(2):
    b = result.x["build"][i]
    f = result.x["flow"][i]
    print(f"Unit {i}: build={b:.0f}, flow={f:.2f}")

Status:    optimal
Profit:    220.0000
Unit 0: build=1, flow=10.00
Unit 1: build=1, flow=50.00

3. Logical Propositions#

GDP models often include logical relationships between binary decisions. discopt provides five proposition methods:

Method

Meaning

Linear encoding

m.implies(y1, y2)

\(y_1 = 1 \Rightarrow y_2 = 1\)

\(y_1 \leq y_2\)

m.iff(y1, y2)

\(y_1 = 1 \Leftrightarrow y_2 = 1\)

\(y_1 = y_2\)

m.at_least(k, [y1, y2, ...])

At least \(k\) are 1

\(\sum y_i \geq k\)

m.at_most(k, [y1, y2, ...])

At most \(k\) are 1

\(\sum y_i \leq k\)

m.exactly(k, [y1, y2, ...])

Exactly \(k\) are 1

\(\sum y_i = k\)

These encode cleanly as linear constraints on binary variables, avoiding the need for manual big-M reasoning.

m = dm.Model("logical_propositions")

# Four project options
projects = [m.binary(f"p{i}") for i in range(4)]
profit = [10.0, 15.0, 8.0, 12.0]

m.maximize(dm.sum(lambda i: profit[i] * projects[i], over=range(4)))

# Budget: select at most 2 projects
m.at_most(2, projects, name="budget")

# Project 2 requires project 0 (precedence)
m.implies(projects[2], projects[0], name="p2_requires_p0")

# Projects 1 and 3 are equivalent (joint venture: both or neither)
m.iff(projects[1], projects[3], name="p1_iff_p3")

result = m.solve()
print(f"Status:    {result.status}")
print(f"Profit:    {result.objective:.2f}")
for i in range(4):
    selected = result.x[f"p{i}"]
    print(f"  Project {i}: {'SELECTED' if selected > 0.5 else 'not selected'} (profit={profit[i]})")

Status:    optimal
Profit:    27.00
  Project 0: not selected (profit=10.0)
  Project 1: SELECTED (profit=15.0)
  Project 2: not selected (profit=8.0)
  Project 3: SELECTED (profit=12.0)

Since projects 1 and 3 are tied together (iff), selecting both counts as 2 of our budget of 2. The solver picks projects 1 and 3 for total profit of 27.0, which beats any other feasible pair.

4. Full Example: Process Superstructure Design#

A common application of GDP is superstructure optimization: given a network of optional processing units, decide which to build and how to operate them.

Our problem has:

  • 3 processing units, each with two operating modes (high-throughput or high-purity)

  • Binary decisions: which units to build, which mode to use

  • Disjunctive constraints for mode selection (either_or)

  • Indicator constraints for unit activation (if_then)

  • Logical propositions for design rules (implies, at_least, at_most)

m = dm.Model("gdp_superstructure")

n_units = 3

# Binary: whether each unit is built
build = m.binary("build", shape=(n_units,))

# Binary: operating mode for each unit (0 = high-throughput, 1 = high-purity)
mode = m.binary("mode", shape=(n_units,))

# Continuous: flow through each unit
flow = m.continuous("flow", shape=(n_units,), lb=0, ub=100)

# Continuous: product quality from each unit
quality = m.continuous("quality", shape=(n_units,), lb=0, ub=1)

# Revenue and cost data
revenue_price = 50.0
fixed_cost = np.array([100.0, 150.0, 120.0])
variable_cost = np.array([2.0, 3.0, 2.5])

# Objective: maximize revenue - cost
total_revenue = revenue_price * dm.sum(lambda i: flow[i] * quality[i], over=range(n_units))
total_fixed = dm.sum(lambda i: fixed_cost[i] * build[i], over=range(n_units))
total_variable = dm.sum(lambda i: variable_cost[i] * flow[i], over=range(n_units))
m.maximize(total_revenue - total_fixed - total_variable)

# --- Disjunctive constraints (GDP) ---
# Each unit has two operating modes:
#   Mode A (high-throughput): flow in [50, 100], quality in [0.3, 0.5]
#   Mode B (high-purity):    flow in [10, 40],  quality in [0.7, 0.9]
for i in range(n_units):
    m.either_or(
        [
            [flow[i] >= 50, flow[i] <= 100, quality[i] >= 0.3, quality[i] <= 0.5],
            [flow[i] >= 10, flow[i] <= 40, quality[i] >= 0.7, quality[i] <= 0.9],
        ],
        name=f"unit{i}_mode",
    )

# --- Indicator constraints ---
# Unit operates only if built
for i in range(n_units):
    m.if_then(build[i], [flow[i] >= 1], name=f"unit{i}_active")
    m.subject_to(flow[i] <= 100 * build[i], name=f"unit{i}_bigM")

# --- Logical propositions ---
# Unit 2 requires unit 0 (precedence)
m.implies(build[2], build[0], name="unit2_requires_unit0")

# Units 0 and 1 cannot both run in high-purity mode
m.at_most(1, [mode[0], mode[1]], name="purity_exclusion")

# At least 2 units must be built
m.at_least(2, [build[i] for i in range(n_units)], name="min_units")

# At most 2 units in high-throughput mode (mode=0 means high-throughput,
# so at_least 1 must be high-purity: mode=1)
m.at_least(1, [mode[i] for i in range(n_units)], name="min_purity_units")

print(m)

Model: gdp_superstructure
  Variables: 12 (6 continuous, 6 integer/binary)
  Constraints: 13
  Objective: maximize (((50 * Σ[3 terms]) - Σ[3 terms]) - Σ[3 terms])
  Parameters: 0

result = m.solve(time_limit=60)

print(f"Status:    {result.status}")
print(f"Profit:    {result.objective:.2f}")
print(f"Nodes:     {result.node_count}")
print(f"Time:      {result.wall_time:.2f}s")
print()

n_units = 3
for i in range(n_units):
    built = result.x["build"][i]
    f = result.x["flow"][i]
    q = result.x["quality"][i]
    mode_val = result.x["mode"][i]
    mode_str = "high-purity" if mode_val > 0.5 else "high-throughput"
    if built > 0.5:
        print(f"Unit {i}: BUILT, mode={mode_str}, flow={f:.1f}, quality={q:.2f}")
    else:
        print(f"Unit {i}: not built")

/Users/jkitchin/Dropbox/projects/discopt/python/discopt/solver.py:315: RuntimeWarning: divide by zero encountered in matmul
  residual = np.asarray(A_eq) @ np.asarray(x_full) - np.asarray(b_eq)
/Users/jkitchin/Dropbox/projects/discopt/python/discopt/solver.py:315: RuntimeWarning: overflow encountered in matmul
  residual = np.asarray(A_eq) @ np.asarray(x_full) - np.asarray(b_eq)
/Users/jkitchin/Dropbox/projects/discopt/python/discopt/solver.py:315: RuntimeWarning: invalid value encountered in matmul
  residual = np.asarray(A_eq) @ np.asarray(x_full) - np.asarray(b_eq)

Status:    optimal
Profit:    4402.22
Nodes:     97
Time:      2.54s

Unit 0: BUILT, mode=high-throughput, flow=60.5, quality=0.50
Unit 1: BUILT, mode=high-purity, flow=40.0, quality=0.90
Unit 2: BUILT, mode=high-purity, flow=40.0, quality=0.90

5. Big-M vs Hull Reformulation#

GDP models must be reformulated into standard MINLP form before solving. The two main approaches are [Grossmann, 2021]:

Big-M reformulation (default)#

For a disjunct \(y_k = 1 \Rightarrow g(x) \leq 0\), big-M adds:

\[g(x) \leq M(1 - y_k)\]

where \(M\) is a sufficiently large constant. This is simple but the LP relaxation can be weak (the feasible region is large), leading to more Branch & Bound nodes [Sawaya and Grossmann, 2005].

Hull reformulation#

The hull reformulation [Lee and Grossmann, 2000] introduces disaggregated variables \(x^k\) for each disjunct \(k\):

\[x = \sum_k x^k, \quad y_k g\left(\frac{x^k}{y_k}\right) \leq 0, \quad x^k_L y_k \leq x^k \leq x^k_U y_k\]

This gives a tighter relaxation (fewer nodes) but introduces more variables and constraints. For linear disjuncts, the hull relaxation is the convex hull of the disjunctive set.

When to use which?#

Criterion

Big-M

Hull

Relaxation quality

Weaker

Tighter

Number of variables

Fewer

More (disaggregated copies)

Number of constraints

Fewer

More

B&B nodes

More

Fewer

Per-node cost

Lower

Higher

Best for

Small models, tight big-M values

Large models, weak big-M values

In discopt, use gdp_method="big-m" (default) or gdp_method="hull" when calling solve().

def build_disjunctive_model():
    """Build a model where hull relaxation should help."""
    m = dm.Model("bigm_vs_hull")
    x = m.continuous("x", lb=0, ub=10)
    y = m.continuous("y", lb=0, ub=10)

    m.minimize((x - 5) ** 2 + (y - 5) ** 2)

    # x must be in [0, 2] or [8, 10]
    m.either_or(
        [[x <= 2], [x >= 8]],
        name="x_range",
    )

    # y must be in [0, 3] or [7, 10]
    m.either_or(
        [[y <= 3], [y >= 7]],
        name="y_range",
    )

    return m


# Solve with big-M (default)
m_bigm = build_disjunctive_model()
result_bigm = m_bigm.solve()
print("=== Big-M ===")
print(f"Objective: {result_bigm.objective:.4f}")
print(f"x = {result_bigm.x['x']:.4f}, y = {result_bigm.x['y']:.4f}")
print(f"Nodes: {result_bigm.node_count}")
print()

# Both methods should find the same optimal solution
# The difference is in the number of nodes explored
print("Both reformulations find the same optimal solution.")
print("The hull reformulation typically explores fewer nodes")
print("due to its tighter LP relaxation.")

=== Big-M ===
Objective: 13.0000
x = 8.0000, y = 3.0000
Nodes: 7

Both reformulations find the same optimal solution.
The hull reformulation typically explores fewer nodes
due to its tighter LP relaxation.

6. Multiple Big-M and LOA Reformulations#

Beyond standard big-M and hull, discopt provides two additional reformulation methods:

Multiple big-M (gdp_method="mbigm")#

Standard big-M uses a single global bound \(M\) per constraint. Multiple big-M solves a small LP relaxation to compute tight, constraint-specific \(M\) values:

\[g(x) \le M_k(1 - y_k) \quad \text{where } M_k = \max_{x \in X_k} g(x)\]

This tightens the continuous relaxation at the cost of one LP solve per constraint:

result = m.solve(gdp_method="mbigm")  # LP-tightened big-M constants

Logic-based Outer Approximation (gdp_method="loa")#

The LOA method [Türkay and Grossmann, 1996] decomposes the GDP into alternating subproblems:

  1. MILP master problem: linear outer approximation of all nonlinear constraints, solved with HiGHS. Provides a valid lower bound and suggests integer assignments.

  2. NLP subproblem: with integer variables fixed to the master solution, solve the continuous NLP to find a feasible solution and generate OA cuts.

The algorithm iterates until the gap \(|UB - LB| / |UB|\) is below tolerance. LOA is often faster than B&B on problems where the NLP subproblem is easy but the combinatorial structure is complex.

Requires highspy for the MILP master:

pip install highspy

# Compare reformulation methods on a disjunctive problem
def build_disjunctive_mbigm():
    m = dm.Model("mbigm_demo")
    x = m.continuous("x", lb=0, ub=10)
    y = m.continuous("y", lb=0, ub=10)

    d1 = m.make_disjunct("low_x")
    d1.subject_to(x <= 3)
    d1.subject_to(y >= 2)

    d2 = m.make_disjunct("high_x")
    d2.subject_to(x >= 7)
    d2.subject_to(y <= 5)

    m.add_disjunction([d1, d2])
    m.minimize(x**2 + y**2)
    return m


print("Solving with big-M:")
r1 = build_disjunctive_mbigm().solve(gdp_method="big-m")
print(f"  Status: {r1.status}  Objective: {r1.objective:.4f}")

print("Solving with multiple big-M (LP-tightened):")
r2 = build_disjunctive_mbigm().solve(gdp_method="mbigm")
print(f"  Status: {r2.status}  Objective: {r2.objective:.4f}")

print("Solving with convex hull:")
r3 = build_disjunctive_mbigm().solve(gdp_method="hull")
print(f"  Status: {r3.status}  Objective: {r3.objective:.4f}")

# LOA requires highspy:
# print("Solving with LOA:")
# r4 = build_disjunctive_mbigm().solve(gdp_method="loa")
# print(f"  Status: {r4.status}  Objective: {r4.objective:.4f}")

Solving with big-M:

  Status: optimal  Objective: 4.0000
Solving with multiple big-M (LP-tightened):
  Status: optimal  Objective: 4.0000
Solving with convex hull:
  Status: optimal  Objective: 4.0000

All four methods should produce the same optimal objective value, though the number of nodes explored and solve time may differ. For problems where the big-M constants are very loose, mbigm or hull typically requires fewer B&B nodes. loa is effective when the linear relaxation is strong but the NLP subproblems are easy to solve.

Summary#

discopt’s GDP support lets you model discrete-continuous optimization problems naturally:

Method

Use case

m.either_or(disjuncts)

Exactly one group of constraints must hold

m.if_then(y, constraints)

Constraints active only when \(y = 1\)

m.implies(y1, y2)

\(y_1 \Rightarrow y_2\) (precedence)

m.iff(y1, y2)

\(y_1 \Leftrightarrow y_2\) (equivalence)

m.at_least(k, binaries)

At least \(k\) selected

m.at_most(k, binaries)

At most \(k\) selected

m.exactly(k, binaries)

Exactly \(k\) selected

These high-level constructs are automatically reformulated into MINLP form and solved via spatial Branch & Bound.

gdp_method=

Description

"big-m" (default)

Standard big-M relaxation

"mbigm"

LP-tightened per-constraint big-M constants

"hull"

Convex hull reformulation (tightest, more variables)

"loa"

Logic-based Outer Approximation (requires highspy)

References#

  • Raman and Grossmann [1994] – Modelling and computational techniques for logic-based integer programming.

  • Lee and Grossmann [2000] – New algorithms for nonlinear generalized disjunctive programming.

  • Trespalacios and Grossmann [2015] – Review of mixed-integer nonlinear and generalized disjunctive programming methods.

  • Sawaya and Grossmann [2005] – Reformulations, relaxations and cutting planes for GDP.

  • Grossmann [2021] – Advanced optimization textbook covering GDP theory.

  • Türkay and Grossmann [1996] – Logic-based Outer Approximation algorithm for MINLP.

  • Anderson et al. [2020] – Strong MILP formulations for trained neural networks.