Two-Phase Simplex Algorithm - Python Code

Today I post code that carries out the well known Two-Phase (or Artificial Variable) Simplex algorithm presented by Papadimitriou and Steiglitz. This method can be used to define an initial basic feasible solution (bfs) in a linear program where not all constraints are given by less-than linear inequalities but instead some linear equations are present. If all constraints were given by linear inequalities we could just introduce a slack variable per inequality, but if some linear equations are present we can use the Two-Phase Simplex algorithm.


-------------------------------------------------------------

def printTableu(tableu):
print '----------------------'
for row in tableu:
print row
print '----------------------'
return


def pivotOn(tableu, row, col):
j = 0
pivot = tableu[row][col]
for x in tableu[row]:
tableu[row][j] = tableu[row][j] / pivot
j += 1
i = 0
for xi in tableu:
if i != row:
ratio = xi[col]
j = 0
for xij in xi:
xij -= ratio * tableu[row][j]
tableu[i][j] = xij
j += 1
i += 1
return tableu


# assuming tablue in standard form with basis formed in last m columns
def phase_1_simplex(tableu):

THETA_INFINITE = -1
opt = False
unbounded = False
n = len(tableu[0])
m = len(tableu) - 2

while ((not opt) and (not unbounded)):
min = 0.0
pivotCol = j = 1
while(j < (n-m)):
cj = tableu[1][j]
if (cj < min):
min = cj
pivotCol = j
j += 1
if min == 0.0:
opt = True
continue
pivotRow = i = 0
minTheta = THETA_INFINITE
for xi in tableu:
if (i > 1):
xij = xi[pivotCol]
if xij > 0:
theta = (xi[0] / xij)
if (theta < minTheta) or (minTheta == THETA_INFINITE):
minTheta = theta
pivotRow = i
i += 1
if minTheta == THETA_INFINITE:
unbounded = True
continue
tableu = pivotOn(tableu, pivotRow, pivotCol)
return tableu

def simplex(tableu):

THETA_INFINITE = -1
opt = False
unbounded = False
n = len(tableu[0])
m = len(tableu) - 1

while ((not opt) and (not unbounded)):
min = 0.0
pivotCol = j = 0
while(j < (n-m)):
cj = tableu[0][j]
if (cj < min) and (j > 0):
min = cj
pivotCol = j
j += 1
if min == 0.0:
opt = True
continue
pivotRow = i = 0
minTheta = THETA_INFINITE
for xi in tableu:
if (i > 0):
xij = xi[pivotCol]
if xij > 0:
theta = (xi[0] / xij)
if (theta < minTheta) or (minTheta == THETA_INFINITE):
minTheta = theta
pivotRow = i
i += 1
if minTheta == THETA_INFINITE:
unbounded = True
continue
tableu = pivotOn(tableu, pivotRow, pivotCol)
return tableu

def drive_out_artificial_basis(tableu):
n = len(tableu[0])
j = n - 1
isbasis = True
while(j > 0):
found = False
i = -1
row = 0
for xi in tableu:
i += 1
if (xi[j] == 1):
if (found):
isbasis = False
continue
elif (i > 1):
row = i
found = True
elif (xi[0] != 0):
isbasis = False
continue
if (isbasis and found):
if (j >= n):
tableu = pivotOn(tableu, row, j)
else:
return tableu
j -= 1
return tableu

def two_phase_simpelx(tableu):
infeasible = False
tableu = phase_1_simplex(tableu)
sigma = tableu[1][0]
if (sigma > 0):
infeasible = True
print 'infeasible'
else:
#sigma is equals to zero
tableu = drive_out_artificial_basis(tableu)
m = len(tableu) - 2
n = len(tableu[0])
n -= m
tableu.pop(1)
i = 0
while (i < len(tableu)):
tableu[i] = tableu[i][:n]
i += 1
tableu = simplex(tableu)
return tableu

