Integer Programming (IP)
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
Table of Contents
- I. 1. Integer Programming (IP)¶
- II. 2. Mixed Integer Linear Programming (MILP)¶
- III. 3. Examples of MILP¶
- IV. 4. Binary Integer Programming¶
- I. 4.1. Knapsack Problems¶
- II. 4.2. Factory and Warehouse¶
- III. 4.3. Shortest Path¶
- IV. 4.4. Vertex Cover Problem¶
- V. 4.5. Coloring Problems¶
- VI. 4.6. Flat Pattern Design of 3D Shapes¶
- VII. 4.7. Traveling Salesman Problem¶
- V. 5. Any Ideas to Solve IP (Optional)¶
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.
maxysubject to−x+y≤13x+2y≤122x+3y≤12x,y≥0x,y∈Z
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. → 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
max(30+x−75)+(90+y−95)subject to50x+24y≤40×6030x+33y≤35×60x≥75−30y≥95−90x,y∈Z
2. Mixed Integer Linear Programming (MILP)¶
A mixed-integer linear program is a problem with
- Linear objective function, cTx, 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
mincTxsubject tobl≤Ax≤bul≤x≤uxi∈Z,i∈Xi
where Xi 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
min[0−1][x1x2]subject to[−113223][x1x2]≤[11212][00]≤[x1x2]x1,x2∈Z(integer)
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.
min[0−1][x1x2]subject to[−113223][x1x2]≤[11212][00]≤[x1x2]
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
minx[−3−2−1][x1x2x3]subject to[111][x1x2x3]≤7[421][x1x2x3]=120≤[x1x2x3]x3∈{0,1}
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
max(30+x−75)+(90+y−95)subject to50x+24y≤40×6030x+33y≤35×60x≥75−30y≥95−90x,y∈Z
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
minx1+x2+x3+x4+x5+x6subject tox1+x2≥1x1+x2+x6 ≥1x3+x4≥1x3+x4+x5≥1x4+x5+x6≥1x5+x6≥1x1 x2 x3 x4 x5 x6∈{0,1}
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)=−9x1−5x2−6x3−4x4
subject to the constraints
[63520011−10100−101][x1x2x3x4]≤[9100],
where x1,x2,x3 and x4 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?
⟹ decision variables xi∈{0,1}i=1,⋯,6
maxx175x1+90x2+20x3+50x4+10x5+200x6subject to10x1+9x2+4x3+2x4+x5+20x6≤20xi∈{0,1}i=1,⋯,6
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.
max7x1+5x2+6x3+4x4subject to5x1+3x2+5x3+2x4≤10x1+x2≥1x3+x4=1x3≤x1x4≤x2x1,x2,x3,x4∈{0,1}⟹min−[7564][x1x2x3x4]subject to[5352−1−100−10100−101][x1x2x3x4]≤[10−100][0011][x1x2x3x4]=1x1,x2,x3,x4∈{0,1}
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
xi={0if not part of shortest path1if part of shortest path
Objective function
minx1+4x2+2x3+4x4+8x5+11x6+3x7
- For each node
∑xin=∑xoutxin:entering edge to node ixout:leaving edge from node i
- At node 0, it has to transit to either node 1 or node 2
x1+x2=1
- At node 4, it has to come from either node 5, 6, or 7
x5+x6+x7=1
- At node 1
x1=x3+x6
- At node 2
x2+x3=x4+x5
- At node 3
x4=x7
In the form of optimization
minx1+4x2+2x3+4x4+8x5+11x6+3x7subject to[110000010−100−10011−1−100000100−10000111][x1x2x3x4x5x6x7]=[10001]x1,x2,x3,x4,x5,x6,x7∈{0,1}
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
xi={1ATM machine allocated at the location i0no ATM machine at the location i
minx1+x2+x3+x4+x5+x6+x7+x8subject tox1+x3≥1x1+x4≥1x1+x5≥1x1+x6≥1x1+x7≥1x2+x4≥1x2+x5≥1x3+x8≥1x5+x7≥1
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},⋯
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
V={V1,V2,V3,V4,V5}E={E12,E13,E14,E23,E24,E34,E35,E45}
- Define ILP decision variables
yk={1if color k is used0if color k is not usedxik={1if vertex i gets color k0otherwise
- the coloring problem can be formulated as follows:
min5∑k=1yksubject to5∑k=1xik=1xik+xjk≤1for Eij∈E and 1≤k≤5xik≤ykfor 1≤i≤5 and 1≤k≤5xik,yk∈{0,1}
## 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
xi={1if two nodes are connected (welded).0
Constraints
No disconnected node.
Only 6-1 edges are required. (no cycle)
maxx1+x2+x3+x4+2x5+2x6+2x7+2x8+3x9+3x10+3x11+3x12subject tox1+x2+x5+x8≥1x2+x4+x10+x12≥1x3+x4+x6+x7≥1x1+x3+x9+x11≥1x5+x6+x9+x10≥1x7+x8+x11+x12≥1x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12=5
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 cij 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.
xij={1if travel from i to j0otherwise
- Optimization
minx4∑i=14∑j=1cijxijsubject to∑jxij=1for i=1,2,3,4∑ixij=1for j=1,2,3,4∑i∈S,j∈¯Sxij≥1for all prtitions(S,¯S)2≤|S|≤m−2
minx5x12+10x13+13x14+5x21+8x23+9x24+10x31+8x32+7x34+13x41+9x42+7x43subject tox11+x12+x13+x14=1x21+x22+x23+x24=1x31+x32+x33+x34=1x41+x42+x43+x44=1x11+x21+x31+x41=1x12+x22+x32+x42=1x13+x23+x33+x43=1x14+x24+x34+x44=1x11=x22=x33=x44=0(no sub-closed loop conditions)x23+x24+x13+x14≥1x31+x32+x41+x42≥1x21+x24+x31+x34≥1x12+x13+x42+x43≥1x21+x23+x41+x43≥1x12+x14+x32+x34≥1
- 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:
minimizef(x)subject tox∈F
Relaxed problem:
minimizef(x)subject tox∈FRwhereF⊆FR
Theorem : Let x∗ and x∗R be optimal solutions of the original and relaxed problems, respectively. Then, f(x∗R)≤f(x∗). If x∗R∈F,x∗R is an optimal solution of the original problem.
Proof: (First part) Suppose f(x∗R)≤f(x∗). Since x∗∈F and F⊆FR,x∗∈FR. This contradicts the fact that x∗R is an optimal solution of the relaxed problem. Therefore, f(x∗R)≤f(x∗).
(Second part) Now we know f(x∗R)≤f(x∗). If x∗R∈F, therefore, x∗R is an optimal solution of the original problem.
5.2. ILP → LP Relaxation¶
Original problem (ILP):
minimizecTxsubject toAx=bx≥0x∈Zn
Relaxed problem (LP):
minimizecTxsubject toAx=bx≥0x∈Rn
Zn⊆Rn⟹ the above theorem applies.
For minimization problem, cTx∗R gives a lower bound to the optimal function value of the original ILP problem.
If x∗R∈Zn, we have solved the ILP problem.
In general, round(x∗R)≠x∗!
5.3. Cutting Plane Algorithms¶
General case
dTx≤e is a valid inequality for X⊆Rn if ∀x∈X,dTx≤e.
If dTx≤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∉F, add an inequality valid for F but invalid for FR to the constraints and solve the relaxed problem again. Repeat until x∗R∈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: xj∈Z→xj∈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 ⟹ 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:
max5x1+8x2subject tox1+x2≤65x1+9x2≤45x1,x2≥0x1,x2∈Z.
Step 1. Solve the IP problem with continous variables instead of discrete ones.
Step 2. If the resulting optimal solution x∗ is integral, stop ⟹ 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 2x1+3x2≤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 ⇒ optimal solution found.
x∗ = (3,3) is integral ⇒ 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 x1,...,xn≥0 and integral. We show how to construct a Gomory Cut for
a1x1+⋯+anxn=b,
where aj,b∈R (not necessarily integral). Note that this can be written as
(⌊a1⌋+[a1−⌊a1⌋]⏟f1)x1+⋯+(⌊an⌋+[an−⌊an⌋]⏟fn)xn=⌊b⌋+[b−⌊b⌋]⏟f,
where ⌊β⌋=max {α∈Z:α≤β} (largest integer smaller than or equal to β).
We separate fractional and integral terms:
(⌊a1⌋+f1)x1+⋯+(⌊an⌋+fn)xn=⌊b⌋+f⟺f1x1+⋯+fnxn−f=⌊b⌋−⌊a1⌋x1−⋯−⌊an⌋xn
Observations.
As xj∈Z for all feasible x, right-hand side is integral.
Thus, for all feasible x, left_hand side must be integral, too.
As 0≤f<1,x≥0 and left-hand side integral, left-hand side must be non-negative.
Consequence
f1x1+...+fnxn−f≥0⟺f1x1+...+fnxn≥ffor every feasiblex.
Suppose Step 1 of our cutting plane algorithm gives non-integral x∗. Then there is row in Simplex tableau with
x∗i+∑j∉Iyijx∗j=yi0
with yi0∉Z.
Geory Cut.
Setting fj:=yij−⌊yij⌋,f:=yi0−⌊yi0⌋, we get:
∑j∉Ifjxj≥f.(∗)
(∗) is fulfilled for all feasible x but not for x∗:∑j∉Ifjx∗j=0<f.
xIB(i)+∑j∈INˉaijxj=ˉai0(Equation 1)⟹xIB(i)+∑j∈IN⌊ˉaij⌋xj≤ˉai0(since⌊x⌋≤x)⟹xIB(i)+∑j∈IN⌊ˉaij⌋xj≤⌊ˉai0⌋(lhs is integer−valid forFbut invalid forFR)⟹∑j∈IN(ˉaij−⌊ˉaij⌋)xj≥ˉai0−⌊ˉai0⌋(substitute Equation 1)⟹∑j∈INfrac(ˉaij)xj≥frac(ˉai0)(frac(x)≡x−⌊x⌋≥0:"fractional part")⟹∑j∈INfrac(ˉaij)xj−s=frac(ˉai0)(s≥0,s∈Z:slack variable)
Works with ILP in standard form ⟹ need A∈Zm×n and b∈Zm in canonical form to ensure all (initial) slack variables are integer ⟹ multiply constraints to make A and b all integer before transforming to standard form.
s=⌊ˉai0⌋−(xIB(i)+∑j∈IN⌊ˉaij⌋xj=(integer)−(integer) = (integer)
Since s∈Z, the problem is IP even after the addition of a Gomory's cut to the constraints ⟹ allow repeated "cuts" in a same manner.
Choice of i is not specified − typically choose i with the largest fractional value of xi form the tableau.
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')