3.1 BIM production with perturbed data

3.1 BIM production with perturbed data#

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

Problem description#

The company BIM realizes that a \(1\%\) fraction of the copper always gets wasted while producing both types of microchips, more specifically \(1\%\) of the required amount. This means that it actually takes \(4.04\) gr of copper to produce a logic chip and \(2.02\) gr of copper to produce a memory chip. If we rewrite the linear problem of the basic BIM problem and modify accordingly the coefficients in the corresponding constraints, we obtain the following problem

\[\begin{split} \begin{align*} \max \quad & 12x_1+9x_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.04 x_1+2.02 x_2 \leq 4800 & \text{(copper with waste)}\\ &x_1, x_2 \geq 0. \end{align*} \end{split}\]

If we solve it, we obtain a different optimal solution than the original one, namely \((x_1,x_2) \approx (626.238,1123.762)\) and an optimal value of roughly \(17628.713\). The new optimal solution is not integral but, in fact, there is no constraint requiring \(x_1\) and \(x_2\) to be integers.

import pyomo.environ as pyo

m = pyo.ConcreteModel("BIM perturbed LO")

m.x1 = pyo.Var(domain=pyo.NonNegativeReals)
m.x2 = pyo.Var(domain=pyo.NonNegativeReals)

m.profit = pyo.Objective(expr=12 * m.x1 + 9 * m.x2, sense=pyo.maximize)

m.silicon = pyo.Constraint(expr=m.x1 <= 1000)
m.germanium = pyo.Constraint(expr=m.x2 <= 1500)
m.plastic = pyo.Constraint(expr=m.x1 + m.x2 <= 1750)
m.copper = pyo.Constraint(expr=4.04 * m.x1 + 2.02 * m.x2 <= 4800)

SOLVER.solve(m)

print(
    "x = ({:.3f}, {:.3f}), optimal value = {:.3f}".format(
        pyo.value(m.x1), pyo.value(m.x2), pyo.value(m.profit)
    )
)

x = (626.238, 1123.762), optimal value = 17628.713

In terms of production, we want to manufacture an integer number of microchips, but it is not clear how to implement the fractional optimal solution \((x_1,x_2) \approx (626.238,1123.762)\). Rounding down to \((x_1,x_2) = (626,1123)\) will intuitively yield a feasible solution, but this could lead to a minor loss of profit and/or an inefficient use of the available material. Rounding up to \((x_1,x_2) = (627,1124)\) could possibly lead to an unfeasible solution for which the available material is not enough. We could, of course, examine all the potential integer solutions by hand. However, if the problem had a more intricate structure or a greater number of decision variables, this would be much more difficult and may not lead to the true optimal solution.

A safer approach is to explicitly require the two decision variables to be non-negative integers, thus transforming the original LO problem into the following mixed-integer linear optimization (MILO) problem:

\[\begin{split} \begin{align*} \max \quad & 12x_1+9x_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.04 x_1+2.02 x_2 \leq 4800 & \text{(copper with waste)}\\ &x_1, x_2 \in \mathbb{N}. \end{align*} \end{split}\]

The optimal solution of this new MILO problem is \((x_1,x_2) = (626,1124)\) with a profit of \(17628\). Note that for this specific problem both the naive rounding strategies outlined above would not have yielded the true optimal solution. The Python code for obtaining the optimal solution using MILO solvers is given below.

m = pyo.ConcreteModel("BIM production planning with perturbed coefficients")

m.x1 = pyo.Var(domain=pyo.NonNegativeIntegers)
m.x2 = pyo.Var(domain=pyo.NonNegativeIntegers)

m.profit = pyo.Objective(expr=12 * m.x1 + 9 * m.x2, sense=pyo.maximize)

m.silicon = pyo.Constraint(expr=m.x1 <= 1000)
m.germanium = pyo.Constraint(expr=m.x2 <= 1500)
m.plastic = pyo.Constraint(expr=m.x1 + m.x2 <= 1750)
m.copper = pyo.Constraint(expr=4.04 * m.x1 + 2.02 * m.x2 <= 4800)

SOLVER.solve(m)

print(
    "x = ({:.3f}, {:.3f}), optimal value = {:.3f}".format(
        pyo.value(m.x1), pyo.value(m.x2), pyo.value(m.profit)
    )
)

x = (626.000, 1124.000), optimal value = 17628.000