04-simple-production-model-gdp.ipynb

Kernel: Python 3 (ipykernel)

3.4 Production model using disjunctions

Preamble: Install Pyomo and a solver

The following cell sets and verifies a global SOLVER for the notebook. If run on Google Colab, the cell installs Pyomo and the HiGHS solver, while, if run elsewhere, it assumes Pyomo and HiGHS have been previously installed. It then sets to use HiGHS as solver via the appsi module and a test is performed to verify that it is available. The solver interface is stored in a global object SOLVER for later use.

import sys
 
if 'google.colab' in sys.modules:
    %pip install pyomo >/dev/null 2>/dev/null
    %pip install highspy >/dev/null 2>/dev/null
 
solver = 'appsi_highs'

import pyomo.environ as pyo
SOLVER = pyo.SolverFactory(solver)

assert SOLVER.available(), f"Solver {solver} is not available."

Disjunctions

Disjunctions appear in applications where there is choice among discrete alternatives. Given two logical propositions $\alpha$ and $\beta$, the "or" disjunction is denoted by $\vee$ and defined by the truth table

$\alpha$	$\beta$	$\alpha \vee \beta$
False	False	False
True	False	True
False	True	True
True	True	True

The "exclusive or" is denoted by $\veebar$ and defined by the truth table

$\alpha$	$\beta$	$\alpha \veebar \beta$
False	False	False
True	False	True
False	True	True
True	True	False

This notebook shows how to express disjunctions in Pyomo models using the Generalized Disjunctive Programming (GDP) extension for a simple production model.

Multi-product factory optimization

A small production facility produces two products, $X$ and $Y$. With current technology $\alpha$, the facility is subject to the following conditions and constraints:

• Product $X$ requires 1 hour of labor A, 2 hours of labor B, and 100$ of raw material. Product $X$ sells for 270$ per unit. The daily demand is limited to 40 units.

• Product $Y$ requires 1 hour of labor A, 1 hour of labor B, and 90$ of raw material. Product $Y$ sells for 210$ per unit with unlimited demand.

• There are 80 hours per day of labor A available at a cost of 50$/hour.

• There are 100 hours per day of labor B available at a cost of 40$/hour.

Using the given data we see that the net profit for each unit of $X$ and $Y$ is 40$ and 30$, respectively. The optimal product strategy is the solution to a linear optimization

$$\begin{align*} \max \quad & 40 x + 30 y\\ \text{s.t.} \quad & x \leq 40 & \text{(demand)}\\ & x + y \leq 80 & \text{(labor A)} \\ & 2 x + y \leq 100 & \text{(labor B)}\\ & x, y \geq 0. \end{align*}$$

m = pyo.ConcreteModel("Multi-Product Factory")

m.production_x = pyo.Var(domain=pyo.NonNegativeReals)
m.production_y = pyo.Var(domain=pyo.NonNegativeReals)


@m.Objective(sense=pyo.maximize)
def maximize_profit(m):
    return 40 * m.production_x + 30 * m.production_y


@m.Constraint()
def demand(m):
    return m.production_x <= 40


@m.Constraint()
def laborA(m):
    return m.production_x + m.production_y <= 80


@m.Constraint()
def laborB(m):
    return 2 * m.production_x + m.production_y <= 100


SOLVER.solve(m)

print(f"Profit = ${pyo.value(m.maximize_profit):.2f}")
print(f"Production X = {pyo.value(m.production_x)}")
print(f"Production Y = {pyo.value(m.production_y)}")

Profit = $2600.00
Production X = 20.0
Production Y = 60.0

Labor B is a relatively high cost for the production of product $X$. Suppose a new technology $\beta$ has been developed with the potential to cut costs by reducing the time required to finish product $X$ to $1.5$ hours, but requires more highly skilled labor with a unit cost of $60$$ per hour.

The net profit for unit of product $X$ with technology $\alpha$ is equal to $270 - 100 - 50 - 2 \cdot 40 = 40$$, while with technology $\beta$ is equal to $270 - 100 - 50 - 1.5 \cdot 40 = 60$$.

We need to assess whether the new technology is beneficial, that is, whether adopting it would lead to higher profits. The decision here is whether to use technology $\alpha$ or $\beta$.

