Instructions
Requirements and Specifications
Source Code
import numpy as np
import re
def read_data(file_name):
"""
This function read the coefficients for a LP problem from a file given by file_name
:param file_name: Name of the file
:return:
x: list with the coefficients of the objective function
c: list of lists with the coefficients of the constraints
b: values of the constants for the constraints
"""
# First, read the file
n_variables = -1
n_constraints = -1
constraints = []
c = []
b = []
r = re.compile(r'([^\t]*)\t*')
with open(file_name, 'r') as f:
lines = f.readlines()
# The first line contains the objective function
obj = lines[0].strip()
objs_coeffs = obj.split()
x = [float(x) for x in objs_coeffs]
# The rest are the constraints
for i in range(1, len(lines)):
constraint = lines[i].strip()
constraints_coeffs = constraint.split()
ci = [float(x) for x in constraints_coeffs[:-1]]
bi = float(constraints_coeffs[-1])
c.append(ci)
b.append(bi)
n_variables = len(objs_coeffs)
n_constraints = len(c) + n_variables # we add n_variables for the constraints of type x_i >= 0
return x, c, b
def getTableau(x, c, b):
"""
This function builds the Tableay for the Simplex Method
:param x: list with the coefficients of the objective function
:param c: list of lists with the coefficients of the constraints
:param b: values of the constants for the constraints
return: 2-D NumPy Array (matrix)
"""
n = len(c) + 1 # the number of rows is equal to the number of constraints plus one for the objective
m = 2*len(x) +1 # the number of coolumns is the double of the number of variable (because of slcak variables) plus one for the equality
nvars = len(x) #number of real variables
nconst = len(c) # number of constraints
# Stack the values of the objective
T = np.array(c) # initializw the table with the values of the coefficients of the constraints
T = np.hstack((T, np.eye(nconst))) # coefficients for slack variables
# Now stack b
T = np.column_stack((T, np.array([b]).T))
# Finally, stack the row of the objective
obj_row = -np.array(x)
# To this row, add the indexes for slack variables
obj_row = np.hstack((obj_row, np.zeros(nconst+1)))
obj_row[-1]= 1
# Now stack
T = np.row_stack((T, obj_row))
return T
def isUnbounded(T):
"""
This function checks if a Tableau T for Simplex Method is unbounded
:param T: Tableau (2D NumPy Array)
:return: boolean
"""
# This function checks the Tableau and checks if the problem is unbounded
# Examine each column
(n, r) = T.shape
for j in range(r-1): # iterate through columns except the last one
column_ok = False
for i in range(n):
if T[i,j] <= 0:
column_ok = True
break
if not column_ok: # column has only positive values, so its unbounded
return True
return False
def getPivotColumnIndex(T):
"""
This functions get the index of the pivot column in T
:param T: 2D NumPy Array
:return: integer
"""
# Initialize with the index and minimum value in the first column
index = 0
min_val = min(T[:,0])
(n, r) = T.shape
# Iterate through columns and get the one with the most negative value
for i in range(1,r):
val =T[-1,i]
if val < min_val:
min_val = val
index = i
return index
def getPivotRowIndex(T, pc):
"""
This function gets the index of the pivt row
:param T: 2D NumPy Array
:param pc: index of the pivot column
:return: integer
"""
(n,r) = T.shape
# Now, divide the pivot column by b and select the index of the min val
pivot_col = T[:,pc]
b = T[:,-1]
# Divide b by the coefficients in the pivot column
result = np.zeros(len(b))
for i in range(len(b)):
if pivot_col[i] > 0:
result[i] = b[i]/pivot_col[i]
else:
result[i] = np.inf
#result = np.divide(b, pivot_col)
# Pick the index of the minimum value
index = np.argmin(result[:-1])
return index
def printTable(T):
"""
this function prints a 2D NumPy Array in a way easier to look for the user
:param T: 2D NumPy Array
:return: None
"""
(n,r) = T.shape
for i in range(n):
print("[ ", end="")
for j in range(r):
print("{:.2f}".format(T[i,j]), end= "\t")
print("]")
if __name__ == '__main__':
import sys
assert len(sys.argv) > 1, "No arguments passed to script!"
file_name = sys.argv[1]
"""
SIMPLEX METHOD
"""
# Read data from file
x, c, b = read_data(file_name)
# Create Table
T = getTableau(x, c, b)
# Number of variables
nvars = len(x)
# Define a max number of iterations
max_iters = 500 # define a max number of iterations
# Shape of table
(n, r) = T.shape
# Iteration counter
iter = 1
# List used to store the index of the columns pivoted
pivots = []
solution_flag = 0 # 0 for no solution found, 1 for solution found, 2 for unbounded problem and 3 for infeasible
# The followig lines prints the Tableau at the initial status
"""
print(f"\n\nInitial table")
printTable(T)
print("\n\n")
"""
# Start algorithm
variables_row = -1*np.ones(T.shape[1])
while True:
# Check unbounded
if isUnbounded(T):
solution_flag = 2 # unbounded
break
# Get pivot column
pc = getPivotColumnIndex(T)
"""
# if the pivot column is for one of the slack variables, then the method is infeasible
if pc > nvars-1:
solution_flag = 3
break
"""
# If the pivot column has only negative values, the problem does not have a solution
if len(np.where(T[:,pc] > 0)[0]) == 0: # no positive values, only negative
solution_flag = 2
break
# Get pivot row
pr = getPivotRowIndex(T, pc)
# Pivot element
pe = T[pr, pc]
# Update pivot row
T[pr, :] = np.divide(T[pr, :], pe)
for i in range(n) : # iterate through rows
if i != pr: # not the pivot row
pivot_col_coeff = T[i, pc]
new_row = T[i,:] - pivot_col_coeff*T[pr,:]
T[i,:] = new_row
# use the following lines to see the Tableau at each iteration
"""
print(f"\n\nITER {iter}")
printTable(T)
print("\n\n")
"""
pivots.append([pr, pc])
variables_row[pc] = pr
# If the last row of the Table has only positive values, it means that the method finished
last_row = T[-1,:]
if np.all((last_row >= 0)):
solution_flag = 1
break
iter += 1
if max_iters == iter: # max number of iterations reached. No solution found. Problem should be infeasible
solution_flag = 3
break
if solution_flag == 1:
#printTable(T)
# Solution found
X = np.zeros(nvars) # Array to store the solutions
vars_found = [] # List to store the index of the variables found in the method
for i in range(nvars):
idx = int(variables_row[i])
if not idx in vars_found and idx < nvars:
X[i] = T[idx, -1]
vars_found.append(idx)
"""
for pivot in pivots:
X[pivot[1]] = T[pivot[0], -1]
vars_found.append(pivot[1])
"""
Z = T[-1, -1] # Value of the objective function
# Now, calculate the value of the variable not found directly.
s = 0
for i in vars_found:
s += x[i] * X[i]
# Calculate that variable from the objective
for i in range(nvars):
if not i in vars_found:
X[i] = Z - s
# Print
print("optimal")
print(Z)
for x in X:
print(x, end=" ")
elif solution_flag == 2:
print("unbounded")
elif solution_flag == 3:
print("infeasible")