[Arjen Markus] (9 june 2004) ''Linear programming'' is the name for a class of optimisation problems. Such problems may arise in operations research, where you have limited resources and other constraints and want to find the parameters for which a gain function is maximum.

Such problems are formulated like:

  Find the values of x and y for which f(x,y) = x + 2y is maximal,
  given that:
     x + y <= 10
  and
    2x - y <= 3

The script below implements the classic ''simplex method'' for solving such problems.
Note that it is not as robust as should be - it does not use a two-phase solution
method, necessary for real-world problems - at least, not yet.

I intend to submit this to Tcllib.

For some more background information, see for instance 
http://www-unix.mcs.anl.gov/otc/Guide/faq/linear-programming-faq.html

[AM] (10 june 2004) The code below is flawed, as it does not work if one or more of the constraints are superfluous for pinning down the solution. I will have to look into this ...

[AM] (25 august 2005) The improved code has been committed to Tcllib. Please use that code.
----

======
 # linprog.tcl --
 #    Script to solve linear programs using the Simplex method
 #

 # math::optimize --
 #    Namespace for the commands
 #
 namespace eval ::math::optimize {
    namespace export solveLinearProgram \
                     linearProgramMaximum
 }

 # solveLinearProgram
 #    Solve a linear program in standard form
 #
 # Arguments:
 #    objective     Vector defining the objective function
 #    constraints   Matrix of constraints (as a list of lists)
 #
 # Return value:
 #    Computed values for the coordinates or "unbounded" or "infeasible"
 #
 proc ::math::optimize::solveLinearProgram { objective constraints } {
    #
    # Check the arguments first and then put them in a more convenient
    # form
    #

    foreach {nconst nvars matrix} \
       [SimplexPrepareMatrix $objective $constraints] {break}

    #SimplexPrintMatrix $matrix

    set solution [SimplexSolve $nconst nvars $matrix]

    if { [llength $solution] > 1 } {
       return [lrange $solution 0 [expr {$nvars-1}]]
    } else {
       return $solution
    }
 }

 # linearProgramMaximum --
 #    Compute the value attained at the optimum
 #
 # Arguments:
 #    objective     The coefficients of the objective function
 #    result        The coordinate values as obtained by solving the program
 #
 # Return value:
 #    Value at the maximum point
 #
 proc ::math::optimize::linearProgramMaximum {objective result} {

    set value    0.0

    foreach coeff $objective coord $result {
       set value [expr {$value+$coeff*$coord}]
    }

    return $value
 }

 # SimplexPrintMatrix
 #    Debugging routine: print the matrix in easy to read form
 #
 # Arguments:
 #    matrix        Matrix to be printed
 #
 # Return value:
 #    None
 #
 proc ::math::optimize::SimplexPrintMatrix {matrix} {
    puts "\nBasis:\t[join [lindex $matrix 0] \t]"
    foreach col [lrange $matrix 1 end] {
       puts "      \t[join $col \t]"
    }
 }

 # SimplexPrepareMatrix
 #    Prepare the standard tableau from all program data
 #
 # Arguments:
 #    objective     Vector defining the objective function
 #    constraints   Matrix of constraints (as a list of lists)
 #
 # Return value:
 #    List of values as a standard tableau and two values
 #    for the sizes
 #
 proc ::math::optimize::SimplexPrepareMatrix {objective constraints} {

    #
    # Check the arguments first
    #
    set nconst [llength $constraints]
    set ncols {}
    foreach row $constraints {
       if { $ncols == {} } {
          set ncols [llength $row]
       } else {
          if { $ncols != [llength $row] } {
             return -code error -errorcode ARGS "Incorrectly formed constraints matrix"
          }
       }
    }

    set nvars [expr {$ncols-1}]

    if { [llength $objective] != $nvars } {
       return -code error -errorcode ARGS "Incorrect length for objective vector"
    }

    #
    # Set up the tableau:
    # Easiest manipulations if we store the columns first
    # So:
    # - First column is the list of variable indices in the basis
    # - Second column is the list of maximum values
    # - "nvars" columns that follow: the coefficients for the actual
    #   variables
    # - last "nconst" columns: the slack variables
    #
    set matrix   [list]
    set lastrow  [concat $objective [list 0.0]]

    set newcol   [list]
    for {set idx 0} {$idx < $nconst} {incr idx} {
       lappend newcol [expr {$nvars+$idx}]
    }
    lappend newcol "?"
    lappend matrix $newcol

    set zvector [list]
    foreach row $constraints {
       lappend zvector [lindex $row end]
    }
    lappend zvector 0.0
    lappend matrix $zvector

    for {set idx 0} {$idx < $nvars} {incr idx} {
       set newcol [list]
       foreach row $constraints {
          lappend newcol [expr {double([lindex $row $idx])}]
       }
       lappend newcol [expr {-double([lindex $lastrow $idx])}]
       lappend matrix $newcol
    }

    return [list $nconst $nvars $matrix]
 }

 # SimplexSolve --
 #    Solve the given linear program using the simplex method
 #
 # Arguments:
 #    nconst        Number of constraints
 #    nvars         Number of actual variables
 #    tableau       Standard tableau (as a list of columns)
 #
 # Return value:
 #    List of values for the actual variables
 #
 proc ::math::optimize::SimplexSolve {nconst nvars tableau} {
    set end 0
    while { !$end } {

       #
       # Find the new variable to put in the basis
       #
       set nextcol [SimplexFindNextColumn $tableau]
       if { $nextcol == -1 } {
          set end 1
          continue
       }

       #
       # Now determine which one should leave
       # TODO: is a lack of a proper row indeed an
       #       indication of the infeasibility?
       #
       set nextrow [SimplexFindNextRow $tableau $nextcol]
       if { $nextrow == -1 } {
          return "infeasible"
       }

       #
       # Make the vector for sweeping through the tableau
       #
       set vector [SimplexMakeVector $tableau $nextcol $nextrow]

       #
       # Sweep through the tableau
       #
       set tableau [SimplexNewTableau $tableau $nextcol $nextrow $vector]
       #SimplexPrintMatrix $tableau
    }

    #
    # Now we can return the result
    #
    SimplexResult $tableau
 }

 # SimplexResult --
 #    Reconstruct the result vector
 #
 # Arguments:
 #    tableau       Standard tableau (as a list of columns)
 #
 # Return value:
 #    Vector of values representing the maximum point
 #
 proc ::math::optimize::SimplexResult {tableau} {
    set result {}

    set firstcol  [lindex $tableau 0]
    set secondcol [lindex $tableau 1]
    set result    {}

    set nvars     [expr {[llength $tableau]-2}]
    for {set i 0} {$i < $nvars } { incr i } {
       lappend result 0.0
    }

    set idx 0
    foreach col [lrange $firstcol 0 end-1] {
       set result [lreplace $result $col $col [lindex $secondcol $idx]]
       incr idx
    }

    return $result
 }

 # SimplexFindNextColumn --
 #    Find the next column - the one with the largest negative
 #    coefficient
 #
 # Arguments:
 #    tableau       Standard tableau (as a list of columns)
 #
 # Return value:
 #    Index of the column
 #
 proc ::math::optimize::SimplexFindNextColumn {tableau} {
    set idx        0
    set minidx    -1
    set mincoeff   0.0

    foreach col [lrange $tableau 2 end] {
       set coeff [lindex $col end]
       if { $coeff < 0.0 } {
          if { $coeff < $mincoeff } {
             set minidx $idx
             set mincoeff $coeff
          }
       }
       incr idx
    }

    return $minidx
 }

 # SimplexFindNextRow --
 #    Find the next row - the one with the largest negative
 #    coefficient
 #
 # Arguments:
 #    tableau       Standard tableau (as a list of columns)
 #    nextcol       Index of the variable that will replace this one
 #
 # Return value:
 #    Index of the row
 #
 proc ::math::optimize::SimplexFindNextRow {tableau nextcol} {
    set idx        0
    set minidx    -1
    set mincoeff   {}

    set bvalues [lrange [lindex $tableau 1] 0 end-1]
    set yvalues [lrange [lindex $tableau [expr {2+$nextcol}]] 0 end-1]

    foreach rowcoeff $bvalues divcoeff $yvalues {
       if { $divcoeff > 0.0 } {
          set coeff [expr {$rowcoeff/$divcoeff}]

          if { $mincoeff == {} || $coeff < $mincoeff } {
             set minidx $idx
             set mincoeff $coeff
          }
       }
       incr idx
    }

    return $minidx
 }

 # SimplexMakeVector --
 #    Make the "sweep" vector
 #
 # Arguments:
 #    tableau       Standard tableau (as a list of columns)
 #    nextcol       Index of the variable that will replace this one
 #    nextrow       Index of the variable in the base that will be replaced
 #
 # Return value:
 #    Vector to be used to update the coefficients of the tableau
 #
 proc ::math::optimize::SimplexMakeVector {tableau nextcol nextrow} {

    set idx      0
    set vector   {}
    set column   [lindex $tableau [expr {2+$nextcol}]]
    set divcoeff [lindex $column $nextrow]

    foreach colcoeff $column {
       if { $idx != $nextrow } {
          set coeff [expr {-$colcoeff/$divcoeff}]
       } else {
          set coeff [expr {1.0/$divcoeff-1.0}]
       }
       lappend vector $coeff
       incr idx
    }

    return $vector
 }

 # SimplexNewTableau --
 #    Sweep through the tableau and create the new one
 #
 # Arguments:
 #    tableau       Standard tableau (as a list of columns)
 #    nextcol       Index of the variable that will replace this one
 #    nextrow       Index of the variable in the base that will be replaced
 #    vector        Vector to sweep with
 #
 # Return value:
 #    New tableau
 #
 proc ::math::optimize::SimplexNewTableau {tableau nextcol nextrow vector} {

    #
    # The first column: replace the nextrow-th element
    # The second column: replace the value at the nextrow-th element
    # For all the others: the same receipe
    #
    set firstcol   [lreplace [lindex $tableau 0] $nextrow $nextrow $nextcol]
    set newtableau [list $firstcol]

    #
    # The rest of the matrix
    #
    foreach column [lrange $tableau 1 end] {
       set yval   [lindex $column $nextrow]
       set newcol {}
       foreach c $column vcoeff $vector {
          set newval [expr {$c+$yval*$vcoeff}]
          lappend newcol $newval
       }
       lappend newtableau $newcol
    }

    return $newtableau
 }

 #
 # Some simple tests
 #
 namespace import ::math::optimize::*

 puts "Linear program:"
 puts "Maximum of 1.7x+1.8y with:"
 puts "   1.1x + 2.2y <= 1.3"
 puts "   2.4x + 1.5y <= 1.6"
 set result [solveLinearProgram \
     { 1.7   1.8  } \
    {{ 1.1   2.2  1.3 }
     { 2.4   1.5  1.6 }}]
 puts "Solution: $result"
 puts "Maximum attained: [linearProgramMaximum {1.7 1.8} $result]"

 puts "Linear program:"
 puts "Maximum of 1.7x+1.8y with:"
 puts "      x - 2y   <= 5"
 puts "     2x - y    <= 5"
 set result [solveLinearProgram \
     { 1.7   1.8  } \
    {{ 1.0   2.0  5.0 }
     { 2.0   1.0  5.0 }}]
 puts "Solution: $result"
 puts "Maximum attained: [linearProgramMaximum {1.7 1.8} $result]"

 puts "Linear program:"
 puts "Maximum of 1.7x+1.8y with:"
 puts "     2x - y    <= 5"
 set result [solveLinearProgram \
     { 1.7   1.8  } \
    {{ 2.0   1.0  5.0 }}]
 puts "Solution: $result"
 puts "Maximum attained: [linearProgramMaximum {1.7 1.8} $result]"
======

----
[Derek Peschel] May 22, 2005 -- George B. Dantzig, the inventor of the simplex method, died on May 13 at the age of 90.

<<categories>> Mathematics | Numerical Analysis