#!/usr/bin/env python # -*- truncate-lines: t -*- # # Written by: # -- # John L. Weatherwax 2007-07-05 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- import numpy as np; from numpy import asarray, any, all, inf, vstack def full_simplex(T,verbose=False): """ Implements the flow diagram for the "full simplex method" which consists of : Figure 39 (The phase I Algorithm) followed by the flow diagram from Figure 18 (the simplex algorithm with nonnegative right-hand-sides). """ T = T.astype(float) m_plus_1, n_plus_1 = T.shape m, n = m_plus_1 - 1, n_plus_1 - 1 # Construct row/column labels for the initial tableau: # top_row, left_column, right_column, bottom_row = construct_initial_tableau_labels( m, n ) if verbose : print "Initial tableau= " print T # Set up the databox putting rows with negative b_i last: # top = T[:-1,:] last_row = T[-1,:] last_column = T[:-1,-1] # the last column (excluding the last row) neg_b_row_ordering = last_column.argsort()[::-1] # sort so negative elements in b are listed last T = vstack( ( top[ neg_b_row_ordering, : ], last_row ) ) if verbose : print "Tableau after putting negative b_i last= " print T # Reorder the right/left columns based on the ordering of b: # left_column = asarray( left_column[:-1] )[ neg_b_row_ordering ].tolist() + [ left_column[-1] ] right_column = asarray( right_column[:-1] )[ neg_b_row_ordering ].tolist() + [ right_column[-1] ] # While there is at least one negative t_{i,n+1} for i = 1, ..., m : # while any( T[:-1,-1] < 0 ) : # Select the row K with b_K < 0 and K the largest i such that b_i < 0: # K = m-1 while( T[K,n] > 0 ): K -= 1 # Are all the entries in row K positive or zero except the last? # if all( T[K,:-1] >= 0 ) : print "HALT: Maximum problem has no feasible solution; Miniumum problem has an unbounded solution" return T, (top_row, left_column, right_column, bottom_row) # Choose a column J <= (n-1) so that t_{K,J} < 0. # J = n-1 while( T[K,J] >= 0 ): J -= 1 # Choose row I >= K so that t_{I,n+1}/t_{I,J} is a minimum value of t_{i,n+1}/t_{i,j} for i=K and all i>K such that t_{i,J} > 0. # I, min_ratio = K, +inf for ii in range(I,m): if T[ii,J] > 0 : ratio = T[ii,n_plus_1-1] / T[ii,J] if ratio < min_ratio : I, min_ratio = ii, ratio if verbose : print "(K,J,I)= (%d,%d,%d); t_{I,J}= %10.6f" % (K,J,I,T[I,J]) # Perform pivoting on the element at t_{I,J}: # T, (top_row, left_column, right_column, bottom_row) = simplex_pivot( I, J, T, (top_row, left_column, right_column, bottom_row), verbose=verbose ) #endwhile # Tableau now satisfies nonnegativity assumptions and can be solved using Figure 18: # T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T, (top_row, left_column, right_column, bottom_row), verbose=verbose ) return T, (top_row, left_column, right_column, bottom_row) def one_or_more_positive_indicators(T): """ Returns true if we have at least one positive indicator. """ indicators = T[-1,:-1] # the last row and all but the last column return any( indicators > 0 ) def construct_initial_tableau_labels(m,n): """ """ top_row = [ ( 'x_%d' % (_+1) ) for _ in range(n) ] + [ '-1' ] left_column = [ ( 'v_%d' % (_+1) ) for _ in range(m) ] + [ '-1' ] right_column = [ ( '= -y_%d' % (_+1) ) for _ in range(m) ] + [ '= f' ] bottom_row = [ ( '= u_%d' % (_+1) ) for _ in range(n) ] + [ '= g' ] return top_row, left_column, right_column, bottom_row def exchange_row_column_labels( I, J, top_row, left_column, right_column, bottom_row ): """ """ tr_j, rc_i = top_row[J], right_column[I] top_row[J] = rc_i.replace('= -','') right_column[I] = '= -' + tr_j br_j, lc_i = bottom_row[J], left_column[I] bottom_row[J] = '= ' + lc_i left_column[I] = br_j.replace('= ','') return top_row, left_column, right_column, bottom_row def simplex_pivot( I, J, T, (top_row, left_column, right_column, bottom_row), verbose=False ): """ """ m_plus_1, n_plus_1 = T.shape # Divide the Ith row of the old tableau by the pivot T_{I,J} and enter the results in the new tableau # T_new = T.copy() T_new[I,:] = T_new[I,:] / T_new[I,J] # For each row i (except the Ith) subtract T_{i,J} times the Ith row of the new tableau from # the ith row of the old tableau and enter the results in the new tableau: # for ii in range(m_plus_1) : if ii==I : continue T_new[ii,:] = T[ii,:] - T[ii,J] * T_new[I,:] if verbose : print "after column elimination (flow diagram box 8)= " print T_new # 1) Replace the I,J entry in the new tableau by 1/T_{I,J}. # 2) Replace all other entries in the pivotal column by -T_{i,J} / T_{I,J} # T_new[I,J] = 1./T[I,J] for ii in range(m_plus_1) : if ii==I : continue T_new[ii,J] = -T[ii,J] / T[I,J] if verbose : print "after pivot column replacement (flow diagram box 9)= " print T_new # Label the new tableau the same as the old tableau except: # # 1) interchange the variable at the top of column J with the variable at the right of row I # 2) interchange the variable at the bottom of column J with the variable at the left of row I # top_row, left_column, right_column, bottom_row = exchange_row_column_labels( I, J, top_row, left_column, right_column, bottom_row ) # Update tableau: T = T_new.copy() return T, (top_row, left_column, right_column, bottom_row) def print_tableau( T, (top_row, left_column, right_column, bottom_row) ): """ """ m_plus_1, n_plus_1 = T.shape m, n = m_plus_1 - 1, n_plus_1 - 1 print "Final tableau= " print "%12s" % "" , for jj in range(n_plus_1) : print "%12s " % top_row[jj] , print "" for ii in range(m_plus_1) : print "%12s " % left_column[ii] , for jj in range(n_plus_1) : print "%12.4f " % T[ii,jj] , print "%-12s " % right_column[ii] print "%12s" % "" , for jj in range(n_plus_1) : print "%12s " % bottom_row[jj] , print "" def nonnegative_rhs_simplex(T,tableau_labels=None,verbose=False): """ Implements the flow diagram from Figure 18 (The simplex algorithm for problems with nonnegative right-hand-sides) Input: T the initial tableau """ assert all( T[:-1,-1] >= 0 ), "this version of the simplex method expects nonnegative right-hand-sides b" T = T.astype(float) m_plus_1, n_plus_1 = T.shape m, n = m_plus_1 - 1, n_plus_1 - 1 # Construct row/column labels for the tableau: # if tableau_labels is not None : top_row, left_column, right_column, bottom_row = tableau_labels else : top_row, left_column, right_column, bottom_row = construct_initial_tableau_labels( m, n ) while one_or_more_positive_indicators(T) : # Select a positive indicator from column J: # we follow the heuristic of selecting the column with most positive indicator # indicators = T[-1,:-1] J = indicators.argmax() Jth_column = T[:,J] if all( Jth_column[:-1] <= 0 ) : print "HALT: Maximum problem has unbounded solution; Miniumum problem has no feasible solution" return T, (top_row, left_column, right_column, bottom_row) # Choose row I so that the ratio ( T_{I,n+1} / T_{I,J} ) is the smallest (over all elements where T_{I,J} > 0) # I, min_ratio = -1, +inf for ii in range(m): if T[ii,J] > 0 : ratio = T[ii,n_plus_1-1] / T[ii,J] if ratio < min_ratio : I, min_ratio = ii, ratio if verbose : print "(J,I)= (%d,%d)" % (J,I) # Perform pivoting on the element at t_{I,J}: # T, (top_row, left_column, right_column, bottom_row) = simplex_pivot( I, J, T, (top_row, left_column, right_column, bottom_row) ) #endwhile # Computation ended. Basic variables have the indicated values and the nonbasis variables are zero. # if verbose : print_tableau( T, (top_row, left_column, right_column, bottom_row) ) return T, (top_row, left_column, right_column, bottom_row) if __name__ == "__main__" : # For Section 5: # print "Duplicate Example 1 from Section 5: " T = asarray( [ [ 5, 2, 180 ], [ 3, 3, 135 ], [300, 200, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Duplicate Example 2 from Section 5: " T = asarray( [ [ 6, 2, 4, 200 ], [ 2, 2, 12, 160 ], [ 12, 8, 24, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Example 3 from Section 3 (Tableau modified from Figure 10 in the book): " T = asarray( [ [ 6, 2, 4, 160 ], [ 2, 2, 12, 200 ], [ 12, 8, 24, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Example 4 from Section 3 (Tableau modified from Figure 10 in the book): Gives one of the infinitely many solutions" T = asarray( [ [ 6, 2, 4, 200 ], [ 2, 2, 12, 200 ], [ 12, 8, 24, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 3 from Section 3:" T = asarray( [ [ 3, 6, 2100 ], [ 6, 5, 2100 ], [ 40, 50, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 4 from Section 3:" T = asarray( [ [ 0.15, 1.67, 8 ], [ 0.1, 3.33, 10 ], [ 0.12, 2.0, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 5 from Section 3:" T = asarray( [ [ 1, 1, 200 ], [ 1, 4, 320 ], [ 10, 20, 2200 ], [ 40, 60, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 6 from Section 3:" T = asarray( [ [ 1, 1, 200 ], [ 1, 4, 320 ], [ 10, 20, 2200 ], [ 40, 90, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 7 from Section 3:" T = asarray( [ [ 1, 3, 12, 1 ], [ 3, 4, 1, 0.25 ], [ 8, 19, 7, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 8 from Section 3:" T = asarray( [ [ 1, 3, 12, 1, 1 ], [ 3, 4, 1, 0, 0.25 ], [ 8, 19, 7, 0.5, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 10 from Section 3:" T = asarray( [ [ 0.01, 0.02, 400 ], [ 0.03, 0.01, 450 ], [ 0.03, 0.015, 480 ], [ 10, 15, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 11 from Section 5:" T = asarray( [ [ 0.3, 0.8, 121 ], [ 0.7, 0.2, 49 ], [ 80, 50, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 12 from Section 5:" T = asarray( [ [ 0.5, 0.25, 0.2, 6 ], [ 0.4, 0.3, 0.25, 5 ], [ 50000, 30000, 10000, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 12 from Section 6:" T = asarray( [ [ 0.5, 0.25, 0.2, 6 ], [ 0.4, 0.3, 0.25, 5 ], [ 36000, 12000, 18000, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = nonnegative_rhs_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) # For Section 6: # print "Exercise 3 (b) from Section 6: " T = asarray( [ [ 5, 2, 180 ], [ 3, 3, 144 ], [ 300, 200, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 4 from Section 6: " T = asarray( [ [ 6, 2, 4, 200 ], [ 2, 2, 12, 160 ], [ 16, 8, 24, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 5 (b) from Section 6: " T = asarray( [ [ 6, 2, 4, 200 ], [ 2, 2, 12, 160 ], [ 16, 10, 24, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 6 from Section 6: " T = asarray( [ [ 6, 2, 4, 200 ], [ 2, 2, 12, 160 ], [ 12, 8, 56, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 7 from Section 6: " T = asarray( [ [ 5, 2, 180 ], [ 3, 3, 153 ], [ 300, 200, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 27 from Section 6: " T = asarray( [ [ -1, 1, 4 ], [ 2, -4, 8 ], [ 2, 3, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) # For Section 7: # print "Duplicate Example 2 from Section 7: " T = asarray( [ [ 1, 1, 20 ], [ 1, 2, 30 ], [ -1, -2, -30 ], [ 2, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 7 (a) from Section 7: " T = asarray( [ [ 1, 1, 10 ], [ -1, -1, -6 ], [ 1, 0, 8 ], [ 2, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 7 (b) from Section 7: " T = asarray( [ [ -1, 0, -2 ], [ 0, -1, -3 ], [ 3, 2, 24 ], [ 1, 0, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 8 from Section 7: " T = asarray( [ [ 1, 1, 20 ], [ -1, -2, -50 ], [ 2, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) # For Section 11: # print "Exercise 2 (a) from Section 11: " T = asarray( [ [ 8, 7, 10, 1 ], [ 5, 10, 7, 1 ], [ 3, 12, 1, 1 ], [ 1, 1, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 2 (b) from Section 11: " T = asarray( [ [ 5, 10, 7, 1, 1 ], [ 10, 3, 12, 14, 1 ], [ 3, 12, 1, 7, 1 ], [ 1, 1, 1, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 3 from Section 11: " T = asarray( [ [ 1, -1, 1 ], [ -1, 1, 1 ], [ 1, 1, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) # For Section 12: # print "The Tasty Nut Company Example from Section 12: " T = asarray( [ [ 0.8, 0.5, 0.3, 82 ], [ 0.15, 0.2, 0.2, 30 ], [ 0.05, 0.2, 0.2, 30 ], [ 0.0, 0.05, 0.15, 16 ], [ 0.0, 0.05, 0.1, 10 ], [ 0.0, 0.0, 0.05, 4 ], [ 0.25, 0.4, 0.75, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "The Tasty Nut Company Example from Section 12 (with an extra 1/2 lb of Brazil nuts): " T = asarray( [ [ 0.8, 0.5, 0.3, 82 ], [ 0.15, 0.2, 0.2, 30 ], [ 0.05, 0.2, 0.2, 30 ], [ 0.0, 0.05, 0.15, 16 ], [ 0.0, 0.05, 0.1, 10 ], [ 0.0, 0.0, 0.05, 4.5 ], [ 0.25, 0.4, 0.75, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) ) print "Exercise 7 from Section 12: " T = asarray( [ [ 0.8, 0.5, 0.3, 82 ], [ 0.15, 0.2, 0.2, 30 ], [ 0.05, 0.2, 0.2, 30 ], [ 0.0, 0.05, 0.15, 16 ], [ 0.0, 0.05, 0.1, 10 ], [ 0.0, 0.0, 0.05, 4.5 ], [ 0.30, 0.50, 0.90, 0 ] ] ) T, (top_row, left_column, right_column, bottom_row) = full_simplex( T ) print_tableau( T, (top_row, left_column, right_column, bottom_row) )