%qtconsole
%load_ext autoreload
%autoreload 2
import numpy as np
Implementing the algorithm seen during the lectures
def simplexit(A,b,c, init, debug=True, expectNullOptimum=False):
"""
Implements the standard matrix form of the simplex algorithm
Hugues Talbot 2018, from an initial version in R ca. 2009.
"""
iter = 0
dim=A.shape[0]
nbvar=A.shape[1]
lastx=0
lastxl=0
sol=None
IBVmem=None # to detect cycles
vars=np.arange(0,nbvar) # indices start at 0
IBV=init # Initial basis vector
while (True):
if debug:
print("*** Iteration: ", iter)
NBV = np.setdiff1d(vars,IBV) # Set complement -> non basis variables
if debug:
print(" IBV = ", IBV)
print(" NBV = ", NBV)
B=A[:,IBV]
E=A[:,NBV]
cb = c[IBV]
ce = c[NBV]
if debug:
print(" B = \n", B)
print(" E = \n", E)
try:
iB = np.linalg.inv(B)
except:
print("B is not invertible, error")
break
if debug:
print(" iB = \n", iB)
bbar = iB.dot(b)
if (np.min(bbar) < -1e-7):
print("bbar is not positive, not a feasible solution")
break
if debug:
print(" bbar= ", bbar)
cebart = ce.T - cb.T.dot(iB).dot(E)
if debug:
print(" reduced costs", cebart)
if (np.min(cebart) >= 0):
print(" optimum found!")
sol = (IBV, bbar)
print("IBV = ", IBV)
print("Bbar= ", bbar)
if ((np.min(cebart) > 0) or expectNullOptimum):
break
if (IBVmem is None):
IBVmem = IBV
else:
if (IBV == IBVmem):
print("Cycle detected\n")
break
l = np.argmin(cebart) # may not be unique
P = iB.dot(A[:,NBV[l]])
if (debug):
print("P = ", P)
if (np.max(P) <= 0):
print("Solution in unbounded!")
break
ratios = bbar/P ## elementwise division
if (debug):
print("Ratios = ", ratios)
## Ignore negative members f P
negP = P < 0
ratios[negP] = np.inf ## this way they will not be selected as the min
nanP = np.isnan(ratios)
ratios[nanP] = np.inf ## nan are eliminated as well
if (debug):
print("Fixed ratios = ", ratios)
xl = np.min(ratios)
j = np.argmin(ratios)
lastx = xl
lastxl = NBV[l]
if (debug):
print("Variable (", IBV[j],") in the basis replaced by variable (",NBV[l],")\n")
print("xl = ",xl)
IBV[j] = NBV[l]
if (debug):
print("IBV = ", IBV)
iter += 1
## end while
return sol
## exercice 2, TD1 (production de A,B,C liƩe)
def exampleII(debug=True):
# A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
A = np.array([[-1, 2, 0, 1, 0, 0],
[ 0, -1, 1, 0, 1, 0],
[ 1, 2, 3, 0, 0, 1]])
b = np.array([ 0, 0, 35])
C = np.array([-10,-36,-44, 0, 0, 0])
IBV = np.array([3,4,5]) # index of the IBV
sol = simplexit(A,b,C,IBV,debug)
# calcul optimum
z = C[sol[0]].dot(sol[1])
return(z)
exampleII(debug=True)
*** Iteration: 0 IBV = [3 4 5] NBV = [0 1 2] B = [[1 0 0] [0 1 0] [0 0 1]] E = [[-1 2 0] [ 0 -1 1] [ 1 2 3]] iB = [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] bbar= [ 0. 0. 35.] reduced costs [-10. -36. -44.] P = [0. 1. 3.] Ratios = [ nan 0. 11.66666667] Fixed ratios = [ inf 0. 11.66666667] Variable ( 4 ) in the basis replaced by variable ( 2 ) xl = 0.0 IBV = [3 2 5] *** Iteration: 1 IBV = [3 2 5] NBV = [0 1 4] B = [[1 0 0] [0 1 0] [0 3 1]] E = [[-1 2 0] [ 0 -1 1] [ 1 2 0]] iB = [[ 1. 0. 0.] [ 0. 1. 0.] [-0. -3. 1.]] bbar= [ 0. 0. 35.] reduced costs [-10. -80. 44.] P = [ 2. -1. 5.] Ratios = [ 0. -0. 7.] Fixed ratios = [ 0. inf 7.] Variable ( 3 ) in the basis replaced by variable ( 1 ) xl = 0.0 IBV = [1 2 5] *** Iteration: 2 IBV = [1 2 5] NBV = [0 3 4] B = [[ 2 0 0] [-1 1 0] [ 2 3 1]] E = [[-1 1 0] [ 0 0 1] [ 1 0 0]] iB = [[ 0.5 0. 0. ] [ 0.5 1. 0. ] [-2.5 -3. 1. ]] bbar= [ 0. 0. 35.] reduced costs [-50. 40. 44.] P = [-0.5 -0.5 3.5] Ratios = [-0. -0. 10.] Fixed ratios = [inf inf 10.] Variable ( 5 ) in the basis replaced by variable ( 0 ) xl = 10.0 IBV = [1 2 0] *** Iteration: 3 IBV = [1 2 0] NBV = [3 4 5] B = [[ 2 0 -1] [-1 1 0] [ 2 3 1]] E = [[1 0 0] [0 1 0] [0 0 1]] iB = [[ 0.14285714 -0.42857143 0.14285714] [ 0.14285714 0.57142857 0.14285714] [-0.71428571 -0.85714286 0.28571429]] bbar= [ 5. 5. 10.] reduced costs [ 4.28571429 1.14285714 14.28571429] optimum found! IBV = [1 2 0] Bbar= [ 5. 5. 10.]
/usr/local/Cellar/ipython/7.3.0/libexec/vendor/lib/python3.7/site-packages/ipykernel_launcher.py:77: RuntimeWarning: invalid value encountered in true_divide
-500.0
## Bakery problem
def bakery2(debug=True):
# A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
A = np.array([[ 5, -1, 1, 1, 0],
[ 1, 0, 0, 0, 1]])
b = np.array([30, 5])
C = np.array([-30, 4, -4, 0, 0])
IBV = np.array([3,4]) # index of the IBV
sol = simplexit(A,b,C,IBV,debug)
# calcul optimum
z = C[sol[0]].dot(sol[1])
return(z)
bakery2(debug=False)
optimum found! IBV = [2 0] Bbar= [5. 5.] Solution in unbounded!
/Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: divide by zero encountered in true_divide
-170.0
## Alloy problem
def alloys(debug=True):
A = np.array([[ 10, 10, 40, 60, 30, 30, 30, 50, 20, 1, 0, 0],
[ 10, 30, 50, 30, 30, 40, 20, 40, 30, 0, 1, 0],
[ 80, 60, 10, 10, 40, 30, 50, 10, 50, 0, 0, 1]])
b = np.array([30, 30, 40])
C = 100*np.array([4.1, 4.3, 5.8, 6.0, 7.6, 7.5, 7.3, 6.9, 7.3, 100, 100, 100])
IBV = np.array([9,10,11]) # index of the IBV
sol = simplexit(A,b,C,IBV,debug)
# calcul optimum
z = C[sol[0]].dot(sol[1])
return(z)
alloys(debug=False)
optimum found! IBV = [3 1 0] Bbar= [4.00000000e-01 6.00000000e-01 1.11022302e-16]
498.0
In the dual, if the solution is unbounded, this means that there is no solution in the primal
## Alloy problem
def collisions(debug=True):
A = np.array([[ 1, 1, -2, -1, -1, 2, 1, 0],
[ 1, -2, 1, 1, -1, 1, 0, 1]])
b = np.array([0, 0])
C = -1*np.array([5, -10, -4, 6, -10, 10, 0, 0])
IBV = np.array([6,7]) # index of the IBV
sol = simplexit(A,b,C,IBV,debug)
# calcul optimum not needed
# z = C[sol[0]].dot(sol[1])
return
collisions(debug=False)
Solution in unbounded!
## Young couple
def youngcouple(debug=False):
A = np.array([[ 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0],
[ 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[ 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
[ 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
[ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
])
b = np.array([2, 2, 1, 1, 1, 1])
C = np.array([4.5, 7.8, 3.6, 2.9, 4.9, 7.2, 4.3, 3.1, 100, 100, 100, 100, 100, 100])
IBV = np.array([8,9,10,11,12,13])
sol = simplexit(A,b,C,IBV,debug)
z= C[sol[0]].dot(sol[1])
return(z)
print(youngcouple())
optimum found! IBV = [ 2 5 4 11 6 3] Bbar= [1. 1. 1. 0. 0. 1.] 18.6
/Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: divide by zero encountered in true_divide /Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: invalid value encountered in true_divide
## tighter problem formulation