Which we can write in standard form as
\begin{align*} \max x_1 + 2x_2\\ -x_1 + 2x_2 \leq 2\\ x_1 + x_2 \leq 3\\ x_1, x_2 \in \mathbb{Z}^+ \end{align*}%run preamble.py
def find_pivots(A):
col_pivot = max((x,i) for i,x in enumerate(A[-1,:-2]))
ab_col = zip(A[:-1,col_pivot[1]],A[:-1,-1])
row_pivot_cands = list((bi/ai, i) for i, (ai, bi) in enumerate(ab_col) if ai > 0)
row_pivot = min(row_pivot_cands)
print("row pivot cands:", row_pivot_cands)
print("col pivot:", col_pivot[1], "row pivot:", row_pivot[1])
return col_pivot[1], min(row_pivot_cands)[1]
A = array([
[-1, 2, 1, 0, 0, 2],
[ 1, 1, 0, 1, 0, 3],
[ 1, 2, 0, 0, 1, 0]
], dtype=float)
tableau(A)
c, r = find_pivots(A)
A[r,:] *= 1/2
A[1,:] -= A[r,:]
A[2,:] -= 2*A[r,:]
tableau(A)
c, r = find_pivots(A)
A[r,:] *= 1/(3/2)
A[0,:] += 1/2*A[r,:]
A[2,:] -= 2*A[r,:]
tableau(A)
We see that it is not an integer solution. Lets choose the first row to cut.
As from the lecture notes we can write the cut as: $$ \sum_{j\in N} f_{uj} x_j \geq f_u $$
#index of non basic variables
N = [2,3,5]
row = A[0,N]
print(row)
from math import floor
#lets write the constraint:
row = [str(f(a - floor(a)).limit_denominator()) for a in row]
print(row)
So we have the constraint $$ 1/3 x_3 + 1/3 x_4 \geq 2/3 $$
We can write the constraint in form of the original variables: \begin{align*} -x_1 + 2x_2 + x_3 = 2 \Rightarrow x_3 = 2 + x_1 - 2x_2\\ x_1 + x_2 + x_4 = 3 \Rightarrow x_4 = 3 - x_1 - x_2 \end{align*}
$$ 1/3 (2 + x_1 - 2x_2) + 1/3 (3 - x_1 - x_2) \geq 2/3 $$which gives: $$ -x_2 \geq -1 \Rightarrow x_2 \leq 1 $$
A1 = A.copy()
A1 = insert(A1, 2,[0, 0, -1/3, -1/3, 0, -2/3], axis = 0)
A1 = insert(A1, 4, [0,0,1,0], axis = 1)
tableau(A1)
def find_dual_pivots(A):
row_pivot = min((x,i) for i,x in enumerate(A[:-1,-1]))
print(row_pivot)
ac_col = zip(A[row_pivot[1],:-2],A[-1,:-2])
col_pivot_cands = list((abs(ci/ai), i) for i, (ai, ci) in enumerate(ac_col) if ai < 0)
col_pivot = min(col_pivot_cands)
print("col pivot cands:", col_pivot_cands)
print("col pivot:", col_pivot[1], "row pivot:", row_pivot[1])
return col_pivot[1], min(col_pivot_cands)[1]
A2 = A1.copy()
c, r = find_dual_pivots(A2)
A2[r,:] *= 1/(-1/3)
A2[0,:] += -1/3 * A2[r,:]
A2[1,:] += 1/3 * A2[r,:]
A2[3,:] += 1/3 * A2[r,:]
tableau(A2)
Hence we get the optimal solution:
$$ p = 3.3; x1 = 1, x2 = 0.5, x3 = 0, x4 = 0, x5 = 1 $$Since x2 is fractional, we can bound it by either setting $x2 =0$ or $x2=1$. By setting $x2=1$ we break the first constraint. However, the solution [1,0,0,0,1] is feasible with the value 2.4. We cannow say it is optimal since 3.3-2.4 is not 0.
Now we can branch. We simply add these constraints to the solution and get the optimal solutions:
SP1, x2=0 $$ p = 3.28; x1 = 1, x2 = 0, x3 = 0.2, x4 = 1, x5 = 1 $$
SP, 2x2=1 $$ p = 3.2; x1 = 0, x2 = 1, x3 = 0, x4 = 0, x5 = 1 $$
We cannot say that either of these are optimal, since the optimality gap is not closed. However, x2 is an integer solution, and gives us a lower bound!
We see that SP1 is active since it has the highest objective value. we check if we can choose x3 or not:
x3=0: $$ p = 3; x1 = 1, x2 = 0, x3 = 0, x4 = 1, x5 = 1 $$
x3=1 $$ p = 3.13333; x1 = 0.333333, x2 = 0, x3 = 1, x4 = 0, x5 = 1 $$
we see that both values are less than our lower bound and can hence be pruned.
from graphviz import *
ps = Digraph(name='tree', node_attr={})
ps.node("SP0", label="3.3")
ps.node("SP1", label="3.28")
ps.edge("SP0", "SP1", label="x2<=0")
ps.node("SP2", label="3.2*")
ps.edge("SP0", "SP2", label="x2>=1")
ps.node("SP3", label="3*")
ps.node("SP4", label="3.13")
ps.edge("SP1", "SP3", label="x3<=0")
ps.edge("SP1", "SP4", label="x3>=1")
ps