Code

!pip install rich

import pandas as pd
import numpy as np

from rich import print

A basic implementation of the Simplex algorithm as described in Chapter 2 of [1] for problems of the form

Roughly speaking the implementation involves iterating through a sequence of dictionaries until an optimal solution is uncovered. A dictionary is a system of linear equations associated with a feasible solution such as:

We will then use the implementation to solve

The data corresponding to this problem and other problem instances can be obtained from the methods in the collapsed code cell below.

Code

def get_problem_1():
  # objective function coefficients
  c = np.array([[5],[4],[3]])
  #c.shape

  # constraint coefficients
  A = np.array([[2,3,1],[4,1,2],[3,4,2]])
  #A.shape

  # right hand sides
  b = np.array([[5],[11],[8]])
  #b.shape

  return c, A, b

def get_problem_2():
  # objective function coefficients
  c = np.array([[5],[5],[3]])
  
  # constraint coefficients
  A = np.array([[1,3,1],[-1,0,3],[2,-1,2],[2,3,-1]])
  
  # right hand sides
  b = np.array([[3],[2],[4],[2]])
  
  return c, A, b

def get_problem_2_1a():
  # objective function coefficients
  c = np.array([[3],[2],[4]])
  
  # constraint coefficients
  A = np.array([[1,1,2],[2,0,3],[2,1,3]])
  
  # right hand sides
  b = np.array([[4],[5],[7]])
  
  return c, A, b

def get_problem_2_1b():
  # objective function coefficients
  c = np.array([[5],[6],[9],[8]])
  
  # constraint coefficients
  A = np.array([[1,2,3,1],[1,1,2,3]])
  
  # right hand sides
  b = np.array([[5],[3]])
  
  return c, A, b

def get_problem_2_1c():
  # objective function coefficients
  c = np.array([[2],[1]])
  
  # constraint coefficients
  A = np.array([[2,3],[1,5],[2,1],[4,1]])
  
  # right hand sides
  b = np.array([[3],[1],[4],[5]])
  
  return c, A, b  

def get_problem_2_2():
  # objective function coefficients
  c = np.array([[2],[3],[5],[4]])
  
  # constraint coefficients
  A = np.array([[1,2,3,1],[1,1,2,3]])
  
  # right hand sides
  b = np.array([[5],[3]])
  
  return c, A, b

def get_unbounded_problem():
  # objective function coefficients
  c = np.array([[1],[-1]])
  
  # constraint coefficients
  A = np.array([[-2,1],[-1,-2]])
  
  # right hand sides
  b = np.array([[-1],[-2]])
  
  return c, A, b  

def get_problem_1_6():
  # objective function coefficients
  c = np.array([[6],[3],[8],[3],[9],[5]])
  #c.shape

  # constraint coefficients
  A = np.array([[1,1,0,0,0,0],[0,0,1,1,0,0],[0,0,0,0,1,1],
                [1,0,1,0,1,0],[0,1,0,1,0,1]])
  #A.shape

  # right hand sides
  b = np.array([[480],[400],[230],[420],[250]])
  #b.shape

  return c, A, b

c, A, b = get_problem_1()

Constants.

m = A.shape[0] # number of constraints
n = A.shape[1] # number of decision variables

FIRST_ROW = 0
ALL_BUT_LAST_ROW = -1

FIRST_COLUMN = 0
SECOND_COLUMN = 1

Z_ROW = m

Identify the variable to enter the basis

Strictly speaking the method identifies the column corresponding to the variable that is going to be entering the basis.

Code

def get_entering_column(chvatal_dict):
  entering_column = None

  # and second column and after since we are using
  # the first column to store the RHSs.
  z_coeffs = chvatal_dict[ Z_ROW , SECOND_COLUMN : ]

  if np.max( z_coeffs ) > 0:
    # the column that enters from the dictionary
    # is one with the largest positive coefficient
    # (since we are maximizing)
    entering_column = np.argmax( z_coeffs ) + 1
    # we need to add one since we looked from the
    # second column onwards

  return entering_column

Identify the variable to leave the basis

The method identifies the row corresponding to the variable that is going to be leaving the basis.

Code