def getTableu(c, eqs, b):
#assume b >= 0 so if there is any b[i] negative make sure to enter
#it possitive by multiplying (-1 * eqs[i]) and (-1 * b[i]) for all i
tableu = []
m = len(eqs)
n = len(c)
c.insert(0, 0.0)
artificial = []
sigma = [0.0]
i = 0
while (i < n):
sigma.append(0.0)
i += 1
i = 0
while (i < m):
artificial.append(0.0)
sigma.append(1.0)
i += 1
c.extend(artificial)
tableu.append(c)
tableu.append(sigma)
i = 0
for eq in eqs:
eq.insert(0, b[i])
eq.extend(artificial)
eq[n+1+i] = 1.0
tableu.append(eq)
i += 1
i = 0
for xi in tableu:
if (i > 1):
j = 0
for xij in xi:
tableu[1][j] -= xij
j += 1
i += 1
return tableu

-------------------------------------------------------------
Given a Linear Program of the form:
Minimize:      c·x
subject to:    Ax≤b


the function getTableu takes the cost vector c, the matrix A, and the vector b. You must make sure to do the necessary operations to ensure that b≥0.


So getTableu will return the tableu ready for the function two_phase_simpelx to do the work and return the final tableu from which we can get the information about the optimal bfs.


So I ran the following example:
Minimize: z = x1 + x2 + x3 + x4 + x5


Subject to: 3x1 + 2x2 + x3               = 1
                 5x1 +  x2 +  x3 + x4        = 3
                 2x1 + 5x2 + x3 +       x5  = 4
as:
----------------------------------------------------------

c = [ 1.0, 1.0, 1.0, 1.0, 1.0,]

eq1 = [ 3.0 ,  2.0 , 1.0 ,  0.0,  0.0]

eq2 = [ 5.0 ,  1.0 , 1.0 ,  1.0,  0.0]

eq3 = [ 2.0 ,  5.0 , 1.0 ,  0.0,  1.0]

b = [1.0 , 3.0 , 4.0]

eqs = []

eqs.append(eq1)

eqs.append(eq2)

eqs.append(eq3)

tableu = getTableu(c,eqs,b)

printTableu(tableu)

tableu = two_phase_simpelx(tableu)

printTableu(tableu)

print 'minimum cost is = {}'.format( -tableu[0][0])

------------------------------------------------------

Here I end up printing the original tableu used by the phase 1 of he algorithm, and then the final tableu produced by the simplex. I also print information telling us if we discovered an optimal solution, or the polytope is unbounded or infeasible.

MArio

Simplex Method in Python

Well today I am going to post some code that carries the Simplex Algorithm in Python.
I have talked about the simplex algorithm before so I am not going to talk much about it. This is the first piece of code I ever write with Python so excuse my style. I did find python to be a rather easy language and I am planning on learning more about it.


Well, the simplex method I present here today takes a tableu in standard form with slack variables already introduced. So here we are assuming that we are solving a set of less-than linear inequalities and we have created a tableu with slack variables already introduced.
See below for example:


let this LP:


Max: 0.5*x + 3*y + z + 4*w
Subject to:  
  x + y + z + w ≤ 40
-2x - y + z + w ≤ 10
             y - w ≤ 10


the corresponding tableu would look like:

[ 0.0 , -0.5 , -3.0 ,-1.0 , -4.0,  0.0,  0.0,  0.0]
[40.0 ,  1.0 ,  1.0 , 1.0 ,  1.0,  1.0,  0.0,  0.0]
[10.0 , -2.0 , -1.0 , 1.0 ,  1.0,  0.0,  1.0,  0.0]
[10.0 ,  0.0 ,  1.0 , 0.0 , -1.0,  0.0,  0.0,  1.0]



all right, so here the code: Python
------------------------------------------------------

def printTableu(tableu):
print '----------------------'
for row in tableu:
print row
print '----------------------'
return


def pivotOn(tableu, row, col):
j = 0
pivot = tableu[row][col]
for x in tableu[row]:
tableu[row][j] = tableu[row][j] / pivot
j += 1
i = 0
for xi in tableu:
if i != row:
ratio = xi[col]
j = 0
for xij in xi:
xij -= ratio * tableu[row][j]
tableu[i][j] = xij
j += 1
i += 1
return tableu


