In [2]:
%qtconsole

In [4]:
import numpy as np


## The simplex algorithm¶

Implementing the algorithm seen during the lectures

In [1]:
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

In [2]:
## 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)

In [6]:
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

Out[6]:
-500.0
In [7]:
## 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)

In [8]:
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

Out[8]:
-170.0

### Alloy problem¶

In [9]:
## 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)

In [10]:
alloys(debug=False)

 optimum found!
IBV =  [3 1 0]
Bbar=  [4.00000000e-01 6.00000000e-01 1.11022302e-16]

Out[10]:
498.0

## Collisions¶

In the dual, if the solution is unbounded, this means that there is no solution in the primal

In [11]:
## 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

In [12]:
collisions(debug=False)

Solution in unbounded!


#### Unbounded solution in the dual, so no solution in the primal, so no collision !¶

In [14]:
## 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)

In [15]:
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

In [ ]:
## tighter problem formulation