Order of exercises: 11, 12, 1, 5, 9
from numpy import *
from fractions import Fraction as f
set_printoptions(precision=3,suppress=True)
def printm(a):
"""Prints the array as strings
:a: numpy array
:returns: prints the array
"""
def p(x):
return str(x)
p = vectorize(p,otypes=[str])
print(p(a))
def tableau(a,W=7):
"""Returns a string for verbatim printing
:a: numpy array
:returns: a string
"""
if len(a.shape) != 2:
raise ValueError('verbatim displays two dimensions')
rv = []
rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
rv+=[r'|'+'|'.join(map(lambda i: '{0:>{width}}'.format("x"+str(i+1)+" ",width=W), range(a.shape[1]-2)) )+"|"+
'{0:>{width}}'.format("-z ",width=W)+"|"
'{0:>{width}}'.format("b ",width=W)+"|"]
rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
for i in range(a.shape[0]-1):
rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
i = a.shape[0]-1
rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
print('\n'.join(rv))
The ILP:
$\begin{aligned} \text{max} & \qquad x_1+2x_2 & && \\ \text{s.t}& \qquad x_1-2x_2\ge -2 \\ & \qquad x_1+x_2 \le 3 \\ & \qquad x_1,x_2\ge 0 \ \ \text{and integer} \end{aligned}$
Retwritten:
$\begin{aligned} \text{max} & \qquad x_1+2x_2 & && \\ \text{s.t}& \qquad -x_1+2x_2\le 2 \\ & \qquad x_1+x_2 \le 3 \\ & \qquad x_1,x_2\ge 0 \ \ \text{and integer} \end{aligned}$
Use Gomory's fractional cutting plane algorithm:
A = array([[-1,f(2,1),1,0,0,2],[1,1,0,1,0,3],[1,2,0,0,1,0]])
tableau(A)
Find pivot, let it be column 2. Check for row 1 and row 2:
print("row 1:",2/2)
print("row 2:",3/1)
Then the pivot is row 1 column 2:
A[0,:]=f(1,2)*A[0,:]
A[1,:]=A[1,:]-A[0,:]
A[2,:]=A[2,:]-2*A[0,:]
tableau(A)
The new pivot is column 1 row 2:
A[1,:]=f(2,3)*A[1,:]
A[0,:]=A[0,:]+f(1,2)*A[1,:]
A[2,:]=A[2,:]-2*A[1,:]
tableau(A)
Let the cut be based on row 1, then the cut is:
$ax_3+bx_4\ge c$
With:
$a = 1/3 - \lfloor 1/3 \rfloor =1/3$
$b = 1/3 - \lfloor 1/3 \rfloor =1/3$
$c = 5/3 - \lfloor5/3\rfloor =2/3$
Hence we have the cut:
$1/3x_3+1/3x_4\ge 2/3 \implies -1/3x_3-1/3x_4 \le -2/3$
import sympy as sy
x1 = sy.Symbol('x_1')
x2 = sy.Symbol('x_2')
sy.init_printing()
1/3*(-x1+2*x2-2)+1/3*(x1+x2-3)<=0
Since $x_2$ is integer, then this would mean: $x_2\le 1$.
And we add the cut with another slack variable $x_5$:
$ -1/3x_3-1/3x_4+x_5 = -2/3$
B = array([[0,0,-f(1,3),-f(1,3),1,0,-f(2,3)],
[0,1,f(1,3),f(1,3),0,0,f(5,3)],
[1,0,-f(1,3),f(2,3),0,0,f(4,3)],
[0,0,-f(1,3),-f(4,3),0,1,-f(14,3)]])
tableau(B)
Use dual simplex with rule $\min|c_j/a_{ij}|$:
print("col 3:",abs((-1/3)/(-1/3)))
print("col 4:", abs((-4/3)/(-1/3)))
We take min of the two, so pivot is row 1 col 3.
B[0,:]=-3*B[0,:]
B[1,:]=B[1,:]-f(1,3)*B[0,:]
B[2,:]=B[2,:]+f(1,3)*B[0,:]
B[3,:]=B[3,:]+f(1,3)*B[0,:]
tableau(B)
Hence the tableau is optimal with integer values: $x_1=2, \ x_2=1$ with $z^*=4$.
Let $I = \{1,2,3,4,5\}$ be the set of projects.
Let $x_i$ be the amount invested in project $i\in I$.
Let $c_i$ be the the cost for project $i\in I$.
Let $a_i$ be the evaluation score for project $i\in I$.
The ILP:
$\begin{aligned} \text{max} & \qquad \sum\limits_{i\in I}a_ix_i & && \\ \text{s.t}& \qquad \sum\limits_{i\in I}c_ix_i \le 20 \\ & \qquad x_2+x_3 \le 1 \\ & \qquad x_i \in \{0,1\} & \forall i \in I \end{aligned}$
The LP relaxation:
The LP:
$\begin{aligned} \text{max} & \qquad \sum\limits_{i\in I}a_ix_i & && \\ \text{s.t}& \qquad \sum\limits_{i\in I}c_ix_i \le 20 \\ & \qquad x_2+x_3 \le 1 \\ & \qquad 0\le x_i \le 1 & \forall i \in I \end{aligned}$
Solved using http://www.zweigmedia.com/simplex/simplex.php:
maximize z = x1 + 1.8x2 + 1.4x3 + 0.6x4 + 1.4x5 subject to
6x1 + 12x2 + 10x3 + 4x4 + 8x5 <= 20
x2 + x3 <= 1
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
x5 <= 1
Optimal solution: z = 3.3; x1 = 1, x2 = 0.5, x3 = 0, x4 = 0, x5 = 1
Solved in Scip:
from pyscipopt import *
import numpy as np
m = Model("activity")
a = [1,1.8,1.4,0.6,1.4]
c = [6,12,10,4,8]
ir = range(0,5)
x = {}
# Variables
for i in ir:
x[i] = m.addVar(vtype = "C", name = "x(%s)" %i)
# Objective
m.setObjective(quicksum(a[i]*x[i] for i in ir),"maximize")
# Constraints
m.addCons(quicksum(c[i]*x[i] for i in ir) <= 20)
m.addCons(x[1]+x[2] <= 1)
for i in ir:
m.addCons(x[i]<=1)
m.addCons(x[i]>=0)
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
print("x[%s] = %s" % (i+1,m.getVal(x[i])))
We get the solution $x=[1,0.5,0,0,1]$ with objective value $z^*=3.3$.
The rounding heuristic: If we round we get either $x_2 = 0$ or $x_2 = 1 $. Letting $x_2 = 1$ we get the solution $x'=[1,1,0,0,1]$ but it is infeasible since the cost is $6+12+8=26$ and the budget is 20.
Letting $x_2 = 0$ we get the solution $x'=[1,0,0,0,1]$ it is feasible with cost $ 6+8=14$ with $z=3$. We don't know if it is optimal or not (so 1. and 2. is not true), which means 3. is true, hence $x'$ might be optimal.
The branch and bound algorithm looks at $x_2$ since this is the only variable that has a fractional value. So it is project 2.
Let $S$ start the search tree being the solution to the LP-problem.
Let $SP1$ be the search tree that does not choose $x_2$, hence we add $x_2\le 0$ and get the following:
maximize z = x1 + 1.8x2 + 1.4x3 + 0.6x4 + 1.4x5 subject to
6x1 + 12x2 + 10x3 + 4x4 + 8x5 <= 20
x2 + x3 <= 1
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
x5 <= 1
x2 <= 0
Optimal solution: z = 3.28; x1 = 1, x2 = 0, x3 = 0.2, x4 = 1, x5 = 1
Using SCIP:
from pyscipopt import *
import numpy as np
m = Model("SP1")
a = [1,1.8,1.4,0.6,1.4]
c = [6,12,10,4,8]
ir = range(0,5)
x = {}
# Variables
for i in ir:
x[i] = m.addVar(vtype = "C", name = "x(%s)" %i)
# Objective
m.setObjective(quicksum(a[i]*x[i] for i in ir),"maximize")
# Constraints
m.addCons(quicksum(c[i]*x[i] for i in ir) <= 20)
m.addCons(x[1]<=0)
m.addCons(x[1]+x[2] <= 1)
for i in ir:
m.addCons(x[i]<=1)
m.addCons(x[i]>=0)
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
print("x[%s] = %s" % (i+1,m.getVal(x[i])))
So node $SP1$ have $x=[1,0,0.2,1,1], \ z^*=3.28$.
Next we have $SP2$ where we choose $x_2$ hence we add instead $x_2\ge 1$ and get the following:
maximize z = x1 + 1.8x2 + 1.4x3 + 0.6x4 + 1.4x5 subject to
6x1 + 12x2 + 10x3 + 4x4 + 8x5 <= 20
x2 + x3 <= 1
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
x5 <= 1
x2 >= 1
Optimal solution: z = 3.2; x1 = 0, x2 = 1, x3 = 0, x4 = 0, x5 = 1
Using SCIP:
from pyscipopt import *
import numpy as np
m = Model("SP2")
a = [1,1.8,1.4,0.6,1.4]
c = [6,12,10,4,8]
ir = range(0,5)
x = {}
# Variables
for i in ir:
x[i] = m.addVar(vtype = "C", name = "x(%s)" %i)
# Objective
m.setObjective(quicksum(a[i]*x[i] for i in ir),"maximize")
# Constraints
m.addCons(quicksum(c[i]*x[i] for i in ir) <= 20)
m.addCons(x[1]>=1)
m.addCons(x[1]+x[2] <= 1)
for i in ir:
m.addCons(x[i]<=1)
m.addCons(x[i]>=0)
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
print("x[%s] = %s" % (i+1,m.getVal(x[i])))
So node $SP2$ have $x=[0,1,0,0,1], \ z^* =3.2$.
The node $SP2$ is not active since it is an integer solution but $SP1$ is still active since it contains a fractional value for $x_3$. This solution also provides us with a lower bound, since it is integer.
The branch and bound ends when there are no more active nodes. We continue with $SP1$ and start with $x_3 \le 0$, call this $SP3$.
maximize z = x1 + 1.8x2 + 1.4x3 + 0.6x4 + 1.4x5 subject to
6x1 + 12x2 + 10x3 + 4x4 + 8x5 <= 20
x2 + x3 <= 1
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
x5 <= 1
x2 <= 0
x3 <= 0
Optimal solution: z = 3; x1 = 1, x2 = 0, x3 = 0, x4 = 1, x5 = 1
from pyscipopt import *
import numpy as np
m = Model("SP3")
a = [1,1.8,1.4,0.6,1.4]
c = [6,12,10,4,8]
ir = range(0,5)
x = {}
# Variables
for i in ir:
x[i] = m.addVar(vtype = "C", name = "x(%s)" %i)
# Objective
m.setObjective(quicksum(a[i]*x[i] for i in ir),"maximize")
# Constraints
m.addCons(quicksum(c[i]*x[i] for i in ir) <= 20)
m.addCons(x[1]<=0)
m.addCons(x[2]<=0)
m.addCons(x[1]+x[2] <= 1)
for i in ir:
m.addCons(x[i]<=1)
m.addCons(x[i]>=0)
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
print("x[%s] = %s" % (i+1,m.getVal(x[i])))
So node $SP3$ have $x=[1,0,0,1,1], \ z^*=3$.
And let $SP4$ be $x_3 \ge 1$:
maximize z = x1 + 1.8x2 + 1.4x3 + 0.6x4 + 1.4x5 subject to
6x1 + 12x2 + 10x3 + 4x4 + 8x5 <= 20
x2 + x3 <= 1
x1 <= 1
x2 <= 1
x3 <= 1
x4 <= 1
x5 <= 1
x2 <= 0
x3 >= 1
Optimal solution: z = 3.13333; x1 = 0.333333, x2 = 0, x3 = 1, x4 = 0, x5 = 1
from pyscipopt import *
import numpy as np
m = Model("SP4")
a = [1,1.8,1.4,0.6,1.4]
c = [6,12,10,4,8]
ir = range(0,5)
x = {}
# Variables
for i in ir:
x[i] = m.addVar(vtype = "C", name = "x(%s)" %i)
# Objective
m.setObjective(quicksum(a[i]*x[i] for i in ir),"maximize")
# Constraints
m.addCons(quicksum(c[i]*x[i] for i in ir) <= 20)
m.addCons(x[1]<=0)
m.addCons(x[2]>=1)
m.addCons(x[1]+x[2] <= 1)
for i in ir:
m.addCons(x[i]<=1)
m.addCons(x[i]>=0)
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
print("x[%s] = %s" % (i+1,m.getVal(x[i])))
So node $SP4$ has $x=[0.333,0,1,0,1], \ z^*=3.133$.
Since the lower bound has value $3.2$ we can prune $SP4$ by bounding.
And the optimal solution is node $SP2$ with $x=[0,1,0,0,1], \ z^* =3.2$
Let $I = \{A,B,C\}$ be the set of source towns.
Let $J = \{D,E,F,G\}$ be the set of sink towns
Let $f_i$ for $i\in I \cup J $ be the current amount of cars in town $i$.
Let $d_i$ for $i\in I \cup J$ be the desired amount of cars in town $i$.
Let $x_{ij}$ be an integer variable that is the number of cars moved from $i\in I$ to $j\in I$.
Let $c_{ij}$ be the cost to move one car from the towns $i\in I$ to the towns $j\in I$.
The IP:
$\begin{aligned} \text{min} & \qquad z=\sum\limits_{i\in I}\sum\limits_{j\in I}c_{ij}x_{ij} & && (1) \\ \text{s.t}& \qquad \sum\limits_{j\in J}x_{ij} = f_i-d_i & \forall i\in I & &(2) \\ & \qquad -\sum\limits_{i\in I}x_{ij} = f_j-d_j & \forall j \in J && (3) \\ & \qquad x_{ij} \ge 0 \ \text{and integer} & \forall i\in I, \ \forall j\in J & &(4) \end{aligned}$
The objective $(1)$ ensures that we minimizes the cost for transporting the cars.
The constraints $(2)$ ensures that sources $i\in I$ sends enough cars to reach their desired amount.
The constraints $(3)$ ensures that the sinks $j\in J$ receives enough cars to reach their desired amount.
The constraints $(4)$ ensures that we dont transport negative amount of cars and that we dont transport fractions of a car.
This can be written as a min cost flow problem.
Balance:
d = [30,40,55,60,80,40,55]
f = [65,90,95,15,60,10,25]
for i in range(len(d)):
print(chr(65+i),":",f[i]-d[i])
from pyscipopt import *
import numpy as np
m = Model("activity")
d = [30,40,55,60,80,40,55]
f = [65,90,95,15,60,10,25]
c = [[5,6,10,9],
[9,11,9,15],
[12,10,14,15]]
ir = range(0,3)
jr = range(3,7)
x = {}
# Variables
for i in ir:
for j in jr:
x[i,j] = m.addVar(vtype = "I", name = "x(%s,%s)" %(i,j))
# Objective
m.setObjective(quicksum(c[i][j-3]*x[i,j] for i in ir for j in jr),"minimize")
# Constraints
for i in ir:
m.addCons(quicksum(x[i,j] for j in jr) ==f[i]-d[i])
for j in jr:
m.addCons(-quicksum(x[i,j] for i in ir) == f[j]-d[j])
m.optimize()
print("obj =", m.getObjVal())
print()
for i in ir:
for j in jr:
if(m.getVal(x[i,j])>0):
print("x[%s,%s] = %s" % (i,j,m.getVal(x[i,j])))
Let $z=X_1\vee X_2$, then:
$\begin{aligned} & \qquad z \ge X_1 \\ & \qquad z \ge X_2 \\ & \qquad z\le X_1+X_2 \\ & \qquad X_1, X_2,z \in \{0,1\} \end{aligned}$
Let $z=X_1 \implies X_2$, then:
$\begin{aligned} & \qquad z\ge 1-X_1-X_2 \\ & \qquad z\ge 2X_2-1 \\ & \qquad z\le 1+2X_2-X_1 \\ & \qquad X_1,X_2,z\in \{0,1\} \end{aligned}$
We can write the table consisting of the values of $X_1, X_2$ and what $z$ should be. Then you can also add the constraints and check how they bound the $z-$value with the different $X_1,X_2-$values.
$X_1$ | $X_2$ | $z$ | $z\ge$ | $z\ge$ | $z\le$ |
---|---|---|---|---|---|
0 | 0 | 1 | 1 | -1 | 1 |
0 | 1 | 1 | 0 | 1 | 3 |
1 | 0 | 0 | 0 | -1 | 0 |
1 | 1 | 1 | -1 | 1 | 2 |
Let $z = X_1 \wedge X_2 $ then:
$\begin{aligned} & \qquad z\le X_1 \\ & \qquad z \le X_2 \\ & \qquad z \ge X_1+X_2-1 \\ & \qquad X_1, X_2,z \in \{0,1\} \end{aligned}$
Let $z = X_1 \iff X_2 $ then:
$\begin{aligned} & \qquad z\ge X_1+X_2-1 \\ & \qquad z \ge 1-X_1-X_2 \\ & \qquad z \le 2-2X_1+X_2 \\ & \qquad z \le 2-2X_2+X_1 \\ & \qquad X_1, X_2,z \in \{0,1\} \end{aligned}$
Let $z = \neg X_1 $ then:
$\begin{aligned} & \qquad z\le 1-X_1 \\ & \qquad z \ge 1-X_1 \\ & \qquad X_1\in \{0,1\} \end{aligned}$
Let $z = X_1 \oplus X_2$ then:
$\begin{aligned} & \qquad z\le X_1+X_2 \\ & \qquad z\le 2-X_1-X_2 \\ & \qquad z \ge X_1-X_2 \\ & \qquad z \ge X_2-X_1 \\ & \qquad X_1,X_2,z \in \{0,1\} \end{aligned}$
We start with $x_1^3$, this can be linearized by introducing $Y_1 $ with:
$\begin{aligned} & \qquad Y_1 \le x_1 \\ & \qquad Y_1 \ge x_1 \\ \end{aligned}$
We can do the same for $x_2^5, \ x_3^2$ with $Y_2, Y_3$.
$\begin{aligned} & \qquad Y_2 \le x_2 \\ & \qquad Y_2 \ge x_2 \\ \\ & \qquad Y_3 \le x_3 \\ & \qquad Y_3 \ge x_3 \\ \end{aligned}$
For $x_1x_2x_3^4$ we can use $Y_{123}$:
$\begin{aligned} & \qquad Y_{123}\le x_1 \\ & \qquad Y_{123}\le x_2 \\ & \qquad Y_{123}\le x_3 \\ & \qquad Y_{123} \ge x_1+x_2+x_3-2 \\ \end{aligned}$
For $x_1x_2, \ x_1x_2^2, \ x_1^2x_2, \ x_1x_2$ we can use $Y_{12}$
$\begin{aligned} & \qquad Y_{12}\le x_1 \\ & \qquad Y_{12}\le x_2 \\ & \qquad Y_{12} \ge x_1+x_2-1 \\ \end{aligned}$
So we can write the 0-1 ILP:
The ILP:
$\begin{aligned} \text{max} & \qquad Y_1+Y_2+Y_3+5Y_{123}-2Y_{12} & && \\ \text{s.t}& \qquad x_1+x_2+x_3 \ge 1 \\ & \qquad 7Y_{12}+3Y_{12}-5Y_{12} \ge 0 \\ & \qquad -2Y_1+3Y_2+Y_3 \le 3 \\ & \qquad Y_1,Y_2,Y_3, Y_{12}, Y_{123} \in \{0,1\} \end{aligned}$
Rewritten:
The ILP:
$\begin{aligned} \text{max} & \qquad Y_1+Y_2+Y_3+5Y_{123}-2Y_{12} & && \\ \text{s.t}& \qquad x_1+x_2+x_3 \ge 1 \\ & \qquad 5Y_{12} \ge 0 \\ & \qquad -2Y_1+3Y_2+Y_3 \le 3 \\ & \qquad Y_1,Y_2,Y_3, Y_{12}, Y_{123} \in \{0,1\} \end{aligned}$
We have: $(x_1 \vee x_2 \vee x_3 ) \wedge ( x_1 \vee \neg x_2 ) \wedge ( x_2 \vee \neg x_3 ) \wedge ( x_3 \vee \neg x_1 ) \wedge ( \neg x_1 \vee \neg x_2 \vee \neg x_3 )$
We have to satisfy all clauses.
For $I=\{1,2,3\}$ we can either assign $x_i=1$ or $\neg x_i = 1$ for all $i\in I$.
The ILP:
$\begin{aligned} \text{max} & \qquad 1 & && \\ \text{s.t}& \qquad x_1+x_2+x_3 \ge 1 \\ & \qquad x_1+(1-x_2) \ge 1 \\ & \qquad x_2+(1-x_3) \ge 1 \\ & \qquad x_3+(1-x_1) \ge 1 \\ & \qquad (1-x_1)+(1-x_2)+(1-x_3) \ge 1 \\ & \qquad x_1,x_2,x_3 \in \{0,1\} \end{aligned}$