def get_leaving_row(chvatal_dict, entering_column):
  cd = chvatal_dict[ FIRST_ROW : ALL_BUT_LAST_ROW, : ]

  cd_first_column = cd[:, FIRST_COLUMN].copy()
  cd_entering_column = -cd[:, entering_column].copy()
  
  # avoid a divide by zero - identify if any
  # coefficients in the entering column are zero
  if sum(cd_entering_column==0):
    # set the numerator to -1
    cd_first_column[cd_entering_column==0] = -1
    # now replace the 0 in the denominators with 1
    cd_entering_column[cd_entering_column==0] = 1
    # This ensures that such columns are not considered
    # when it comes time to decide the leaving variable
    # Also note that after this step
    # all entries in cd_entering_column==0 will be False

  leaving_row_candidates = cd_first_column / cd_entering_column
  # print(leaving_row_candidates)

  # Identify the leaving variable
  leaving_row = None
  if np.min( leaving_row_candidates ) < 0:
    # Replace candidates violating non-negativity by infinity
    # so these will be ignored when choosing the leaving variable
    leaving_row_candidates[ leaving_row_candidates < 0 ] = np.inf

  # the row that leaves from the dictionary
  # is the one that will impose the strictest
  # constraint on the requirement that the
  # variable remain non-negative
  leaving_row = np.argmin( leaving_row_candidates )

  #print(leaving_row)
  return leaving_row

Construct the system for the next iteration

Code

def get_pivot_row(chvatal_dict, leaving_row, entering_column):
  # Need to find the coefficient of the non-basic variable
  # that is going to be entering the basis.
  pivot = -1 * chvatal_dict[leaving_row, entering_column]
  
  # the pivot_array expresses the newly arrived 
  # basic variable in terms of the non-basic variables
  next_pivot_row = chvatal_dict[leaving_row,:]/pivot
  
  # the coefficent of the newly arrived non-basic variable
  # it is moving from the LHS to the RHS hence the -1
  next_pivot_row[entering_column] = -1./pivot
  
  return pivot, next_pivot_row   

def get_next_dictionary(chvatal_dict, leaving_row, entering_column,
                           row_lookup, column_lookup):
  _, pivot_array = get_pivot_row(chvatal_dict,
                                 leaving_row,
                                 entering_column)

  # Update the remaining rows in the dictionary
  # so they are now expressed in terms of the new arrived
  # non-basic variable.
  # m+1 since we also need to update the row for z
  cd_candidate_rows = []
  for j in range(m+1):
    if j == leaving_row:
      cd_candidate_rows.append( pivot_array )
    else:
      # the coefficient of the leaving non-basic variable
      multiplier = chvatal_dict[j, entering_column]
      
      multiplier_times_pivot_array = multiplier * pivot_array
      updated_row_array = multiplier_times_pivot_array +  chvatal_dict[j, : ]
      
      # correct the multiplier for the newly entered non-basic variable
      updated_row_array[entering_column] = multiplier_times_pivot_array[entering_column]
      
      # print(updated_row_array)
      cd_candidate_rows.append( updated_row_array )

  cd_next_it = np.vstack( cd_candidate_rows )

  # assemble the updated list of basic variables
  basic_next_it = row_lookup.copy()
  # update the leaving row with the entering new basic variable
  basic_next_it[leaving_row] = column_lookup[entering_column]

  # assemble the updated list of nonbasic variables
  nonbasic_next_it = column_lookup.copy()
  # update the entering column with the new nonbasic variable
  nonbasic_next_it[entering_column] = row_lookup[leaving_row]

  #print(cd_next_it)

  return basic_next_it, nonbasic_next_it, cd_next_it

Helper functions

Code

def get_row_variable_names(row_lookup):
  row_names = {}
  for row_index in row_lookup:
    row_names[row_index] = variable_id_to_name[row_lookup[row_index]]
    # print(f'row_index {row_index} variable {variable_to_name[row_lookup[row_index]]}')
  # Last row is for z 
  row_names[Z_ROW] = 'z'
  return row_names

def get_column_variable_names(col_lookup):
  col_names = {}

  # The first column contains RHSs
  col_names[0] = 'RHS'
  
  for col_index in col_lookup:
    col_names[col_index] = variable_id_to_name[col_lookup[col_index]]
    # print(f'column_index {col_index} variable {variable_to_name[col_lookup[col_index]]}')
  return col_names

