# Optimization: principles and algorithms¶

Bierlaire, M. (2015). Optimization: Principles and Algorithms. EPFL Press.
The numbering of the algorithms, tables and page refer to the book.

## Algorithm 16.1: Vertex enumeration¶

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

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

In [3]:
A = np.matrix([[1,1,1,0],[1,-1,0,1]])
b = np.array([[1],[1]])
c = np.array([[-1],[-2],[0],[0]])

In [4]:
sol,optvalue,iters = vertexEnumeration(A,b,c)
sol

Out[4]:
array([[0.],
[1.],
[0.],
[2.]])

Table p. 368

In [5]:
for k in range(len(iters)):
print("{}\t{}\t{}".format(k+1,iters[k][0].T,iters[k][1]))

1	[[ 1. -0.  0.  0.]]	-1.0
2	[[ 1.  0. -0.  0.]]	-1.0
3	[[1. 0. 0. 0.]]	-1.0
4	[[ 0. -1.  2.  0.]]	inf
5	[[0. 1. 0. 2.]]	-2.0
6	[[0. 0. 1. 1.]]	0.0


### Algorithm 16.2: Simplex method¶

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

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

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

Iteration 1
***********
Basis: [3, 4]
B=[[1 0]
[0 1]]
x=[[0.]
[0.]
[1.]
[1.]]
cbar=[[-1. -2.  0.  0.]]
Objective function: 0.0
Index entering the basis: 1
db=[-1. -1.]
d=[[ 1.]
[ 0.]
[-1.]
[-1.]]
alpha=1.0
Index exiting the basis: 3
Iteration 2
***********
Basis: [1, 4]
B=[[1 0]
[1 1]]
x=[[1.]
[0.]
[0.]
[0.]]
cbar=[[ 0. -1.  1.  0.]]
Objective function: -1.0
Index entering the basis: 2
db=[-1.  2.]
d=[[-1.]
[ 1.]
[ 0.]
[ 2.]]
alpha=1.0
Index exiting the basis: 1
Iteration 3
***********
Basis: [2, 4]
B=[[ 1  0]
[-1  1]]
x=[[0.]
[1.]
[0.]
[2.]]
cbar=[[1. 0. 2. 0.]]
Objective function: -2.0


Example 16.7, page 375

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

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

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


0|3 4 5|-10.0 -12.0 -12.0| 0.0  0.0  0.0 20.0 20.0 20.0|-1.0 -2.0 -2.0|10|1|5|    0
1|3 0 5| -7.0  -2.0   5.0|10.0  0.0  0.0 10.0  0.0  0.0|-1.5 -0.5 -1.0| 0|2|6| -100
2|3 0 1| -9.0  -2.0   7.0|10.0  0.0  0.0 10.0  0.0  0.0|-2.5 -1.5  1.0| 4|3|4| -100
3|2 0 1|  3.6   1.6   1.6| 4.0  4.0  4.0  0.0  0.0  0.0|              |  | | | -136


### Algorithm 16.3: Tableau pivoting¶

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


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

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

  0.00	  1.50	  1.00	  1.00	 -0.50	  0.00| 10.00
1.00	  0.50	  1.00	  0.00	  0.50	  0.00| 10.00
0.00	  1.00	 -1.00	  0.00	 -1.00	  1.00|  0.00
------	------	------	------	------	------	------
0.00	 -7.00	 -2.00	  0.00	  5.00	  0.00|100.00


In [15]:
Apivoted = pivoting(A,2,1)

In [16]:
print(prettyTableau(Apivoted))

  0.00	  0.00	  2.50	  1.00	  1.00	 -1.50| 10.00
1.00	  0.00	  1.50	  0.00	  1.00	 -0.50| 10.00
0.00	  1.00	 -1.00	  0.00	 -1.00	  1.00|  0.00
------	------	------	------	------	------	------
0.00	  0.00	 -9.00	  0.00	 -2.00	  7.00|100.00



### Algorithm 16.4: Simplex algorithm¶

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

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

  0.00	  0.00	  1.00	  0.40	  0.40	 -0.60|  4.00
1.00	  0.00	  0.00	 -0.60	  0.40	  0.40|  4.00
0.00	  1.00	  0.00	  0.40	 -0.60	  0.40|  4.00
------	------	------	------	------	------	------
0.00	  0.00	  0.00	  3.60	  1.60	  1.60|136.00


In [19]:
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")

