## Linear programming

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.

lpsolve-tcl : Tcl extension for LPSolve , a Mixed Integer Linear Programming (MILP) solver.

 Category Mathematics Category Numerical Analysis