def get_chvatal_df(cd, row_lookup, col_lookup):
  row_names = get_row_variable_names(row_lookup)
  col_names = get_column_variable_names(col_lookup)
  return pd.DataFrame(cd,
             columns=[col_names[k] for k in sorted(col_names)],
             index = [row_names[k] for k in sorted(row_names)])

    # get the keys of nonbasic_next_iteration in the sorted order of values
    # values here are nonbasic variables. Each key is a column index 
    # ordered_columns = sorted(col_lookup, 
    #                       key=col_lookup.__getitem__)

The implementation

next_iteration = True

# iteration counter
it = 0

variable_id_to_name = dict(zip(range(n+m),[f'x_{i+1}' for i in range(n+m)]))

# variable ids that will take on zero values
nonbasic = [i for i in range(n)]

# variable ids that will take on nonzero values
basic = [i for i in range(n,n+m)]

row_lookup_it = {}
# row_lookup_it[i] is a dictionary telling us the id of the 
# decision variable corresponding to any row in the Chvátal 
# dictionary for iteration i
row_lookup_it[it] = dict( zip(range(m), basic) )

col_lookup_it = {}
# col_lookup_it[i] is a dictionary telling us the id of the 
# decision variable corresponding to any column (strictly greater than 0) 
# in the Chvátal dictionary for iteration i
# We ignore column 0 since we use it to store the right hand sides.
col_lookup_it[it] = dict( zip(range(1,n+1), nonbasic) )

z_star = 0

cd_it = {}
cd_it[it] = np.vstack( [ np.hstack([b, -A]), np.vstack( [ [z_star] , c ]).T ] )

cd_df = get_chvatal_df(cd_it[it],
                       row_lookup_it[it],
                       col_lookup_it[it])
print(f'Iteration {it}\n Chvátal dictionary \n {cd_df}')

if np.min( cd_it[it][ FIRST_ROW : ALL_BUT_LAST_ROW, FIRST_COLUMN ] ) < 0:
  # This means that some basic variable takes on a negative
  # value in the starting solution.
  print('Stopping. Infeasible starting solution.')
else:
  while next_iteration:
    entering_column = get_entering_column(cd_it[it])

    if entering_column is None:
      next_iteration = False
    else:
      # get the identifier for the variable that corresponds to
      # the entering column
      entering_var_id = col_lookup_it[it][entering_column]
      entering_var_name = variable_id_to_name[entering_var_id]
      print(f' Entering variable {entering_var_name}')

      leaving_row = get_leaving_row(cd_it[it], entering_column)

      # get the identifier for the variable that corresponds to
      # the leaving row
      leaving_var_id = row_lookup_it[it][leaving_row]
      leaving_var_name = variable_id_to_name[leaving_var_id]
      print(f' Leaving variable {leaving_var_name}\n')

      basic, nonbasic, cd = get_next_dictionary(cd_it[it],
                                                leaving_row,
                                                entering_column,
                                                row_lookup_it[it],
                                                col_lookup_it[it])
      
      # Update
      it += 1

      row_lookup_it[it] = basic
      col_lookup_it[it] = nonbasic
      cd_it[it] = cd

      # Pretty print the updated dictionary
      cd_df = get_chvatal_df(cd,
                            basic,
                            nonbasic)
      print(f'Iteration {it}\n Chvátal dictionary \n {cd_df}')

Iteration 0
 Chvátal dictionary 
      RHS  x_1  x_2  x_3
x_4    5   -2   -3   -1
x_5   11   -4   -1   -2
x_6    8   -3   -4   -2
z      0    5    4    3

 Entering variable x_1

 Leaving variable x_4

Iteration 1
 Chvátal dictionary 
       RHS  x_4  x_2  x_3
x_1   2.5 -0.5 -1.5 -0.5
x_5   1.0  2.0  5.0  0.0
x_6   0.5  1.5  0.5 -0.5
z    12.5 -2.5 -3.5  0.5

 Entering variable x_3

 Leaving variable x_6

Iteration 2
 Chvátal dictionary 
       RHS  x_4  x_2  x_6
x_1   2.0 -2.0 -2.0  1.0
x_5   1.0  2.0  5.0 -0.0
x_3   1.0  3.0  1.0 -2.0
z    13.0 -1.0 -3.0 -1.0

Next steps

We will continue to expand on this basic implementation (and rectify it’s many deficiences) in future posts.

References

[1]

V. Chvátal, Linear programming. Macmillan, 1983.