# assuming tablue in standard form with basis formed in last m columns
def simplex(tableu):

THETA_INFINITE = -1
opt = False
unbounded = False
n = len(tableu[0])
m = len(tableu) - 1

while ((not opt) and (not unbounded)):
min = 0.0
pivotCol = j = 0
while(j < (n-m)):
cj = tableu[0][j]
# we will simply choose the most negative column
#which is called: "the nonbasic gradient method"
#other methods as "all-variable method" could be used
#but the nonbacis gradient method is simpler
#and all-variable method has only been shown to be superior for some
#extensive experiments by Kuhn and Quandt, the random tests used
#by Kuhn and Quandt might not really represent "typical" LP's for
#certain users.
if (cj < min) and (j > 0):
min = cj
pivotCol = j
j += 1
if min == 0.0:
#we cannot profitably bring a column into the basis
#which means that we've found an optimal solution
opt = True
continue
pivotRow = i = 0
minTheta = THETA_INFINITE
for xi in tableu:
# Bland's anticycling algorithm  is theoretically a better option than
#this implementation because it is guaranteed finite while this policy can produce cycling.
#Kotiath and Steinberg (1978) reported that cylcing does occur in practice
#contradicting previous reports. For simplicity, I don't use Bland's algorithm here
#so I just choose smallest xij for pivot row
if (i > 0):
xij = xi[pivotCol]
if xij > 0:
theta = (xi[0] / xij)
if (theta < minTheta) or (minTheta == THETA_INFINITE):
minTheta = theta
pivotRow = i
i += 1
if minTheta == THETA_INFINITE:
#bringing pivotCol into the basis
#we can move through that vector indefinetly without
#becoming unfesible so the polytope is not bounded in all directions
unbounded = True
continue

#now we pivot on pivotRow and pivotCol
tableu = pivotOn(tableu, pivotRow, pivotCol)
print 'opt = {}'.format(opt)
print 'unbounded = {}'.format(unbounded)
print 'Final Tableu'
printTableu(tableu)
return tableu

------------------------------------------------------
To test it let's use the same tableu we showed above:

------------------------------------------------------

z  = [ 0.0 , -0.5 , -3.0 ,-1.0 , -4.0,  0.0,  0.0,  0.0]

x1 = [40.0 ,  1.0 ,  1.0 , 1.0 ,  1.0,  1.0,  0.0,  0.0]

x2 = [10.0 , -2.0 , -1.0 , 1.0 ,  1.0,  0.0,  1.0,  0.0]

x3 = [10.0 ,  0.0 ,  1.0 , 0.0 , -1.0,  0.0,  0.0,  1.0]

tableu = []

tableu.append(z)

tableu.append(x1)

tableu.append(x2)

tableu.append(x3)

tableu = simplex(tableu)

------------------------------------------------------

simplex returns the final tableu from where we can retrieve the optimal basic feasible solution and its corresponding global optimal cost found at tableu[0][0].

Coming up soon a two-phase simplex algorithm that can help us when the linear program has some strict equations instead of only less than inequalities. Having only less-than-inequalities makes it easy to find a basis for the LP because of the introduced slack variables. When we have equalities too it becomes a little bit more difficult. But with the two-phase simplex algorithm we can take care of that.

MArio

Maximum Subsequence Sum Problem. C++

Today I am loading some code to solve a simple problem called "The Maximum Subsequence Sum Problem". OK, maybe the "official" name is different, but you get the concept. The problem is described as follows:


Problem:
Let us have an sequence/array of numbers [-1 , -2, 6 , 0, -4, 10, -3, -13, 1]. We want to find a subsequence of this array such that the sum of the numbers contained in the subsequence is no smaller than any other subsequence of numbers, meaning its sum is the maximum sum attainable by any subsequence in the array. For example, in the case above, the maximum subsequence is: [-1 , -2, (6 , 0, -4, 10), -3, -13, 1]  , from the third to the sixth items for a total sum of 12.