Iteration 0
1.00	  2.00	  2.00	  1.00	  0.00	  0.00| 20.00
2.00	  1.00	  2.00	  0.00	  1.00	  0.00| 20.00
2.00	  2.00	  1.00	  0.00	  0.00	  1.00| 20.00
------	------	------	------	------	------	------
-10.00	-12.00	-12.00	  0.00	  0.00	  0.00|  0.00

Pivot in row 1 and column 0

Iteration 1
0.00	  1.50	  1.00	  1.00	 -0.50	  0.00| 10.00
1.00	  0.50	  1.00	  0.00	  0.50	  0.00| 10.00
0.00	  1.00	 -1.00	  0.00	 -1.00	  1.00|  0.00
------	------	------	------	------	------	------
0.00	 -7.00	 -2.00	  0.00	  5.00	  0.00|100.00

Pivot in row 2 and column 1

Iteration 2
0.00	  0.00	  2.50	  1.00	  1.00	 -1.50| 10.00
1.00	  0.00	  1.50	  0.00	  1.00	 -0.50| 10.00
0.00	  1.00	 -1.00	  0.00	 -1.00	  1.00|  0.00
------	------	------	------	------	------	------
0.00	  0.00	 -9.00	  0.00	 -2.00	  7.00|100.00

Pivot in row 0 and column 2

Iteration 3
0.00	  0.00	  1.00	  0.40	  0.40	 -0.60|  4.00
1.00	  0.00	  0.00	 -0.60	  0.40	  0.40|  4.00
0.00	  1.00	  0.00	  0.40	 -0.60	  0.40|  4.00
------	------	------	------	------	------	------
0.00	  0.00	  0.00	  3.60	  1.60	  1.60|136.00

Pivot in row -1 and column 0



### Algorithm 16.5: Simplex algorithm in two phases¶

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

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

In [22]:
x,copt,phaseTwoTableau,feasible,bounded,itersPhaseOne,itersPhaseTwo = twoPhasesSimplex(A,b,c)


Phase I

In [23]:
for k in range(len(itersPhaseOne)):
print(prettyTableau(itersPhaseOne[k][0]))

  1.00	  3.00	  0.00	  4.00	  1.00	  1.00	  0.00	  0.00|  2.00
1.00	  2.00	  0.00	 -3.00	  1.00	  0.00	  1.00	  0.00|  2.00
-1.00	 -4.00	  3.00	  0.00	  0.00	  0.00	  0.00	  1.00|  1.00
------	------	------	------	------	------	------	------	------
-1.00	 -1.00	 -3.00	 -1.00	 -2.00	  0.00	  0.00	  0.00| -5.00

1.00	  3.00	  0.00	  4.00	  1.00	  1.00	  0.00	  0.00|  2.00
0.00	 -1.00	  0.00	 -7.00	  0.00	 -1.00	  1.00	  0.00|  0.00
0.00	 -1.00	  3.00	  4.00	  1.00	  1.00	  0.00	  1.00|  3.00
------	------	------	------	------	------	------	------	------
0.00	  2.00	 -3.00	  3.00	 -1.00	  1.00	  0.00	  0.00| -3.00

1.00	  3.00	  0.00	  4.00	  1.00	  1.00	  0.00	  0.00|  2.00
0.00	 -1.00	  0.00	 -7.00	  0.00	 -1.00	  1.00	  0.00|  0.00
0.00	 -0.33	  1.00	  1.33	  0.33	  0.33	  0.00	  0.33|  1.00
------	------	------	------	------	------	------	------	------
0.00	  1.00	  0.00	  7.00	  0.00	  2.00	  0.00	  1.00|  0.00



Phase II

In [24]:
for k in range(len(itersPhaseTwo)):
print(prettyTableau(itersPhaseTwo[k][0]))

  1.00	  0.00	  0.00	-17.00	  1.00|  2.00
-0.00	  1.00	 -0.00	  7.00	 -0.00| -0.00
0.00	  0.00	  1.00	  3.67	  0.33|  1.00
------	------	------	------	------	------
0.00	  0.00	  0.00	  3.00	 -5.00| -7.00

1.00	  0.00	  0.00	-17.00	  1.00|  2.00
0.00	  1.00	  0.00	  7.00	  0.00|  0.00
-0.33	  0.00	  1.00	  9.33	  0.00|  0.33
------	------	------	------	------	------
5.00	  0.00	  0.00	-82.00	  0.00|  3.00

1.00	  2.43	  0.00	  0.00	  1.00|  2.00
0.00	  0.14	  0.00	  1.00	  0.00|  0.00
-0.33	 -1.33	  1.00	  0.00	  0.00|  0.33
------	------	------	------	------	------
5.00	 11.71	  0.00	  0.00	  0.00|  3.00