Integer Programming (IP)
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
Table of Contents
1. Integer Programming (IP)¶
An integer programming problem is a mathematical optimization or feasibility program in which some or all of the variables are restricted to be integers.
$$ \begin{align*} \max \quad &y \\ \text{subject to} \quad &-x + y \leq 1 \\ &3x + 2y \leq 12 \\ & 2x + 3y \leq 12 \\ & x,y \geq 0 \\ & \color{red}{x,y \in \mathbb{Z}} \end{align*} $$
There are two main reasons for using integer variables when modeling problems as a linear program:
- The integer variables represent quantities that can only be integer. For example, it is not possible to build 3.7 cars.
- The integer variables represent decisions (e.g., whether to include an edge in a graph) and so should only take on the value 0 or 1.
These considerations occur frequently in practice and so integer linear programming can be used in many applications areas, some of which are briefly described below.
Manufacture of television sets
LP might give a production plan of 205.7 sets per week. In such a model, most people would have no trouble stating that production should be 205 sets per week (or even ''roughly 200 sets per week'').
Warehouses
Purchasing 0.7 warehouse at some location and 0.6 somewhere else would be of little value. However, warehouses come in integer quantities.
Integer programming is not a part of convex optimization problem. $\rightarrow$ we cannot use cvx.
IP Modeling Example
- Manufacturer produces two parts P and Q with machine A or B.
- Time of producing each unit of P
- in machine A: 50 min
- in machine B: 30 min
- Time of producing each unit of Q
- in machine A: 24 min
- in machine B: 33 min
- Working plan for a week
- 40 hrs of work on machine A
- 35 hrs of work on machine B
- The week starts with
- stock of 30 units of P
- stock of 90 units of Q
- demand of 75 units of P
- demand of 95 units of Q
Question: how to plan the production to end the week with the maximum stock?
Solution
Define decision variables
- $x$ = unit of P to be produced
- $y$ = unit of Q to be produced
Optimization Problem
$$\begin{align*} \max \quad &(30 + x - 75) + (90 + y - 95) \\\\ \text{subject to} \quad & 50x + 24y \leq 40 \times 60 \\ & 30x + 33y \leq 35 \times 60 \\ & x \geq 75 - 30 \\ & y \geq 95 - 90 \\ & \color{red}{x,y \in \mathbb{Z}} \end{align*}$$
2. Mixed Integer Linear Programming (MILP)¶
A mixed-integer linear program is a problem with
- Linear objective function, $c^T x$, where $c$ is a column vector of constants, and $x$ is the column vector of unknowns
- Bounds and linear constraints, but no nonlinear constraints
- Restrictions on some components of $x$ to have integer values
milp in Python
$$ \begin{align*} \min \quad &c^Tx\\ \text{subject to} \quad & b_l \leq Ax \leq b_u\\ & l \leq x \leq u \\ & x_i \in \mathbb{Z}, \; i \in X_i \end{align*} $$
where $X_i$ is the set of indices of decision variables that must be integral
By default, $l = 0$ and $u=$ np.inf unless specified with bounds.
Integrality indicates the type of integrality constraint on each decision variable.
0: Continuous variable, no integrality constraint.
1: Integer variable, decision variable must be an integer within bounds.
3. Examples of MILP¶
Example 1
$$\begin{align*} \min \quad & \begin{bmatrix} 0 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\\\ \text{subject to} \quad & \begin{bmatrix} -1 & 1 \\ 3 & 2 \\ 2 & 3\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 12 \\ 12 \end{bmatrix} \\\\ & \begin{bmatrix} 0 \\ 0 \end{bmatrix} \leq \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\\\ & \color{red}{x_1, x_2 \in \mathbb{Z} \; \text{(integer)}} \end{align*} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = np.array([0, -1])
A = np.array([[-1, 1],
[3, 2],
[2, 3]])
b_l = np.array([-np.inf, -np.inf, -np.inf])
b_u = np.array([1, 12, 12])
constraints = LinearConstraint(A, b_l, b_u)
bounds = Bounds(0, )
integrality = np.array([1, 1])
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
If the integer constraints are removed, the original problem becomes a linear programming problem.
$$\begin{align*} \min \quad & \begin{bmatrix} 0 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\\\ \text{subject to} \quad & \begin{bmatrix} -1 & 1 \\ 3 & 2 \\ 2 & 3\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 12 \\ 12 \end{bmatrix} \\\\ & \begin{bmatrix} 0 \\ 0 \end{bmatrix} \leq \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\\\ \end{align*} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = -np.array([0, 1])
A = np.array([[-1, 1],
[3, 2],
[2, 3]])
b_l = np.array([-np.inf, -np.inf, -np.inf])
b_u = np.array([1, 12, 12])
constraints = LinearConstraint(A, b_l, b_u)
bounds = Bounds(0, )
integrality = np.array([0, 0])
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
The optimal solutions of the integer problem are the points $(1, 2)$ and $(2, 2)$ which both have an objective value of 2. The unique optimum of the relaxation is $(1.8, 2.8)$ with objective value of 2.8. If the solution of the relaxation is rounded to the nearest integers, it is not feasible for the ILP.
Example 2
Carefully observe how equality constraints are managed in milp
$$\begin{align*} \min_{x} \quad & \begin{bmatrix} -3 & -2 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \\ \text{subject to} \quad & \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \leq 7 \\ & \begin{bmatrix} 4 & 2 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = 12 \\ & 0 \leq \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \\ & x_3 \in \{0,1\} \end{align*} $$
c = np.array([-3, -2, -1])
A = np.array([[1, 1, 1],
[4, 2, 1]])
b_l = np.array([-np.inf, 12])
b_u = np.array([7, 12])
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.array([0, 0, 1])
bounds = Bounds(lb = [0, 0, 0],
ub = [np.inf, np.inf, 1])
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(res.fun)
Example 3
Back to the first manufacturing problem
$$\begin{align*} \max \quad &(30 + x - 75) + (90 + y - 95) \\\\ \text{subject to} \quad & 50x + 24y \leq 40 \times 60 \\ & 30x + 33y \leq 35 \times 60 \\ & x \geq 75 - 30 \\ & y \geq 95 - 90 \\ & \color{red}{x,y \in \mathbb{Z}} \end{align*}$$
import numpy as np
from scipy.optimize import linprog
c = [-1, -1]
A = [[50, 24], [30, 33]]
integrality = np.array([1, 1])
b = [2400, 2100]
x_bounds = (45, None)
y_bounds = (5, None)
result = linprog(c,
A_ub = A,
b_ub = b,
integrality = integrality,
bounds = [x_bounds, y_bounds])
x_opt, y_opt = result.x
# Calculate the final stock of P and Q
final_stock_P = x_opt - 45
final_stock_Q = y_opt - 5
total_stock = final_stock_P + final_stock_Q
print(f"Optimal units of P to produce: {x_opt:.2f}")
print(f"Optimal units of Q to produce: {y_opt:.2f}")
print(f"Final stock of P: {final_stock_P:.2f}")
print(f"Final stock of Q: {final_stock_Q:.2f}")
print(f"Total final stock: {total_stock:.2f}")
4. Binary Integer Programming¶
Binary integer programming is a specific subset of integer programming, or more generally, mixed-integer linear programming.
In binary integer programming, each variable is restricted to take on only the values 0 or 1. These binary variables are often used to represent discrete decisions, such as the selection or rejection of an option, the activation or deactivation of a switch, or a binary response (e.g., yes/no).
Example 1
$$ \begin{align*} \min \quad & x_1 + x_2 + x_3 + x_4 + x_5 + x_6 \\ \\ \text{subject to} \quad & x_1 + x_2 \geq 1 \\ & x_1 + x_2 + x_6\ \geq 1 \\ & x_3 + x_4 \geq 1 \\ & x_3 + x_4 + x_5 \geq 1 \\ & x_4 + x_5 + x_6 \geq 1 \\ & x_5 + x_6 \geq 1 \\ & \color{red}{x_1 \ x_2 \ x_3 \ x_4 \ x_5 \ x_6 \in \{0, 1\}} \\ \end{align*}$$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = np.array([1, 1, 1, 1, 1, 1])
A = np.array([[1, 1, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 1],
[0, 0, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0],
[0, 0, 0, 1, 1, 1],
[0, 0, 0, 0, 1, 1]])
b_l = np.array([1, 1, 1, 1, 1, 1])
b_u = np.array([np.inf, np.inf, np.inf, np.inf, np.inf, np.inf])
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.array([1, 1, 1, 1, 1, 1])
bounds = Bounds(0, 1) # 0 <= x_i <= 1
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(res.fun)
Example 2
To minimize the function
$$
f(x) = -9x_1 - 5x_2 - 6x_3 - 4x_4
$$
subject to the constraints
$$ \begin{bmatrix} 6 & 3 & 5 & 2 \\ 0 & 0 & 1 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} \leq \begin{bmatrix} 9 \\ 1 \\ 0 \\ 0 \\ \end{bmatrix}, $$
where $x_1, x_2, x_3$ and $x_4$ are binary integers.
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
# Define the coefficients of the objective function
f = np.array([-9, -5, -6, -4])
# Define the coefficients of the constraints
A = np.array([[6, 3, 5, 2],
[0, 0, 1, 1],
[-1, 0, 1, 0],
[0, -1, 0, 1]])
# Define the lower and upper bounds of the constraints
b_l = np.array([-np.inf, -np.inf, -np.inf, -np.inf])
b_u = np.array([9, 1, 0, 0])
constraints = LinearConstraint(A, b_l, b_u)
# Define binary integers for the variables using the bounds and integrality
integrality = np.array([1, 1, 1, 1])
bounds = Bounds(0, 1) # 0 <= x_i <= 1
# Solve the integer linear programming problem
res = milp(c = f,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
4.1. Knapsack Problems¶
There is a container (knapsack) of capacity $C = 20$. Furthermore, there is a set 6 of objects. Each object has a weight and a value.
Determine the number of each item to include in a collection so that the total size is less than or equal to 20 and the total value is as large as possible.
|items|1|2|3|4|5|6 |:-:|:-:|:-:|:-:|:-:|:-:|:-: |weights|10|9|4|2|1|20 |values|175|90|20|50|10|200
Question: Can we formulate this knapsack problem into a binary integer program?
$\implies$ decision variables $\quad x_{i} \in \{0,1\} \quad i = 1,\cdots,6$
$$\begin{align*} \max_{x} \quad & 175x_1 + 90x_2 + 20x_3 + 50x_4 + 10x_5 + 200x_6 \\\\ \text{subject to} \quad & 10x_1 + 9x_2 + 4x_3 + 2x_4 + x_5 + 20x_6 \leq 20 \\ \quad & x_{i} \in \{0,1\} \quad i = 1,\cdots,6 \end{align*} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
values = np.array([175, 90, 20, 50, 10, 200])
weights = np.array([10, 9, 4, 2, 1, 20])
c = -values
A = weights
b_l = 0
b_u = 20
constraints = LinearConstraint(A, b_l, b_u)
# integrality = np.array([1, 1, 1, 1, 1, 1])
# integrality = np.array([1]*6)
integrality = np.full_like(values, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = -np.array([175, 90, 20, 50, 10, 200])
A = np.array([10, 9, 4, 2, 1, 20])
b_l = 0
b_u = 20
constraints = LinearConstraint(A, b_l, b_u)
# integrality = np.array([1, 1, 1, 1, 1, 1])
# integrality = np.array([1]*6)
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
4.2. Factory and Warehouse¶
A campany wants to
- build new factory in Ulsan, Busan, or both
- build new warehouse (only one has to be built)
- warehouse must be built in a city of a new factory
- Available capital to build is 10.
We want to maximize total profit.
| | |Profit|Cost to build |:-:|:-|:-:|:-: |1|Build a factory in Ulsan|7|5 |2|Build a factory in Busan|5|3 |3|Build a warehouse in Ulsan|6|5 |4|Build a warehouse in Busan|4|2
Form a binary integer programming to find the optimal locations of new factory and warehouse and solve it using milp in Python.
$$ \begin{array}{Icl} \max \quad & 7x_1 + 5x_2 + 6x_3 + 4x_4 \\\\ \text{subject to} \quad & 5x_1 + 3x_2 + 5x_3 + 2 x_4 \leq 10\\ & x_1 + x_2 \geq 1 \\ &x_3 + x_4 = 1 \\ &x_3 \leq x_1\\ &x_4 \leq x_2\\ &x_1, x_2, x_3, x_4 \in \{0, 1\} \end{array} \quad \implies \quad \begin{array} \quad \min \quad & -\begin{bmatrix} 7 & 5 & 6 & 4 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} \\\\ \text{subject to} \quad & \begin{bmatrix} 5 & 3 & 5 & 2\\-1 & -1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} \leq \begin{bmatrix} 10\\ -1 \\ 0 \\ 0 \end{bmatrix}\\ &\begin{bmatrix} 0 & 0 & 1 & 1\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = 1 \\ &x_1, x_2, x_3, x_4 \in \{0, 1\} \end{array} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = -np.array([7, 5, 6, 4])
A_ub = np.array([[5, 3, 5, 1],
[-1, -1, 0, 0],
[-1, 0, 1, 0],
[0, -1, 0, -1]])
b_ub = np.array([10, -1, 0, 0])
A_eq = np.array([0, 0, 1, 1])
b_eq = 1
constraints = [LinearConstraint(A = A_ub, lb = -np.inf, ub = b_ub),
LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)]
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = -np.array([7, 5, 6, 4])
A = np.array([[5, 3, 5, 1],
[-1, -1, 0, 0],
[-1, 0, 1, 0],
[0, -1, 0, -1],
[0, 0, 1, 1]])
b_l = np.array([-np.inf, -np.inf, -np.inf, -np.inf, 1])
b_u = np.array([10, -1, 0, 0, 1])
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(-res.fun)
4.3. Shortest Path¶
Binary integer programming can be applied to finding the shortest path in the graph.
Define
$$ x_i = \begin{cases} 0 & \text{if not part of shortest path}\\ 1 & \text{if part of shortest path}\end{cases} $$
Objective function
$$\min \quad x_1 + 4x_2 + 2x_3 + 4x_4 + 8x_5 + 11x_6 + 3x_7$$
- For each node
$$ \begin{align*} \sum x_{\text{in}} &= \sum x_{\text{out}} \\\\ x_{\text{in}}&: \text{entering edge to node } i \\ x_{\text{out}}&: \text{leaving edge from node } i \end{align*} $$
- At node 0, it has to transit to either node 1 or node 2
$$x_1 + x_2 = 1$$
- At node 4, it has to come from either node 5, 6, or 7
$$x_5 + x_6 + x_7 = 1$$
- At node 1
$$x_1 = x_3 + x_6$$
- At node 2
$$x_2 + x_3 = x_4 + x_5$$
- At node 3
$$x_4 = x_7$$
In the form of optimization
$$ \begin{align*} \min \quad & x_1 + 4x_2 + 2x_3 + 4x_4 + 8x_5 + 11x_6 + 3x_7\\\\ \text{subject to} \quad & \begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0\\ 1 & 0 & -1 & 0 & 0 & -1 & 0\\0 & 1 & 1 & -1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & -1\\ 0 & 0 & 0 & 0 & 1 & 1 & 1\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7\end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 1\end{bmatrix}\\ &x_1, x_2, x_3, x_4, x_5, x_6, x_7 \in \{0, 1\} \end{align*} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = np.array([1, 4, 2, 4, 8, 11, 3])
A_eq = np.array([[1, 1, 0, 0, 0, 0, 0],
[1, 0, -1, 0, 0, -1, 0],
[0, 1, 1, -1, -1, 0, 0],
[0, 0, 0, 1, 0, 0, -1],
[0, 0, 0, 0, 1, 1, 1]])
b_eq = np.array([1, 0, 0, 0, 1])
constraints = LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(res.fun)
4.4. Vertex Cover Problem¶
In graph theory, a vertex cover (sometimes node cover) of a graph is a set of vertices that includes at least one endpoint of every edge of the graph.
In computer science, the problem of finding a minimum vertex cover is a classical optimization problem
Example
Suppose you want to have an ATM in every street (block, district) of a city of Ulsan. However, each additional ATM machine costs more money. Therefore, we want to have the minimum number of ATM machines at the right location while they are accessible from every street.
- Decision variable
$$x_i = \begin{cases} 1 & \text{ATM machine allocated at the location $i$}\\ 0 & \text{no ATM machine at the location $i$}\end{cases}$$
$$ \begin{align*} \min \quad &x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7 + x_8\\\\ \text{subject to} \quad &x_1 + x_3 \geq 1\\ \quad &x_1 + x_4 \geq 1\\ \quad &x_1 + x_5 \geq 1\\ \quad &x_1 + x_6 \geq 1\\ \quad &x_1 + x_7 \geq 1\\ \quad &x_2 + x_4 \geq 1\\ \quad &x_2 + x_5 \geq 1\\ \quad &x_3 + x_8 \geq 1\\ \quad &x_5 + x_7 \geq 1 \end{align*} $$
c = np.array([1, 1, 1, 1, 1, 1, 1, 1])
A = np.array([[1, 0, 1, 0, 0, 0, 0, 0],
[1, 0, 0, 1, 0, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0],
[1, 0, 0, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 0, 0, 1, 0],
[0, 1, 0, 1, 0, 0, 0, 0],
[0, 1, 0, 0, 1, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 1, 0, 1, 0]])
b_l = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1])
b_u = np.array([np.inf]*9)
constraints = LinearConstraint(A = A, lb = b_l, ub = b_u)
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
print(res.fun)
Not a unique solution
$\{1, 2, 3, 5\}, \{1, 2, 3, 7\}, \{1, 2, 7, 8\}, \cdots$
4.5. Coloring Problems¶
- Coloring problem: color the vertex of a graph in such a way that adjacent vertices are differently colored, while minimizing the number of colors.
Example
Find the minimum number of colors that can cover the below Ulsan map without conflicting the adjacent rule.
- Graph representation
$$ \begin{align*} V &= \{ V_1, V_2, V_3, V_4, V_5\} \\ E &= \{E_{12}, E_{13}, E_{14}, E_{23}, E_{24}, E_{34}, E_{35}, E_{45} \} \end{align*} $$
- Define ILP decision variables
$$\begin{align*} y_k &= \begin{cases} 1 \quad \text{if color $k$ is used} \\ 0 \quad \text{if color $k$ is not used} \end{cases} \\\\ x_{ik} &= \begin{cases} 1 \quad \text{if vertex $i$ gets color $k$} \\ 0 \quad \text{otherwise} \end{cases} \end{align*} $$
- the coloring problem can be formulated as follows:
$$\begin{align*} \min \quad & \sum \limits_{k=1}^5 y_k \\ \text{subject to} \quad & \sum \limits_{k=1}^5 x_{ik} = 1 \\ & x_{ik} + x_{jk} \leq 1 \quad \; \text{for } E_{ij} \in E \text{ and } 1 \leq k \leq 5 \\ & x_{ik} \leq y_k \quad \quad \quad \text{for } 1 \leq i \leq 5 \text{ and } 1\leq k \leq 5 \\ & x_{ik}, y_{k} \in \{0,1\} \end{align*} $$
## coloring problem using binary integer programming
# Actually number of colours is at most 5
# Therefore we restrict our number of colours to be 5
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
n_vertices = 5
n_colors = 5
edges = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3), (2, 4), (3, 4)]
c = np.array([0] * 25 + [1] * 5)
A_ub = []
b_ub = []
for u, v in edges:
for k in range(n_colors):
row = np.zeros(n_vertices * n_colors + n_colors)
row[u * n_colors + k] = 1
row[v * n_colors + k] = 1
A_ub.append(row)
b_ub.append(1)
for i in range(n_vertices):
for k in range(n_colors):
row = np.zeros(n_vertices * n_colors + n_colors)
row[i * n_colors + k] = 1
row[n_vertices * n_colors + k] = -1
A_ub.append(row)
b_ub.append(0)
A_eq = np.zeros((n_vertices, n_vertices * n_colors + n_colors))
b_eq = np.ones(n_vertices)
for i in range(n_vertices):
A_eq[i, i * n_colors: (i + 1) * n_colors] = 1
A_ub = np.array(A_ub)
b_ub = np.array(b_ub)
A_eq = np.array(A_eq)
b_eq = np.array(b_eq)
integrality = np.array([1] * (30))
constraints = [LinearConstraint(A = A_ub, ub = b_ub),
LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)]
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
# Print the results
print("Optimal value:", res.fun)
# print("Optimal solution:", res.x)
x = res.x[:25]
X = np.array(x).reshape(5, 5)
y = res.x[25:]
print("\nX:", X)
print("\ny:", y)
The minimum number of colors are four.
4.6. Flat Pattern Design of 3D Shapes¶
A box with the dimensions specified is to be made by bending and welding a flat pattern cut out of sheet metal.
Design the at pattern with minimum welding length that can be made to the box.
Origami
- Japanese art of paper folding
- In modern usage, the word "origami" is often used as an inclusive term for all folding practices.
- The goal is to transform a flat square sheet of paper into a finished sculpture through folding and sculpting techniques.
- Graph Representation: Example 1
- Graph Representation: Example 2
- Modeling
$$ x_i = \begin{cases} 1 \quad \text{if two nodes are connected (welded).}\\ 0 \\ \end{cases} $$
Constraints
No disconnected node.
Only 6-1 edges are required. (no cycle)
$$\begin{align*} \max \quad & x_1 + x_2 + x_3 + x_4 + 2x_5 + 2x_6 + 2x_7 + 2x_8 + 3x_9 + 3x_{10} + 3x_{11} + 3x_{12} \\ \text{subject to} \quad & x_1 + x_2 + x_5 + x_8 \geq 1 \\ & x_2 + x_4 + x_{10} + x_{12} \geq 1 \\ & x_3 + x_4 + x_6 + x_7 \geq 1 \\ & x_1 + x_3 + x_9 + x_{11} \geq 1 \\ & x_5 + x_6 + x_9 + x_{10} \geq 1 \\ & x_7 + x_8 + x_{11} + x_{12} \geq 1 \\ & x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7 + x_8 + x_9 + x_{10} + x_{11} + x_{12} = 5 \end{align*} $$
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = -np.array([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3])
A_ub = np.array([[1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1],
[0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0],
[1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0],
[0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]])
b_lb = np.array([1, 1, 1, 1, 1, 1])
A_eq = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
b_eq = 5
constraints = [LinearConstraint(A = A_ub, lb = b_lb, ub = np.inf),
LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)]
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
- One possible solution
- Note that this is not a unique solution
4.7. Traveling Salesman Problem¶
Starting from his home, a salesman wishes to visit each of $(n-1)$ other cities and return home at minimal cost.
He must visit each city exactly once and it cost $c_{ij}$ to travel from city $i$ to city $j$.
What route should he select?
Starting from node 1, visit 2,3,4 once and return node 1 at minimal cost.
$$ x_{ij} = \begin{cases} 1 \quad \text{if travel from $i$ to $j$} \\ 0 \quad \text{otherwise} \end{cases} $$
- Optimization
$$\begin{align*} \min_{x} \quad & \sum \limits_{i=1}^4 \sum \limits_{j=1}^4 c_{ij} x_{ij} \\ \text{subject to} \quad & \sum \limits_{j} x_{ij} = 1 \quad \text{for } i=1,2,3,4 \\ & \sum \limits_{i} x_{ij} = 1 \quad \text{for } j=1,2,3,4 \\ & \sum \limits_{i \in S, j \in \overline{S}} x_{ij} \geq 1 \quad \text{for all prtitions} \; \left(S,\overline{S}\right) \quad 2 \leq \vert S \vert \leq m-2 \\ \end{align*} $$
$$\begin{align*} \\ \min_{x} \quad & 5x_{12} + 10x_{13} + 13x_{14} + 5x_{21} + 8x_{23} + 9x_{24} \\ & + 10x_{31} + 8x_{32} + 7x_{34} + 13x_{41} + 9x_{42} + 7x_{43} \\ \\ \text{subject to} \quad & x_{11} + x_{12} + x_{13} + x_{14} = 1 \\ & x_{21} + x_{22} + x_{23} + x_{24} = 1 \\ & x_{31} + x_{32} + x_{33} + x_{34} = 1 \\ & x_{41} + x_{42} + x_{43} + x_{44} = 1 \\ & x_{11} + x_{21} + x_{31} + x_{41} = 1 \\ & x_{12} + x_{22} + x_{32} + x_{42} = 1 \\ & x_{13} + x_{23} + x_{33} + x_{43} = 1 \\ & x_{14} + x_{24} + x_{34} + x_{44} = 1 \\ & x_{11} = x_{22} = x_{33} = x_{44} = 0 \\ \\ & \text{(no sub-closed loop conditions)} \\ \\ & x_{23} + x_{24} + x_{13} + x_{14} \geq 1 \\ & x_{31} + x_{32} + x_{41} + x_{42} \geq 1 \\ \\ & x_{21} + x_{24} + x_{31} + x_{34} \geq 1 \\ & x_{12} + x_{13} + x_{42} + x_{43} \geq 1 \\ \\ & x_{21} + x_{23} + x_{41} + x_{43} \geq 1 \\ & x_{12} + x_{14} + x_{32} + x_{34} \geq 1 \\ \end{align*} $$
- No sub-closed loop conditions
import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
from scipy.optimize import Bounds
c = np.array([5, 10, 13, 5, 8, 9, 10, 8, 7, 13, 9, 7])
A_eq = np.array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1],
[0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
[0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0]])
b_eq = np.array([1, 1, 1, 1, 1, 1, 1, 1])
constraints = LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)
integrality = np.full_like(c, True)
bounds = Bounds(0, 1)
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
# print(-res.fun)
What is wrong?
A_ub = np.array([[0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0]])
b_lb = np.array([1, 1])
constraints = [LinearConstraint(A = A_ub, lb = b_lb, ub = np.inf),
LinearConstraint(A = A_eq, lb = b_eq, ub = b_eq)]
res = milp(c = c,
constraints = constraints,
integrality = integrality,
bounds = bounds)
print(res.x)
- Possible solution
5. Any Ideas to Solve IP (Optional)¶
The bleow has not been completed yet
The naive way to solve an ILP is
- to simply remove the constraint that x is integer,
- solve the corresponding LP (called the LP relaxation of the ILP), and
- then round the entries of the solution to the LP relaxation.
But, not only may this solution not be optimal, it may not even be feasible
Rounding Does Not Give Any Useful Result
5.1. Relaxation¶
Original problem:
$$\begin{align*}
\text{minimize} \quad & f(x) \\
\text{subject to} \quad & x \in F
\end{align*}$$
Relaxed problem:
$$\begin{align*} \text{minimize} \quad & f(x) \\ \text{subject to} \quad & x \in F_R \\ \quad &\text{where} \; F \subseteq F_R \end{align*} $$
Theorem : Let $x^*$ and $x_R^*$ be optimal solutions of the original and relaxed problems, respectively. Then, $f(x_R^*) \leq f(x^*).$ If $x_R^* \in F, x_R^* $ is an optimal solution of the original problem.
Proof: (First part) Suppose $f(x_R^*) \leq f(x^*)$. Since $x^* \in F$ and $F \subseteq F_R,\; x^* \in F_R.$ This contradicts the fact that $x_R^*$ is an optimal solution of the relaxed problem. Therefore, $f(x_R^*) \leq f(x^*)$.
(Second part) Now we know $f(x_R^*) \leq f(x^*).$ If $x_R^* \in F$, therefore, $x_R^*$ is an optimal solution of the original problem.
5.2. ILP $\rightarrow$ LP Relaxation¶
Original problem (ILP):
$$\begin{align} \text{minimize}\quad & c^Tx \\ \text{subject to} \quad &Ax = b \\ & x \geq 0 \\ & x \in \mathbb{Z}^n \end{align}$$
Relaxed problem (LP):
$$\begin{align} \text{minimize}\quad & c^Tx \\ \text{subject to} \quad &Ax = b \\ & x \geq 0 \\ & x \in \mathbb{R}^n \end{align}$$
$Z^n \subseteq R^n \implies$ the above theorem applies.
For minimization problem, $c^Tx_R^*$ gives a lower bound to the optimal function value of the original ILP problem.
If $x_R^* \in Z^n,$ we have solved the ILP problem.
In general, round$\left(x_R^* \right) \neq x^* !$
5.3. Cutting Plane Algorithms¶
General case
$d^Tx \leq e$ is a valid inequality for $X \subseteq \mathbb{R}^n$ if $\forall x \in X, \; d^Tx \leq e.$
If $d^Tx \leq e$ is not a valid inequality for $X$ (i.e., invalid for $X$), it "cuts" $X$.
Basic idea: solve the relaxed problem. If $x_R^* \notin F$, add an inequality valid for $F$ but invalid for $F_R$ to the constraints and solve the relaxed problem again. Repeat until $x_R^* \in F$.
Need a way to automatically generate the inequalities that "cuts" well and efficiently solve the relaxed problems.
5.4. Ralph E. Gomory¶
- "Outline of an Algorithm for Integer Solutions to Linear Programs," Bulletin of the American Mathematical Society, Vol. 64, pp. 275-278, 1958.
5.5. Big Picture of Cutting Plane Algorithm¶
Want to solve a pure IP problem:
- Don`t know how to solve IP problems.
- Know how to solbe continuous problems (Simplex).
Outline of a solution procedure:
- Solve a continuous relaxation.
- Contains all originally feasible solutions, plus others.
- If optimal solution is integral, we are done.
- Otherwise, tighten the relaxation and repeat.
Continuous relaxation: $x_j \in \mathbb{Z} \rightarrow x_j \in \mathbb{R}$
Tightening:
- Add cutting planes (cut off current optimum).
5.6. Cutting Plane Algorithm¶
A cutting plane algorithm to solve pure integer programming problems works as follows.
- Solve the IP problem with continuous variables instead of discrete ones.
- If the resulting optimal solution $x^*$ is integral, stop $\implies$ optimal solution found.
- Generate a cut, i.e., a constraint which is satisfied by all feasible integer solutions but not by $x^*$.
- Add this new constraint, resolve problem, and go back to (2).
Terminates after finite number of iterations in (2). Resulting $x^*$ is integral and optimal
5.7. Example¶
Consider the following problem:
$$\begin{align*} \max \quad & 5x_1 + 8x_2\\ \text{subject to} \quad &x_1 + x_2 \leq 6 \\ &5x_1 + 9x_2 \leq 45 \\ &x_1, x_2 \geq 0 \\ &x_1, x_2 \in \mathbb{Z}. \end{align*}$$
Step 1. Solve the IP problem with continous variables instead of discrete ones.
Step 2. If the resulting optimal solution $x^*$ is integral, stop $\implies$ optimal solution found.
Resulting solution is $x^*$ = (2.25, 3.75) and hence not integral.
Step 3. Generate a cut, i.e., a constraint which is satisfied by all feasible integer solutions but not by $x^*$.
We generate cut $2x_1 + 3x_2 \leq 15.$
Step 4. Add this new constraint, resolve problem, and go back to (2).
New optimal solution is $x^*$ = (3,3).
Note: previous $x^*$ is not feasible anymore.
Step 2. If the resulting optimal solution $x^*$ is integral, stop $\Rightarrow$ optimal solution found.
$x^*$ = (3,3) is integral $\Rightarrow$ optimal solution found.
Remark. The cut we introduced only removed non-integral solutions. Cuts never cut off feasible solutions of the original IP problem!
5.8. Gomory cut¶
Assume $x_1,...,x_n \geq 0$ and integral. We show how to construct a Gomory Cut for
$$ a_1x_1 + \cdots + a_nx_n = b,$$
where $a_j,b \in \mathbb{R}$ (not necessarily integral). Note that this can be written as
$$(\lfloor a_1 \rfloor + \underbrace{[a_1 - \lfloor a_1 \rfloor]}_{f_1})x_1 + \cdots + (\lfloor a_n \rfloor + \underbrace{[a_n - \lfloor a_n \rfloor]}_{f_n})x_n = \lfloor b \rfloor + \underbrace{[b - \lfloor b \rfloor]}_{f},$$
where $ \lfloor \beta \rfloor = \max$ {$\alpha \in \mathbb{Z} : \alpha \leq \beta$} (largest integer smaller than or equal to $\beta$).
We separate fractional and integral terms:
$$\begin{align} &(\lfloor a_1 \rfloor + f_1)x_1 + \cdots + (\lfloor a_n \rfloor + f_n)x_n = \lfloor b \rfloor + f \\ \Longleftrightarrow \quad &f_1x_1 + \cdots + f_nx_n - f = \lfloor b \rfloor - \lfloor a_1 \rfloor x_1 - \cdots - \lfloor a_n \rfloor x_n \end{align}$$
Observations.
As $x_j \in \mathbb{Z}$ for all feasible x, right-hand side is integral.
Thus, for all feasible $x$, left_hand side must be integral, too.
As $ 0 \leq f < 1, x \geq 0 $ and left-hand side integral, left-hand side must be non-negative.
Consequence
$$f_1x_1 + ... + f_nx_n - f \geq 0 \quad \Longleftrightarrow \quad f_1x_1 + ... + f_nx_n \geq f \quad \text{for every feasible} \; x.$$
Suppose Step 1 of our cutting plane algorithm gives non-integral $x^*$. Then there is row in Simplex tableau with
$$ x_i^* + \sum_{j \notin I} y_{ij} x_j^* = y_{i0}$$
with $y_{i0} \notin \mathbb{Z}.$
Geory Cut.
Setting $f_j := y_{ij} - \lfloor y_{ij} \rfloor,\; f := y_{i0} - \lfloor y_{i0} \rfloor$, we get:
$$\sum_{j \notin I} f_jx_j \geq f. \qquad (*) $$
$(*)$ is fulfilled for all feasible $x$ but not for $x^*: \sum_{j \notin I} f_jx_j^* = 0 < f.$
\begin{align} &\qquad\quad x_{I_{B^{(i)}}} + \sum_{j \in I_N} \bar{a}_{ij}x_j = \bar{a}_{i0} \qquad\quad\;\; (\text{Equation 1}) \\ &\Longrightarrow \quad x_{I_{B^{(i)}}} + \sum_{j \in I_N} \lfloor \bar{a}_{ij} \rfloor x_j \leq \bar{a}_{i0} \qquad\quad (\text{since} \lfloor x \rfloor \leq x) \\ &\Longrightarrow \quad x_{I_{B^{(i)}}} + \sum_{j \in I_N} \lfloor \bar{a}_{ij} \rfloor x_j \leq \lfloor \bar{a}_{i0} \rfloor \qquad\; (\text{lhs is integer} - \text{valid for} F \text{but invalid for} F_R) \\ &\Longrightarrow \;\; \displaystyle\sum_{j \in I_N} (\bar{a}_{ij} - \lfloor \bar{a}_{ij} \rfloor)x_j \geq \bar{a}_{i0} - \lfloor \bar{a}_{i0} \rfloor \;\;\; (\text{substitute Equation 1}) \\ &\Longrightarrow \;\; \displaystyle\sum_{j \in I_N} frac(\bar{a}_{ij})x_j \geq frac(\bar{a}_{i0}) \qquad\;(frac(x) \equiv x - \lfloor x \rfloor \geq 0:\;\;\; \text{"fractional part"}) \\ &\Longrightarrow \;\; \displaystyle\sum_{j \in I_N} frac(\bar{a}_{ij})x_j - s = frac(\bar{a}_{i0}) \;\;\; (s \geq 0, s \in Z: \text{slack variable}) \end{align}
Works with ILP in standard form $\Longrightarrow$ need $A \in \mathbb{Z}^{m \times n}$ and $ b \in \mathbb{Z}^m$ in canonical form to ensure all (initial) slack variables are integer $\Longrightarrow$ multiply constraints to make $A$ and $b$ all integer before transforming to standard form.
$s = \lfloor \bar{a}_{i0} \rfloor - (x_{I_{B^{(i)}}} + \sum_{j \in I_N} \lfloor \bar{a}_{ij} \rfloor x_j = \text{(integer)} - \text{(integer) = (integer)} $
Since $ s \in \mathbb{Z}$, the problem is IP even after the addition of a Gomory's cut to the constraints $\Longrightarrow$ allow repeated "cuts" in a same manner.
Choice of $i$ is not specified $-$ typically choose $i$ with the largest fractional value of $x_i$ form the tableau.
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')