I want you to notice that if all numbers in the array are positive then obviously the maximum subsequence is the whole array itself. If all numbers are negative then an empty subsequence will have sum of zero and would be maximum. But suppose that you are not allowed to select an empty subsequence. So if all subsequences have negative sum then you have to select the one with the less negative sum. Well, today I am pasting code here that solves this problem.


Here you'll see three basic algorithms. One is the Brute Force algorithm that works in O(N^2) time, one is the Divide and Conquer algorithm that works in O(N logN), and the linear time algorithm that of course runs in O(N) time.


Here they go:


Brute Force:



int BruteForce(vector<int> A, int & start, int & end)
{
int max = 0;
bool first = true;
for(int i = 0; i < A.size(); i++)
{
int thisSubVal = 0;
for(int j = i; j < A.size(); j++)
{
thisSubVal += A[j];
if(thisSubVal > max || first)
{
max = thisSubVal;
start = i;
end = j;
first = false;
}
}
}
return max;
}

------------------------------


Divide and Conquer:



vector<int> FIND_MAX_CROSSING_SUBARRAY(vector<int> & A, int low, int mid, int high)
{
vector<int> tuple(3, 1);
int MAX_LEFT = 0;
int MAX_RIGHT = 1;
int MAX_SUM = 2;

int left_sum = 0;
int sum = 0;
bool first = true;
for(int i = mid; i >= low; i--)
{
sum += A[i];
if(sum > left_sum || first)
{
left_sum = sum;
tuple[MAX_LEFT] = i;
first = false;
}
}


int right_sum = 0;
sum = 0;
first = true;
for(int i = mid + 1; i <= high; i++)
{
sum += A[i];
if(sum > right_sum || first)
{
right_sum = sum;
tuple[MAX_RIGHT] = i;
first = false;
}
}
tuple[MAX_SUM] = left_sum + right_sum;
return tuple;
}


vector<int> FIND_MAXIMUM_SUBARRAY(vector<int> & A, int low, int high)
{
vector<int> tuple(3, 1);
int LEFT = 0;
int RIGHT = 1;
int MAX_SUM = 2;


if(low == high)
{
tuple[LEFT] = low;
tuple[RIGHT] = high;
tuple[MAX_SUM] = A[low];
return tuple;
}
else
{
int mid = (low + high) / 2;
vector<int> tupleLeft = FIND_MAXIMUM_SUBARRAY(A, low, mid);
vector<int> tupleRight = FIND_MAXIMUM_SUBARRAY(A, mid + 1, high);
vector<int> tupleCross = FIND_MAX_CROSSING_SUBARRAY(A, low, mid, high);


if(tupleLeft[MAX_SUM] >= tupleRight[MAX_SUM] && tupleLeft[MAX_SUM] >= tupleCross[MAX_SUM])
return tupleLeft;
else if(tupleRight[MAX_SUM] >= tupleLeft[MAX_SUM] && tupleRight[MAX_SUM] >= tupleCross[MAX_SUM])
return tupleRight;
else
return tupleCross;
}
}


int DivideAndConquer(vector<int> A, int & start, int & end)
{
vector<int> tuple = FIND_MAXIMUM_SUBARRAY(A, 0, A.size() - 1);
start = tuple[0];
end = tuple[1];
return tuple[2];
}

------------------------------
Linear:



int Linear(vector<int> A, int & start, int & end)
{
int size = A.size();
vector<int> P(size, 0);
vector<int> P_start(size, 0);
vector<int> P_end(size, 0);
P_start[0] = 0;
P_end[0] = 0;
int max = P[0] = A[0];
start = end = 0;
int sum = 0;
for(int j = 1; j < size; j++)
{
if(A[j] >= P[j-1] + A[j])
{
P[j] = A[j];
P_start[j] = j;
P_end[j] = j;
}
else
{
P[j] = P[j-1] + A[j];
P_start[j] = P_start[j-1];
P_end[j] = j;
}
if(P[j] >= max)
{
max = P[j];
start = P_start[j];
end = P_end[j];
}
}
return max;
}

------------------------------


So for all these functions just pass start and end as references and pass the array as an integer vector.


MArio