import numpy as np
import pandas as pd
This is Algorithm 16.4 in Bierlaire (2015) Optimization: principles and algorithms, EPFL Press.
First, the routine that performs one iteration of the simplex using the tableau.
def simplexTableau(tableau):
"""
:param tableau: the first simplex tableau
:type tableau: pandas dataframe
:return: p, q, opt, bounded where
- p is the column of the variable that must enter the basis, or None,
- q is the row of the variable that must leave the basis, or None,
- opt is True if the tableau is optimal (in this case, p and q are None)
- bounded is False if basic direction is unbounded (in this case, p and q are None)
:rtype: int, int, bool, bool
"""
mtab, ntab = tableau.shape
m = mtab - 1
n = ntab - 1
reducedCost = tableau[-1, : -1]
# Identify the negative reduced costs
negativeReducedCost = reducedCost < 0
if not negativeReducedCost.any():
# The tableau is optimal
return None, None, True, True
# In Python, True is larger than False. The next statement returns the
# index of a True entry in the array, that is the index of a negative reduced cost.
# It is the index of the variable that will enter the basis.
p = np.argmax(negativeReducedCost)
# Calculate the maximum step that can be done along the basic direction d[p]
xb = tableau[:-1, -1]
minusd = tableau[:-1, p]
steps = np.array([xb[k] / minusd[k] if minusd[k] > 0 else np.inf for k in range(m)])
q = np.argmin(steps)
step = steps[q]
if step == np.inf:
# The tableau is unbounded
return None, None, False, False
else:
return p, q, False, True
Then, the routine to pivot the tableau.
def pivoting(tableau, p, q):
"""
:param tableau: valid simplex tableau
:type tableau: numpy.array 2D
:param p: columns of the pivot
:type p: int
:param q: row of the pivot
:type q: int
:return: pivoted tableau
:rtype: numpy.array 2D
"""
m, n = tableau.shape
if q >= m:
raise Exception(f'The row of the pivot ({q}) must be between 0 and {m - 1})')
if p >= n:
raise Exception(f'The column of the pivot ({p}) must be between 0 and {n - 1})')
thepivot = tableau[q][p]
if np.abs(thepivot) < np.finfo(float).eps:
raise Exception(f'The pivot is too close to zero: {thepivot}')
thepivotrow = tableau[q, :]
newtableau = np.empty(tableau.shape)
newtableau[q, :] = tableau[q, :] / thepivot
for i in set(range(m)) - {q}:
newtableau[i, :] = tableau[i, :] - tableau[i][p] * thepivotrow / thepivot
return newtableau
Finally, we put everything together.
def simplexAlgorithmTableau(tableau):
"""
:param tableau: valid simplex tableau
:type tableau: numpy.array 2D
:return: tableau, optimal, unbounded, where tableau is the tableau from the last iteration,
optimal is True if the last tableau is optimal,
unbounded is True if the problem is unbounded.
:rtype: numpy.array 2D, bool, bool
"""
while True:
colPivot, rowPivot, optimal, bounded = simplexTableau(tableau)
if optimal:
return tableau, True, False
if not bounded:
return tableau, False, True
tableau = pivoting(tableau, colPivot, rowPivot)
The following function is designed to print the results.
def printResults(res):
lastTableau, optimal, unbounded = res
if optimal:
print('Optimal tableau')
print(pd.DataFrame(lastTableau).to_string(index = False, header = False))
elif unbounded:
print('Unbounded problem')
else:
print('Incorrect output')
Consider the following initial tableau: $$ \begin{array}{|cccccc|c|} \hline 1& 2& 2& 1& 0& 0& 20 \\ 2& 1& 2& 0& 1& 0& 20 \\ 2& 2& 1& 0& 0& 1& 20 \\ \hline -10& -12& -12& 0& 0& 0& 0 \\ \hline \end{array} $$
T = 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]])
res = simplexAlgorithmTableau(T)
printResults(res)
Optimal tableau 0.0 0.0 1.0 0.4 0.4 -0.6 4.0 1.0 0.0 0.0 -0.6 0.4 0.4 4.0 0.0 1.0 0.0 0.4 -0.6 0.4 4.0 0.0 0.0 0.0 3.6 1.6 1.6 136.0
Consider the following initial tableau: $$ \begin{array}{|rrrrrrrr|r|} \hline 1& 2& 3& 0& 1& 0& 0& 0& 3 \\ -1& 2& 6& 0& 0& 1& 0& 0& 2 \\ 0& 4& 9& 0& 0& 0& 1 & 0& 5 \\ 0& 0& 3& 1& 0& 0& 0& 1& 1 \\ \hline 0& -8& -21& -1& 0& 0& 0& 0& -11 \\ \hline \end{array} $$
T = np.array([[1, 2, 3, 0, 1, 0, 0, 0, 3],
[-1, 2, 6, 0, 0, 1, 0, 0, 2],
[0, 4, 9, 0, 0, 0, 1 , 0, 5],
[0, 0, 3, 1, 0, 0, 0, 1, 1],
[0, -8, -21, -1, 0, 0, 0, 0, -11]])
res = simplexAlgorithmTableau(T)
printResults(res)
Optimal tableau 1.0 0.0 0.0 0.500000 0.50 -0.50 0.0 0.500000 1.000000 0.0 1.0 0.0 -0.750000 0.25 0.25 0.0 -0.750000 0.500000 0.0 0.0 0.0 0.000000 -1.00 -1.00 1.0 0.000000 0.000000 0.0 0.0 1.0 0.333333 0.00 0.00 0.0 0.333333 0.333333 0.0 0.0 0.0 0.000000 2.00 2.00 0.0 1.000000 0.000000
Consider the linear optimization problem $$ \min_{x \in \mathbb{R}^2} 3 x_1 - 2 x_2 $$ subject to $$ \begin{array}{rl} x_1 &= 1, \\ x_1, x_2 & \geq 0. \end{array} $$
The problem is clearly unbounded, as $x_1$ is fixed by the constraints, and increasing $x_2$ decreases the objective function. We will verify that it is detected by the algorithm.
It is a problem in standard form with $$A = (1 \; 0), \; b = 1, \; c = \left(\begin{array}{r} 3 \\ -2 \end{array}\right). $$
We select $x_1$ to be in the basis, so that $B=1$ and $c_B=3$. The tableau is $$ \begin{array}{| c| c |} \hline B^{-1}A & B^{-1}b\\ \hline c^T -c_B^TB^{-1}A & -c_B^TB^{-1}b\\ \hline \end{array} $$ that is $$ \begin{array}{|cc| c |} \hline 1 & 0 & 1\\ \hline 0 & -2 & -3\\ \hline \end{array} $$
Note that $x_2$ cannot be in the basis. Indeed, it would mean that $B=0$, which is not valid.
We now apply the algorithm on this tableau
T = np.array([[1, 0, 1],
[0, -2, -3]])
res = simplexAlgorithmTableau(T)
printResults(res)
Unbounded problem