2.6 Dual of the BIM production problem

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."

Derivation of the dual problem

In a previous notebook, we introduce the BIM production problem and showed that it can be modeled as the following LO problem:

\[\begin{split} \begin{align*} \max \quad & 12 x_1 + 9 x_2 \\ \text{s.t.} \quad & x_1 \leq 1000 &\text{(silicon)}\\ & x_2 \leq 1500 &\text{(germanium)}\\ & x_1 + x_2 \leq 1750 &\text{(plastic)}\\ & 4 x_1 + 2 x_2 \leq 4800 &\text{(copper)}\\ & x_1, x_2 \geq 0. \end{align*} \end{split}\]

In this notebook, we will derive step by step its dual problem.

One can construct bounds for the value of objective function of the original problem by multiplying the constraints by non-negative numbers and adding them to each other so that the left-hand side looks like the objective function, while the right-hand side is the corresponding bound.

Let \(\lambda_1,\lambda_2,\lambda_3,\lambda_4\) be non-negative numbers. If we multiply each of these variables by one of the four constraints of the original problem and sum all of them side by side to obtain the inequality

\[ \begin{align*} (\lambda_1+\lambda_3+4\lambda_4) x_1 + (\lambda_2+\lambda_3+2 \lambda_4) x_2 \leq 1000 \lambda_1 + 1500 \lambda_2 + 1750 \lambda_3 + 4800 \lambda_4. \end{align*} \]

It is clear that if \(\lambda_1,\lambda_2,\lambda_3,\lambda_4 \geq 0\) satisfy

\[\begin{split} \begin{align*} \lambda_1+\lambda_3+4\lambda_4 & \geq 12,\\ \lambda_2+\lambda_3+2 \lambda_4 & \geq 9, \end{align*} \end{split}\]

then we have the following:

\[\begin{split} \begin{align*} 12 x_1 + 9 x_2 & \leq (\lambda_1+\lambda_3+4\lambda_4) x_1 + (\lambda_2+\lambda_3+2 \lambda_4) x_2 \\ & \leq 1000 \lambda_1 + 1500 \lambda_2 + 1750 \lambda_3 + 4800 \lambda_4, \end{align*} \end{split}\]

where the first inequality follows from the fact that \(x_1, x_2 \geq 0\), and the most right-hand expression becomes an upper bound on the optimal value of the objective.

If we seek \(\lambda_1,\lambda_2,\lambda_3,\lambda_4 \geq 0\) such that the upper bound on the RHS is as tight as possible, that means that we need to minimize the expression \(1000 \lambda_1 + 1500 \lambda_2 + 1750 \lambda_3 + 4800 \lambda_4\). This can be formulated as the following LO, which we name the dual problem:

\[\begin{split} \begin{align*} \min \quad & 1000 \lambda_1 + 1500 \lambda_2 + 1750 \lambda_3 + 4800 \lambda_4 \\ \text{s.t.} \quad & \lambda_1+\lambda_3+4\lambda_4 \geq 12,\\ & \lambda_2+\lambda_3+2 \lambda_4 \geq 9,\\ & \lambda_1,\lambda_2,\lambda_3,\lambda_4 \geq 0. \end{align*} \end{split}\]

It is easy to solve and find the optimal solution \((\lambda_1,\lambda_2,\lambda_3,\lambda_4)=(0,0,6,1.5)\), for which the objective functions takes the value \(17700\). Such a value is (the tightest) upper bound for the original problem.

The Pyomo code that implements and solves the dual problem is given below.

model = pyo.ConcreteModel("BIM production planning dual problem")

# Decision variables and their domain
model.y1 = pyo.Var(domain=pyo.NonNegativeReals)
model.y2 = pyo.Var(domain=pyo.NonNegativeReals)
model.y3 = pyo.Var(domain=pyo.NonNegativeReals)
model.y4 = pyo.Var(domain=pyo.NonNegativeReals)

# Objective function
model.obj = pyo.Objective(
    sense=pyo.minimize,
    expr=1000 * model.y1 + 1500 * model.y2 + 1750 * model.y3 + 4800 * model.y4,
)

# Constraints
model.x1 = pyo.Constraint(expr=model.y1 + model.y3 + 4 * model.y4 >= 12)
model.x2 = pyo.Constraint(expr=model.y2 + model.y3 + 2 * model.y4 >= 9)

# Solve and print solution
SOLVER.solve(model)
print(
    f"y = ({model.y1.value:.2f}, {model.y2.value:.2f}, {model.y3.value:.2f}, {model.y4.value:.2f})"
)
print(f"optimal value = {pyo.value(model.obj):.2f}")

y = (0.00, 0.00, 6.00, 1.50)
optimal value = 17700.00

Note that since the original LO is feasible and bounded, strong duality holds and the optimal value of the primal problem coincides with the optimal value of the dual problem.

If we are interested only in the optimal value of the dual variables, we can solve the original problem and ask Pyomo to return us the optimal values of the dual variables.

model = pyo.ConcreteModel("BIM production planning with decorators")

# Decision variables and their domains
model.x1 = pyo.Var(domain=pyo.NonNegativeReals)
model.x2 = pyo.Var(domain=pyo.NonNegativeReals)


# Objective function defined using a decorator
@model.Objective(sense=pyo.maximize)
def profit(m):
    return 12 * m.x1 + 9 * m.x2


# Constraints defined using decorators
@model.Constraint()
def silicon(m):
    return m.x1 <= 1000


@model.Constraint()
def germanium(m):
    return m.x2 <= 1500


@model.Constraint()
def plastic(m):
    return m.x1 + m.x2 <= 1750


@model.Constraint()
def copper(m):
    return 4 * m.x1 + 2 * m.x2 <= 4800


model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
SOLVER.solve(model)

for i, c in enumerate(model.component_objects(pyo.Constraint)):
    print(f"Constraint name: {c}")
    print(
        f"Optimal value corresponding dual variable: y_{i+1} = {model.dual[c]:0.2f}\n"
    )

Constraint name: silicon
Optimal value corresponding dual variable: y_1 = -0.00

Constraint name: germanium
Optimal value corresponding dual variable: y_2 = -0.00

Constraint name: plastic
Optimal value corresponding dual variable: y_3 = 6.00

Constraint name: copper
Optimal value corresponding dual variable: y_4 = 1.50