## Instructions

**Objective**

## 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")