In this situation we have an `either-or' structure for both the objective and for Labor B constraint:

$$ \underbrace{p = 40x + 30y, \ 2 x + y \leq 100}_{\text{$\alphaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 12: technology}̲} \quad \text{ …\betaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 12: technology}̲}. $

There are several commonly used techniques for embedding disjunctions into mixed-integer linear optimization problems, which we will explore in this notebook.

MILO implementation

The first approach is using the "big-M" technique introduces a single binary decision variable $z$ associated with choosing technology $\alpha$ ($z=0$) or technology $\beta$ ($z=1$). Using MILO, we can formulate this problem as follows:

$$ \begin{align*} \max \quad & \text{profit}\\ \text{s.t.} \quad & x \leq 40 & \text{(demand)}\\ & x + y \leq 80 & \text{(labor A)} \\ & \text{profit} \leq 40x + 30y + M z & \text{(profit with technology $\alphaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ \\ & \text…\betaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲\\ & 2 x + …\alphaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ \\ & 1.5 x…\betaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ \\ & z \in…$

where the variable $z \in \{ 0, 1\}$ "activates" the constraints related to the old or new technology, respectively, and $M$ is a large enough constant. It can be implemented in Pyomo as follows.

m = pyo.ConcreteModel("Multi-Product Factory - MILO formulation")

m.profit = pyo.Var(domain=pyo.NonNegativeReals)
m.production_x = pyo.Var(domain=pyo.NonNegativeReals)
m.production_y = pyo.Var(domain=pyo.NonNegativeReals)

m.z = pyo.Var(domain=pyo.Binary)
M = 10000


@m.Objective(sense=pyo.maximize)
def maximize_profit(m):
    return m.profit


@m.Constraint()
def profit_constr_1(m):
    return m.profit <= 40 * m.production_x + 30 * m.production_y + M * m.z


@m.Constraint()
def profit_constr_2(m):
    return m.profit <= 60 * m.production_x + 30 * m.production_y + M * (1 - m.z)


@m.Constraint()
def demand(m):
    return m.production_x <= 40


@m.Constraint()
def laborA(m):
    return m.production_x + m.production_y <= 80


@m.Constraint()
def laborB_1(m):
    return 2 * m.production_x + m.production_y <= 100 + M * m.z


@m.Constraint()
def laborB_2(m):
    return 1.5 * m.production_x + m.production_y <= 100 + M * (1 - m.z)


SOLVER.solve(m)

print(f"Profit = ${pyo.value(m.maximize_profit):.2f}")
print(f"Production X = {pyo.value(m.production_x)}")
print(f"Production Y = {pyo.value(m.production_y)}")

Profit = $3600.00
Production X = 40.0
Production Y = 40.0

Disjunctive programming implementation

Alternatively, we can formulate our problem using a disjunction, preserving the logical structure, as follows:

$$\begin{align*} \max \quad & \text{profit}\\ \text{s.t.} \quad & x \leq 40 & \text{(demand)}\\ & x + y \leq 80 & \text{(labor A)} \\ & \begin{bmatrix} \text{profit} = 40x + 30y\\ 2 x + y \leq 100 \end{bmatrix} \veebar \begin{bmatrix} \text{profit} = 60x + 30y\\ 1.5 x + y \leq 100 \end{bmatrix}\\ & x, y \geq 0. \end{align*}$$

This formulation, should the software be capable of handling it, has the benefit that the solver can intelligently partition the problem's solution into various sub-cases, based on the given disjunction. Pyomo natively supports disjunctions, as illustrated in the following implementation.

m = pyo.ConcreteModel("Multi-Product Factory - Disjunctive Programming")

m.profit = pyo.Var(bounds=(-1000, 10000))
m.x = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0, 1000))
m.y = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0, 1000))


@m.Objective(sense=pyo.maximize)
def maximize_profit(m):
    return m.profit


@m.Constraint()
def demand(m):
    return m.x <= 40


@m.Constraint()
def laborA(m):
    return m.x + m.y <= 80


# Define a disjunction using Pyomo's Disjunction component
# The 'xor=True' indicates that only one of the disjuncts must be true
@m.Disjunction(xor=True)
def technologies(m):
    # The function returns a list of two disjuncts
    # each containing a profit and a constraint
    return [
        [m.profit == 40 * m.x + 30 * m.y, 2 * m.x + m.y <= 100],
        [m.profit == 60 * m.x + 30 * m.y, 1.5 * m.x + m.y <= 100],
    ]


# Transform the Generalized Disjunctive Programming (GDP) model using
# the big-M method into a MILO problem and solve it
pyo.TransformationFactory("gdp.bigm").apply_to(m)
SOLVER.solve(m)

print(f"Profit = ${pyo.value(m.maximize_profit):.2f}")
print(f"Production X = {pyo.value(m.x)}")
print(f"Production Y = {pyo.value(m.y)}")

Profit = $3600.00
Production X = 40.0
Production Y = 40.0