Bierlaire, M. (2015). Optimization: Principles and Algorithms. EPFL Press.
The numbering of the algorithms, tables and page refer to the book.
import numpy as np
from scipy import linalg
# Define a function simple to call for the inner product
def scalarInner(a,b):
return(np.asscalar(a.T.dot(b)))
from optimizationExceptions import *
import itertools as itt
def vertexEnumeration(A,b,c):
m,n = A.shape
if b.shape != (m,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, b is {}x{}, and should be {}x1".format(m,n,b.shape[0],b.shape[1],m))
if c.shape != (n,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, c is {}x{}, and should be {}x1".format(m,n,c.shape[0],c.shape[1],n))
optvalue = np.inf
iters = list()
for k in itt.combinations(range(n),m):
B = A[:,k]
try:
x = linalg.solve(B,b)
xfull = np.zeros((n,1))
for i,ik in enumerate(k):
xfull[[ik],:] = x[[i],:]
iteroutput = [xfull]
if (x >= 0).all():
f = scalarInner(xfull,c)
iteroutput += [f]
if f < optvalue:
optvalue = f
solution = xfull
else:
iteroutput += [np.inf]
except linalg.LinAlgError:
iteroutput += [np.inf]
pass
iters.append(iteroutput)
return solution,optvalue,iters
Example 16.4, page 367
A = np.matrix([[1,1,1,0],[1,-1,0,1]])
b = np.array([[1],[1]])
c = np.array([[-1],[-2],[0],[0]])
sol,optvalue,iters = vertexEnumeration(A,b,c)
sol
Table p. 368
for k in range(len(iters)):
print("{}\t{}\t{}".format(k+1,iters[k][0].T,iters[k][1]))
def simplex(A,b,c,basis):
def getFullSolution(basis):
# For printing only
B = A[:,basis]
xb = linalg.solve(B,b)
xfull = np.zeros((n,1))
for i,ik in enumerate(basis):
xfull[[ik],:] = xb[[i],:]
return xfull,scalarInner(c,xfull)
m,n = A.shape
if b.shape != (m,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, b is {}x{}, and should be {}x1".format(m,n,b.shape[0],b.shape[1],m))
if c.shape != (n,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, c is {}x{}, and should be {}x1".format(m,n,c.shape[0],c.shape[1],n))
optvalue = np.inf
iters = list()
stop = False
while not stop:
B = A[:,basis]
# Each column of the matrix d is the dB vector of the corresponding index
d = linalg.solve(B,A)
reducedCost = c.T - c[basis].T.dot(d)
negativeReducedCost = reducedCost < 0
if not negativeReducedCost.any():
iters.append([basis.copy(),reducedCost.copy(),getFullSolution(basis)[0],None,np.inf,None,None,getFullSolution(basis)[1]])
optimalbasis = basis
unbounded = False
stop = True
else:
# In Python, True is larger than False. The next statement return the index of a True entry in the array, that is the index of a negative reduced cost
p = np.argmax(negativeReducedCost)
xb = linalg.solve(B,b)
xfull = np.zeros((n,1))
for i,ik in enumerate(basis):
xfull[[ik],:] = xb[[i],:]
step = np.inf
q = -1
for k in range(m):
if d[k][p] > 0 and step > 0:
astep = np.asscalar(xb[k]) / d[k][p]
if astep < step:
step = astep
q = k
iters.append([basis.copy(),reducedCost.copy(),getFullSolution(basis)[0],-d[:,p],step,p,basis[q],getFullSolution(basis)[1]])
if q == -1:
optimalbasis = np.zeros((m,1))
unbounded = True
stop = True
else:
basis[q] = p
return optimalbasis,unbounded,iters
# Python counts from 0. In the book, indices are counted from 1.
basis = [2,3]
sol,unbounded,iters = simplex(A,b,c,basis)
Example 16.6, page 373. We have renumbered the indices to be as in the book
for k in range(len(iters)):
m,n = A.shape
basis = iters[k][0]
B = A[:,basis]
print("Iteration {}".format(k+1))
print("***********")
print("Basis: {}".format([k+1 for k in basis]))
print("B={}".format(B))
x = iters[k][2]
print("x={}".format(x))
print("cbar={}".format(iters[k][1]))
print("Objective function: {}".format(scalarInner(c,x)))
p = iters[k][5]
if p is not None:
print("Index entering the basis: {}".format(1+p))
db = iters[k][3]
print("db={}".format(db))
dfull = np.zeros((n,1))
for i,ik in enumerate(basis):
dfull[[ik]] = db[[i]]
dfull[p] = 1.0
print("d={}".format(dfull))
print("alpha={}".format(iters[k][4]))
q = iters[k][6]
print("Index exiting the basis: {}".format(1+q))
Example 16.7, page 375
A = np.matrix([[1,2,2,1,0,0],[2,1,2,0,1,0],[2,2,1,0,0,1]])
b = np.array([[20],[20],[20]])
c = np.array([[-10],[-12],[-12],[0],[0],[0]])
# Python counts from 0. In the book, indices are counted from 1.
basis = [3,4,5]
sol,unbounded,iters = simplex(A,b,c,basis)
Table 16.1, page 376
for k in range(len(iters)):
m,n = A.shape
basis = iters[k][0]
cbar = iters[k][1].reshape(n,)
x = iters[k][2].reshape(n,)
obj = scalarInner(c,x)
s = "{}|{} {} {}|{:5.1f} {:5.1f} {:5.1f}|{:4.1f} {:4.1f} {:4.1f} {:4.1f} {:4.1f} {:4.1f}".format(k,*basis,*cbar[np.nonzero(cbar)],*x)
if iters[k][5] is not None:
p = iters[k][5] + 1
q = iters[k][6] + 1
db = iters[k][3]
alpha = iters[k][4]
s += "|{:4.1f} {:4.1f} {:4.1f}|{:2.0f}|{}|{}|".format(*db,alpha,p,q)
else:
s += "| | | | |"
s += "{:5.0f}".format(obj)
print(s)
import sys
def pivoting(tableau,l,j):
m,n = tableau.shape
if l >= m:
raise OptimizationInputError("The row of the pivot ({}) exceeds the size of the tableau ({})".format(l,m))
if j >= n:
raise OptimizationInputError("The column of the pivot ({}) exceeds the size of the tableau ({})".format(j,n))
thepivot = tableau[l][j]
if abs(thepivot) < sys.float_info.epsilon:
raise OptimizationValueOutOfRange("The pivot is too close to zero: {}".format(thepivot))
thepivotrow = tableau[l,:]
newtableau = np.empty(tableau.shape)
newtableau[l,:] = tableau[l,:] / thepivot
for i in range(m):
if i != l:
mult = -tableau[i][j] / thepivot
newtableau[i,:] = tableau[i,:] + mult * thepivotrow
return newtableau
def prettyTableau(tableau):
m,n = tableau.shape
s = ''
for i in range(m-1):
formattedRow = ['{:6.2f}' for k in tableau[i,:-1]]
s += '\t'.join(formattedRow).format(*tableau[i,:-1])
s += '|{:6.2f}'.format(tableau[i,-1])
s += '\n'
for i in range(n):
s += '------\t'
s += '\n'
formattedRow = ['{:6.2f}' for k in tableau[m-1,:-1]]
s += '\t'.join(formattedRow).format(*tableau[m-1,:-1])
s += '|{:6.2f}'.format(tableau[m-1,-1])
s += '\n'
return s
Example 16.12, page 380
A = np.array([[0,1.5,1,1,-0.5,0,10],[1,0.5,1,0,0.5,0,10],[0,1,-1,0,-1,1,0],[0,-7,-2,0,5,0,100]])
print(prettyTableau(A))
Apivoted = pivoting(A,2,1)
print(prettyTableau(Apivoted))
def simplexTableau(tableau,rowindex):
"""
:param tableau: the first simplex tableau
:rowindex: for each row, index of the corresponding basic variable
"""
mtab,ntab = tableau.shape
m = mtab - 1
n = ntab - 1
iters = list()
while True:
colpivot = np.argmax(tableau[-1,:] < 0)
if tableau[-1,colpivot]>= -sys.float_info.epsilon:
# The tableau is optimal
iters.append([tableau,-1,colpivot])
return tableau,rowindex,True,iters
rowpivot = -1
beststep = np.inf
for j in range(m):
if tableau[j,colpivot] > sys.float_info.epsilon:
step = tableau[j,-1] / tableau[j,colpivot]
if step < beststep:
beststep = step
rowpivot = j
if rowpivot == -1:
# The problem is unbounded
iters.append([tableau,rowpivot,colpivot])
return tableau,rowindex,False,iters
iters.append([tableau,rowpivot,colpivot])
newtableau = pivoting(tableau,rowpivot,colpivot)
tableau = newtableau.copy()
rowindex[rowpivot] = colpivot
Example 16.13, page 382
A = np.array([[1,2,2,1,0,0,20],[2,1,2,0,1,0,20],[2,2,1,0,0,1,20],[-10,-12,-12,0,0,0,0]])
optimalTableau,rowindex,bounded,iters = simplexTableau(A,[3,4,5])
print(prettyTableau(optimalTableau))
for k in range(len(iters)):
print("Iteration {}".format(k))
print(prettyTableau(iters[k][0]))
print("Pivot in row {} and column {}".format(iters[k][1],iters[k][2]))
print("\n")
def twoPhasesSimplex(A,b,c):
m,n = A.shape
if b.shape != (m,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, b is {}x{}, and should be {}x1".format(m,n,b.shape[0],b.shape[1],m))
if c.shape != (n,1):
raise OptimizationInputError("Incompatible sizes: A is {}x{}, c is {}x{}, and should be {}x1".format(m,n,c.shape[0],c.shape[1],n))
# Make sure that b >= 0
b[b < 0.0] *= -1.0
for k in range(m):
if b[k] < 0:
A[k,:] *= -1.0
# First tableau for Phase I
tableau = np.empty((m+1,n+m+1))
tableau[:m,:n] = A
tableau[:m,n:n+m] = np.eye(m)
tableau[:m,n+m:n+m+1] = b
# The last row
tableau[-1,:n] = np.array([-np.sum(tableau[:m,k]) for k in range(n)]).copy()
tableau[-1,n:n+m] = 0.0
tableau[-1,-1] = -np.sum(tableau[:m,-1])
rowIndex = np.array([n+k for k in range(m)])
# Solve the Phase I auxiliary problem
phaseOneTableau,rowIndex,bounded,itersPhaseOne = simplexTableau(tableau,rowIndex)
if not bounded:
return None,None,None,True,False
if phaseOneTableau[-1,-1] > np.sqrt(sys.float_info.epsilon):
return None,None,None,False,True
feasible = True
# Remove the auxiliary variables from the basis
clean = False
rowsToRemove = []
while not clean:
# Check if some auxiliary variables are in the basis
tobeCleaned = set(rowIndex).intersection(set(range(n,n+m)))
if tobeCleaned == set():
clean = True
else:
rowpivotIndex = list(rowIndex).index(tobeCleaned.pop())
rowpivot = phaseOneTableau[rowpivotIndex,:]
nonbasic = list(set(range(n+m))-set(rowIndex))
nonzerosPivots = abs(rowpivot[nonbasic]) > sys.float_info.epsilon
if nonzerosPivots.any():
colpivot = nonbasic[np.argmax(nonzerosPivots)]
newtableau = pivoting(phaseOneTableau,rowpivotIndex,colpivot)
phaseOneTableau = newtableau.copy()
rowIndex[rowpivotIndex] = colpivot
else:
rowsToRemove.append(rowpivotIndex)
rowIndex[rowpivotIndex] = -1
# Delete rows
tmp = np.delete(phaseOneTableau,rowsToRemove,0)
rowIndex = rowIndex[rowIndex >= 0]
m -= len(rowsToRemove)
# Delete columns
startPhaseTwo = np.delete(tmp,range(n,n+m),1)
# Calculate last row
cb = c[rowIndex]
nonbasic = set(range(n))-set(rowIndex)
for k in nonbasic:
startPhaseTwo[-1,k] = c[k] - scalarInner(cb,startPhaseTwo[:m,k])
startPhaseTwo[-1,-1] = -scalarInner(cb,startPhaseTwo[:m,-1])
# Phase II
phaseTwoTableau,rowIndex,bounded,itersPhaseTwo = simplexTableau(startPhaseTwo,rowIndex)
solution = np.zeros(n)
solution[rowIndex] = phaseTwoTableau[:m,-1]
copt = -phaseTwoTableau[-1,-1]
return x,copt,phaseTwoTableau,True,bounded,itersPhaseOne,itersPhaseTwo
Example 16.16, p. 392
A = np.matrix([[1,3,0,4,1],[1,2,0,-3,1],[-1,-4,3,0,0]],dtype=float)
b = np.array([[2],[2],[1]],dtype=float)
c = np.array([[2],[3],[3],[1],[-2]],dtype=float)
x,copt,phaseTwoTableau,feasible,bounded,itersPhaseOne,itersPhaseTwo = twoPhasesSimplex(A,b,c)
Phase I
for k in range(len(itersPhaseOne)):
print(prettyTableau(itersPhaseOne[k][0]))
Phase II
for k in range(len(itersPhaseTwo)):
print(prettyTableau(itersPhaseTwo[k][0]))