"""
Author: Niklaus Eggenberg
Created: 23.02.2023
Simplex Class - contains the code to solve an LP via simplex method
"""
# run command:
# cd
# python c:\HEPIA\neg-teaching-material\Simplexe\main.py C:\HEPIA\neg-teaching-material\2022-2023\2e\lp_test.txt
import numpy as np
import os.path
from os import path
from src.constants import Constants, OptimizationType, OptStatus
import time
import sys
class Simplexe:
def __init__(self, A, B, C, objType):
self.IsPhaseI = False
self.__isRemovingAuxVariables = False
self.__PivotCount = 0
self.__PhaseIPivotCount = 0
self.__PrintDetails = False
self.Costs = C
self.AMatrix = A
self.RHS = B
self.ObjectiveType = objType
self.OptStatus = OptStatus.Unknown
# some calculated stuff
self.NumRows = self.AMatrix.shape[0]
self.NumCols = self.AMatrix.shape[1]
self.hasCycle = False
# self.initBasicVariable = self.__basicVariables.copy()
def SolveProblem(self):
self.__start = time.time()
self.__initTableau()
# make sure we have a feasible solution - go to PhaseI
# phase I runs only if
if (self.__isSolutionFeasible() != True):
print("Initial solution is unfeasible - starting Phase I...")
self.__solvePhaseI()
if self.OptStatus == OptStatus.Infeasible:
self.__end = time.time()
print(
"-------------------------------------------------------------------------------")
print(
"Pase I has non-zero optimum : STOP - problem has not feasible solution!")
print(
"-------------------------------------------------------------------------------")
print("Execution statistics")
self.PrintSolution()
print("------------------------------------------------")
print("LP is INFEASIBLE - no solution exists!")
print("------------------------------------------------")
sys.exit()
else:
self.PrintTableau("Initial feasible solution")
else:
self.OptStatus = OptStatus.Feasible
self.__performPivots()
#
self.__end = time.time()
if self.OptStatus == OptStatus.NotBounded:
print(
"-------------------------------------------------------------------------------")
print("Found unbounded column => solution has no optimal solution")
print(
"-------------------------------------------------------------------------------")
print("Execution statistics")
self.PrintSolution()
if self.OptStatus == OptStatus.NotBounded:
print("------------------------------------------------")
print("LP is UNBOUNDED!")
print("------------------------------------------------")
def __performPivots(self):
# we DO have a feasible solution - go on with the simplex
while (self.OptStatus == OptStatus.Feasible):
self.__pivotTableau(self.__selectPivot())
# self.PrintTableau("After pivot")
def PrintTableau(self, header):
print("------------------------------------------------")
print("Simplex tableau: ", header, " Opt status: ", self.OptStatus)
print("------------------------------------------------")
# make sure to print all in one line - note that the width of the command prompt might wrap automatically
# but output to file works just perfect!
varNames = [self.__varName(item) for item in self.__basicVariables]
varNames.append(self.__padStr("* OBJ *"))
tmpArray = np.array([varNames]).T
tableau = np.append(tmpArray, self.__tableau.copy(), axis=1)
with np.printoptions(precision=3, suppress=True, linewidth=np.inf):
print(tableau)
def ObjValue(self):
objVal = self.__tableau[-1, -1]
# phase I has no associated LP => we always have a min !
if self.IsPhaseI:
return -objVal
return objVal if self.ObjectiveType == OptimizationType.Max else -objVal
def PrintSolution(self):
if self.OptStatus == OptStatus.Optimal:
if self.__PrintDetails:
self.PrintTableau("FINAL TABLEAU")
print("------------------------------------------------")
print("Non-zero variables")
nonZeros = []
print(self.__basicVariables)
for rowId, baseColId in enumerate(self.__basicVariables):
# NOTE: in Phase 1, the LAST row is original objective => skip it:
if (not self.IsPhaseI) & (baseColId < 0) & (rowId < len(self.__basicVariables) - 1):
print("ERROR: basic variable with negative index at position ",
rowId, " value " + baseColId)
sys.exit()
if baseColId >= 0:
nonZeros.append([baseColId, "{} = {}".format(self.__varName(
baseColId), "%.4f" % self.__getBasicVariableValue(rowId, baseColId))])
# sort array
nonZeros = sorted(nonZeros, key=lambda tup: tup[0])
for val in nonZeros:
print(val[1])
else:
print("------------------------------------------------")
print("No solution found...")
#
print("------------------------------------------------")
print("Num rows / cols : ", self.NumRows, " / ", self.NumCols)
print("Tableau dimensions: : ", self.__tableau.shape)
print("Optimiser status : ", self.OptStatus)
print("Objective Value : ", self.ObjValue())
print("Nbr pivots Phase I : ", self.__PhaseIPivotCount)
print("Nbr pivots Phase II : ", self.__PivotCount)
print("Total exec time [s] : ", self.__end - self.__start)
print("------------------------------------------------")
def __padStr(self, str):
return str.ljust(8, " ")
def __varName(self, colId):
if colId < 0:
return "* OBJ_0 *"
if colId < self.NumCols:
return self.__padStr("x[{}]".format(colId))
return self.__padStr("z[{}]".format(colId - self.NumCols))
def __initTableau(self):
# initializes the tableau from the previously loaded LP
# tableau is
# -------------
# | A | I | b |
# |-----------|
# | c | 0 | 0 |
# -------------
# NOTE: to make things easier for Phase I and signs, always make sure b is >= 0 (if originally we have b < 0, simply change the signs of the entire row)
# make sure to perform a deep copy of the coefficients - also make sure all are 2D arrays (i.d. add arrays as [])
tmpA = self.AMatrix.copy()
tmpI = np.identity(self.NumRows)
tmpB = np.array([self.RHS.copy()]).T
# now append an identity matrix to the right and the b vector
self.__tableau = np.append(tmpA, tmpI, axis=1)
self.__tableau = np.append(self.__tableau, tmpB, axis=1)
# last row
tmpC = np.array(
[np.concatenate((self.Costs.copy(), np.zeros(self.NumRows + 1)))])
self.__tableau = np.append(self.__tableau, tmpC, axis=0)
self.__basicVariables = np.arange(
self.NumCols, self.NumCols + self.NumRows, 1, dtype=int)
self.__initBasicVariable = self.__basicVariables.copy()
self.TableauRowCount, self.TableauColCount = self.__tableau.shape
for rowId in range(self.NumRows):
if self.__tableau[rowId, -1] < -Constants.EPS:
self.__tableau[rowId, :] = self.__tableau[rowId, :] * -1.
if self.__PrintDetails:
self.PrintTableau("Initial tableau")
# returns next pivot as an array of leaving row id and entering col Id
# return None if there is pivot (sets optimisation status accoringly before returning)
def __selectPivot(self):
# print(self.__basicVariables)
colId = self.__selectEnteringColumn()
if ((not colId) | colId < 0):
# no more entering column => optimiser status is OPTIMAL (must be - we never select a pivot if the tableau is unfeasible)
self.OptStatus = OptStatus.Optimal
return None
# get leaving row id
rowId = self.__selectLeavingRow(colId)
if (rowId < 0):
# there is a negative reduced cost column that never hits any constraint => problem is NotBounded!
self.OptStatus = OptStatus.NotBounded
return None
# normal case : we do have a standard pivot!
# print(f"Before : {self.__basicVariables}")
# print(f"RowID = {rowId}")
# idx_to_replace = np.argwhere(self.__basicVariables == rowId)[0][0]
self.__basicVariables[rowId] = colId
# print(f"After : {self.__basicVariables}")
# print(idx_to_replace)
return [rowId, colId]
def __getBasicVariableValue(self, rowId, baseColId):
# the value of the basic variable is RHS / coeff in basic variable's colum
return self.__tableau[rowId, -1] / self.__tableau[rowId, baseColId]
def __solvePhaseI(self):
# initializes and solves the Phase I simplexe automatically
phaseISplx = Simplexe(self.AMatrix, self.RHS,
self.Costs, self.ObjectiveType)
phaseISplx.initPhaseI(self)
# if optimal solution of Phase I is NOT 0 => problem has no feasible solution
if phaseISplx.ObjValue() > Constants.EPS:
self.OptStatus = OptStatus.Infeasible
return
# make sure we eliminate ALL auxiliary variables from the current base (auxiliary vars are determined by their id : the all have id >= num cols + num rows)
if self.__PrintDetails:
phaseISplx.PrintTableau("Optimal solution Phase I")
phaseISplx.RemoveAuxiliaryBasicVariables()
phaseISplx.PrintSolution()
# we have a feasible solution => extract feasible solution to proceed with Phase II
# -> first part of the tableau (original matrix)
self.__tableau = phaseISplx.__tableau[0:self.NumRows,
0:self.NumCols+self.NumRows]
# -> append the objective row - previous last row of PhaseI tableau
self.__tableau = np.append(self.__tableau, np.array(
[phaseISplx.__tableau[-2, 0:self.NumCols+self.NumRows]]), axis=0)
# -> append the RHS for all except the last row
self.__tableau = np.append(self.__tableau, np.array(
[phaseISplx.__tableau[0:-1, -1]]).T, axis=1)
# finally the basic variables indices
self.__basicVariables = phaseISplx.__basicVariables[0:-1]
self.__PhaseIPivotCount = phaseISplx.__PivotCount
self.OptStatus = OptStatus.Feasible
def RemoveAuxiliaryBasicVariables(self):
# after resolving PhaseI, we might still have some auxiliary varaibles in the current basis
# if so, make sure we remove them - they always have index >= initial number of columns (i.e. self.NumCols here)
maxId = np.argmax(self.__basicVariables)
while (self.__basicVariables[maxId] >= self.NumCols):
# leaving row is the one corresponding to the auxiliary variable, i.e. row with index maxId
# we now have to retrieve which is the ORIGINAL SLACK COLUMN corresponding to the SAME row as the AUXILIARY variable (there is always one)
# in Phase I, we have 1 additional row (corresponding to the original objective row) => the id of the original SLACK is id of AUXILIARY - (numRows - 1)
originSlackVarId = self.__basicVariables[maxId] - (
self.NumRows - 1)
# pivot ROW ID is current index of auxiliary variable (i.e. maxId), and entering column is the original slack variable, i.e. originSlackVarId
# - perform pivot
self.__isRemovingAuxVariables = True
self.__pivotTableau([maxId, originSlackVarId])
self.__isRemovingAuxVariables = False
# update the max Id
maxId = np.argmax(self.__basicVariables)
# stop measuring time
self.__end = time.time()
###################################################################
########## METHODS TO BE IMPLEMENTED ##########
###################################################################
# returns cheks if the current solution is feasible or not (return True if feasible)
def __isSolutionFeasible(self):
# we MUST have that all BASIC variables are >= 0 to have a FEASIBLE solution
# to iterate over basic variables do:
for rowId, baseColId in enumerate(self.__basicVariables):
if self.__getBasicVariableValue(rowId, baseColId) < 0:
return False
return True
# returns index of entering col Id
# return None or -1 if there is pivot
def __selectEnteringColumn(self):
# find an entering column for pivot => return -1 if none is found!
# print(self.__tableau[0])
# for rowId, baseColId in enumerate(self.__basicVariables):
# print(self.__tableau[rowId, baseColId])
#
if self.__PivotCount != 0 and np.array_equal(self.__basicVariables,
self.__initBasicVariable):
self.hasCycle = True
if self.hasCycle:
for idx, val in enumerate(self.__tableau[-1, :-1]):
if val < 0:
return idx
return -1
else:
if np.size(self.__tableau[-1, :-1][np.where(self.__tableau[-1, :-1]
< 0)]) == 0:
return -1
return np.argmin(self.__tableau[-1, :-1])
# returns leaving row ID
# return None or -1 if there is pivot (sets optimisation status accoringly before returning)
def __selectLeavingRow(self, pivotColId):
# given the entering column, find the leaving row - return the index of the row or -1 if none is found
# iterate using :
# rowCount = self.TableauRowCount - 2 if self.IsPhaseI else self.TableauRowCount - 1
# for index in range(rowCount):
# rowCount = self.TableauRowCount - 2 if self.IsPhaseI else \
# self.TableauRowCount - 1
#
# lhs = self.__tableau[:rowCount, pivotColId]
# rhs = self.__tableau[:rowCount, -1]
#
# ratios = rhs / lhs
#
# candidates = np.where(ratios > 0, ratios, np.inf)
#
# if len((np.unique(candidates))) == 0:
# return -1
#
# return candidates.argmin()
minRatio = float('inf')
leavingRow = -1
rowCount = self.TableauRowCount - 2 if self.IsPhaseI else \
self.TableauRowCount - 1
for index in range(rowCount):
# Si le ration est positif
if self.__tableau[index][pivotColId] > 0:
ratio = self.__tableau[index][-1] / \
self.__tableau[index][pivotColId]
if ratio < minRatio:
minRatio = ratio
leavingRow = index
# print("leaving row :", leavingRow) # vérification ok
return leavingRow
def __pivotTableau(self, pivotIDs):
if pivotIDs is None or pivotIDs[0] < 0 or pivotIDs[1] < 0:
# no pivot => optimiser status is updated in pivot selection => return!
return
# ---- update stats
self.__PivotCount += 1
#
pivotRowId = pivotIDs[0]
pivotColId = pivotIDs[1]
# print("Pivot is ", pivotRowId, " ", pivotColId)
pivotVal = self.__tableau[pivotRowId, pivotColId]
# double-check pivot is actually positive
if pivotVal < Constants.EPS:
# pivoting around a negative value is allowed ONLY if we're doing it while removing auxiliary variables from PhaseI
if not self.__isRemovingAuxVariables:
print(
"ERROR: pivot on non-positive value {} at [{},{}]".format(pivotVal, pivotRowId, pivotColId))
self.OptStatus = OptStatus.ERROR
sys.exit()
# perform pivot (above code ensures the pivot is possible!)
# iterate using
# pivotRow = self.__tableau[pivotRowId, :] / pivotVal
# for rowId in range(self.TableauRowCount):
# Calculer la ligne du pivot normalisée (avoir la valeur 1 à la colonne du pivot)
pivotRow = self.__tableau[pivotRowId, :] / pivotVal
# Pour chaque ligne autre que le pivot
for rowId in range(self.TableauRowCount):
if rowId != pivotRowId:
# Mettre la colonne du pivot à 0
multiplier = self.__tableau[rowId, pivotColId]
self.__tableau[rowId, :] -= multiplier * pivotRow
self.__tableau[pivotRowId, :] = pivotRow
# raise Exception(
# 'Not implemented', 'Simplexe.__pivotTableau: missing code to be implemented.')
# For a given unfeasible initial problem, this method creates the corresponding Tableau of the Phase I, then launches the optimization of Phase I
def initPhaseI(self, simplexe):
# create augmented tableau for Phase I => current
# -- Original ----- # -- Auxiliary -----
# | A | I | b | | A | I | I | b |
# |-----------| |---------------|
# | c | 0 | 0 | | c | 0 | 0 | 0 |
# |-----------| |---------------|
# | d |-1 | 0 | S |
# d = - sum(coeff in column) and S is the -sum vector b (NOTE: we already made sure all elements in b are non-negative!)
# store the fact we are actually using a Phase I - to make sure we don't select a pivot in the "original" objective row!
self.IsPhaseI = True
self.__start = time.time()
# new dimensions (note : all original variables are considered as "normal" variables, whereas all auxiliary ones are now slacks!)
self.NumRows = simplexe.NumRows + 1
self.NumCols = simplexe.NumCols + simplexe.NumRows
self.OptStatus = OptStatus.Feasible
# copy tableau excep LAST row and last column (we will add them last)
tmpCopy = simplexe.__tableau.copy()
self.__tableau = tmpCopy[0:-1, 0:-1]
# add the auxiliary variables and their objective row
self.__tableau = np.append(
self.__tableau, np.identity(simplexe.NumRows), axis=1)
# add the RHS column
self.__tableau = np.append(
self.__tableau, np.array([tmpCopy[0:-1, -1]]).T, axis=1)
# create the tableau for the PhaseI from above initialized arrays!!!
raise Exception('Not implemented',
'Simplexe.initPhaseI: missing code to be implemented.')
# set dimensions
self.TableauRowCount, self.TableauColCount = self.__tableau.shape
if self.__PrintDetails:
self.PrintTableau("Phase 1 Tableau")
self.__performPivots()