BIM production#

# install dependencies and select solver
%pip install -q amplpy matplotlib pandas

SOLVER = "highs"

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

General LO formulation#

The simplest and most scalable class of optimization problems is the one where the objective function and the constraints are formulated using the simplest possible type of functions - linear functions. A linear optimization (LO) is a problem of the form

\[\begin{split} \begin{align*} \min \quad & c^\top x\\ \text{s.t.} \quad & A x \leq b\\ & x \geq 0, \nonumber \end{align*} \end{split}\]

where the \(n\) (decision) variables are grouped in a vector \(x \in \mathbb{R}^n\), \(c \in \mathbb{R}^n\) are the objective coefficients, and the \(m\) linear constraints are described by the matrix \(A \in \mathbb{R}^{m \times n}\) and the vector \(b \in \mathbb{R}^m\).

Of course, linear problems could also (i) be maximization problems, (ii) involve equality constraints and constraints of the form \(\geq\), and (iii) have unbounded or non-positive decision variables \(x_i\)ā€™s. In fact, any LP problem with such features can be easily converted to the ā€˜canonicalā€™ LP form by adding/removing variables and/or multiplying specific inequalities by \(-1\).

The microchip production problem#

The company BIM (Best International Machines) produces two types of microchips, logic chips (1 g silicon, 1 g plastic, 4 g copper) and memory chips (1 g germanium, 1 g plastic, 2 g copper). Each of the logic chips can be sold for a 12 ā‚¬ profit, and each of the memory chips for a 9 ā‚¬ profit. The current stock of raw materials is as follows: 1000 g silicon, 1500 g germanium, 1750 g plastic, 4800 g of copper. How many microchips of each type should be produced to maximize the profit while respecting the raw material stock availability?

Let \(x_1\) denote the number of logic chips and \(x_2\) that of memory chips. This decision can be reformulated as an optimization problem of the following form:

\[\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}\]

The problem has \(n=2\) decision variables and \(m=4\) constraints. Using the standard notation introduced above, denote the vector of decision variables by \(x = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\) and define the problem coefficients as

\[\begin{split} \begin{align*} c = \begin{pmatrix} 12 \\ 9 \end{pmatrix}, \qquad A = \begin{bmatrix} 1 & 0\\ 0 & 1\\ 1 & 1\\ 4 & 2\\ \end{bmatrix}, \quad \text{ and } \quad b = \begin{pmatrix} 1000 \\ 1500 \\ 1750 \\ 4800 \end{pmatrix}. \end{align*} \end{split}\]

This model can be implemented and solved in AMPL as follows.

%%ampl_eval

var x1 >= 0;
var x2 >= 0;

maximize profit: 12 * x1 + 9 * x2;

s.t. silicon: x1 <= 1000;
s.t. germanium: x2 <= 1500;
s.t. plastic: x1 + x2 <= 1750;
s.t. copper: 4 * x1 + 2 * x2 <= 4800;

ampl.option["solver"] = SOLVER
ampl.solve()

print(f'x = ({ampl.var["x1"].value():.1f}, {ampl.var["x2"].value():.1f})')
print(f'optimal value = {ampl.obj["profit"].value():.2f}')

HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 17700
2 simplex iterations
0 barrier iterations
x = (650.0, 1100.0)
optimal value = 17700.00

Dual problem#

One can construct bounds for the value of objective function 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{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*} \]

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 LP, 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. Here, we present the AMPL code for this example.

%%writefile bim_dual.mod

var y1 >= 0;
var y2 >= 0;
var y3 >= 0;
var y4 >= 0;

minimize obj: 1000 * y1 + 1500 * y2 + 1750 * y3 + 4800 * y4;

s.t. x1: y1 + y3 + 4 * y4 >= 12;
s.t. x2: y2 + y3 + 2 * y4 >= 9;

Writing bim_dual.mod

ampl = AMPL()
ampl.read("bim_dual.mod")
ampl.option["solver"] = SOLVER
ampl.solve()

print(
    f'y = ({ampl.var["y1"].value():.1f}, {ampl.var["y2"].value():.1f}, {ampl.var["y3"].value():.1f}, {ampl.var["y4"].value():.1f})'
)
print(f'optimal value = {ampl.obj["obj"].value():.2f}')

HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 17700
3 simplex iterations
0 barrier iterations
y = (0.0, 0.0, 6.0, 1.5)
optimal value = 17700.00

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