## Operations research - Simplex Algorithm

FF - 2007-06-10 - Simplex Algorithm, created by the North American mathematician George Dantzig in 1947, is a popular technique for numerical solution of the linear programming problem. This implementation deals with problems expressed in standard form, that means: minimization of a linear function respecting a number of constraints (linear), and that all variables are positive.

We could write the problem (in standard form) in a matrix way:

``` min f(x)
Ax = b
x >= 0```

where f(x) is a linear function R^n -> R and Ax = b is a linear system of m equations (in matrix form) in n variables, representing the constraints.

The code below implements Phase 2 and mostly Phase 1 of the algorithm. What I miss is keep going from Phase 1 to Phase 2, when Phase 1 result is not optimal.

``` namespace import ::math::linearalgebra::*

proc trace {arg} {if {0} {puts "\$arg"}}
proc error {arg} {if {1} {puts "\$arg"}}

proc simplex_phase2 {l} {
## prepare matrix (fix float numbers)
lassign \$l cVect cMatrix bVect
lassign [shape \$cMatrix] rows cols
set nvars [llength \$cVect]
set nconstr [llength \$bVect]
set l [list [scale_vect 1.0 \$cVect] [scale_mat 1.0 \$cMatrix] [scale_vect 1.0 \$bVect]]
trace "cVect:\n[show \$cVect]\ncMatrix:\n[show \$cMatrix]\nbVect:\n[show \$bVect]"
if {![expr {\$nvars==\$cols&&\$nconstr==\$rows}]} {
error "problem is malformed (unmatching rows/cols count)."
exit
}

set fix_b_sign 0
for {set row 0} {\$row < \$rows} {incr row} {
set bVect_i [lindex \$bVect \$row]
if {\$bVect_i < 0} {
set cMatrix_i [getrow \$cMatrix \$row]
lset cMatrix \$row [scale -1.0 \$cMatrix_i]
lset bVect \$row [expr -1.0*\$bVect_i]
set fix_b_sign 1
}
}

if \$fix_b_sign {trace "fixed cMatrix/bVect:\n[show \$cMatrix]\n[show \$bVect]"}

## find identity matrix (if any), that is:
## choose indices for in-base vars, and out-base vars
lassign [shape \$cMatrix] nc n
set I [mkIdentity \$nc]
set iB {}
set iN {}
set cMt [transpose \$cMatrix]
foreach Ir \$I {
set idx [lsearch -exact \$cMt \$Ir]
if {\$idx != -1} {lappend iB \$idx}
}
for {set iout 0} {\$iout < \$n} {incr iout} {
if {[lsearch -exact \$iB \$iout] == -1} {
lappend iN \$iout
}
}
unset cMt I Ir

trace "iB=\$iB, iN=\$iN"

if {[llength \$iB] < \$rows} {
trace "### not ready for startup; entering phase 1"
set new_cVect {}
set n_artificial_vars 0
set new_cMatrix {}
set cMt [transpose \$cMatrix]
set cnt 0
foreach cMatrix_i \$cMatrix {
set e_k {}
for {set j 0} {\$j < [llength \$cMatrix]} {incr j} {
lappend e_k [expr \$j==\$cnt?1.0:0.0]
}
incr cnt
if {[lsearch -exact \$cMt \$e_k] == -1} {
incr n_artificial_vars
lappend cMt \$e_k
}
}
set new_cMatrix [transpose \$cMt]
for {set j 0} {\$j < \$cols} {incr j} {
lappend new_cVect 0.0
}
for {set j 0} {\$j < \$n_artificial_vars} {incr j} {
lappend new_cVect 1.0
}
set result [simplex_phase2 [list \$new_cVect \$new_cMatrix \$bVect]]
array set _ [set result]
trace "### end phase1: \$result"
for {set j \$cols} {\$j < [expr \$cols+\$n_artificial_vars]} {incr j} {
if [expr [lindex \$_(result:) \$j]!=0.0?1:0] {
#TODO: continue swapping artificial vars with original vars
trace "### unfeasible problem - artificial vars remain in base"
return -1
}
}
return \$result
}

## algorithm startup
set Bt {}
set cb {}
foreach ib \$iB {
lappend Bt [getcol \$cMatrix \$ib]
lappend cb [lindex \$cVect \$ib]
}
set B [transpose \$Bt]
# inversion not required for first step, cause we have an identity matrix
#set Binv [invert \$B]
set Binv \$B
set Nt {}
set cn {}
foreach in \$iN {
lappend Nt [getcol \$cMatrix \$in]
lappend cn [lindex \$cVect \$in]
}
set N [transpose \$Nt]
set BinvN [matmul \$Binv \$N]
set Binvb [matmul \$Binv \$bVect]

set iteration 0
set iB_history {}
### LOOP ENTER POINT
while {1} {
trace "*********** Iteration #[incr iteration] ******************"

## compute gamma (reduced costs) vector
set gamma [axpy -1 [matmul \$cb \$BinvN] \$cn]
trace "gamma:\n[show \$gamma]"

set h 0
foreach gamma_i \$gamma {if {\$gamma_i < 0} {break} {incr h}}
## check for optimal solution
if {\$h >= [llength \$gamma]} {
## we found an optimal solution
## check uniqueness:
set unique 1
foreach gamma_i \$gamma {if {\$gamma_i == 0} {set unique 0}}

## now let's build our x vector (x_star)
set x_star {}
for {set ix 0} {\$ix < \$cols} {incr ix} {
set ib_idx [lsearch -exact \$iB \$ix]
if {\$ib_idx != -1} {
lappend x_star [lindex \$Binvb \$ib_idx]
} else {
lappend x_star 0.0
}
}
#warn "found optimal solution: \$x_star"
#warn "uniqueness: \$unique"
#warn "value at optimum: [dotproduct \$x_star \$cVect]"
return [list result: \$x_star iB: \$iB gamma: \$gamma]
} else {
## or unbound problem:
for {set hi 0} {\$hi < [llength \$gamma]} {incr hi} {
## foreach negative component of gamma, check if the B^(-1)N matching
## column is all < 0
if {[lindex \$gamma \$hi] < 0} {
set BinvN_h [getcol \$BinvN \$hi]
set BinvN_neg_cnt 0
foreach BinvN_h_i \$BinvN_h {if {\$BinvN_h_i < 0} {incr BinvN_neg_cnt}}

if {\$BinvN_neg_cnt == [llength \$BinvN_h]} {
#warn "found unbound problem"
return -1
}
}
}
}

## - ingoing var and outgoing var -
## find h and k (using Bland's anti-loop rule)
set hmin [lindex \$gamma 0]
set h 0
for {set j 1} {\$j < [llength \$gamma]} {incr j} {
set newhmin [lindex \$gamma \$j]
if {\$newhmin < [lindex \$hmin 0]} {
set h \$j
set hmin \$newhmin
} elseif {\$newhmin < [lindex \$hmin 0]} {
lappend h \$j
lappend hmin \$newhmin
}
}
set h [lindex \$h 0]
set hmin [lindex \$hmin 0]

set pi_h [getcol \$BinvN \$h]

set kmin [expr [lindex \$Binvb 0]/[lindex \$pi_h 0]]
set k 0
for {set j 1} {\$j < [llength \$pi_h]} {incr j} {
set newkmin [expr [lindex \$Binvb \$j]/[lindex \$pi_h \$j]]
if {(\$newkmin < \$kmin && \$newkmin >= 0) || \$kmin < 0} {
set k \$j
set kmin \$newkmin
}
}

trace "h=\$h, k=\$k, rho=\$kmin"

## compute new iB,iN (indices for in-base/out-base vars)
## swapping h-th in-base var with k-th out-base var
set iB_old \$iB
set iN_old \$iN
lset iB \$k [lindex \$iN_old \$h]
lset iN \$h [lindex \$iB_old \$k]

trace "iB=\$iB, iN=\$iN"

## pivot operation - prepare M matrix

set M {}
lappend M [getcol \$BinvN \$h]
lassign [shape \$BinvN] BinvN_rows BinvN_cols
for {set ij 0} {\$ij < \$BinvN_cols} {incr ij} {
if {\$ij == \$h} {
## e_k here
set e_k [getcol \$BinvN \$h]
set e_k_i 0
foreach e_k___ \$e_k {
lset e_k \$e_k_i [expr (\$e_k_i==\$k)?1.0:0.0]
incr e_k_i
}
lappend M \$e_k
} else {
set BinvN_col [getcol \$BinvN \$ij]
lappend M \$BinvN_col
}
}
lappend M \$Binvb
set M [transpose \$M]

trace "before pivot:\nM:\n[show \$M]"

## pivot operation:
set row_k [getrow \$M \$k]
lset M \$k [scale [expr 1/[lindex [getcol \$M 0] \$k]] \$row_k]
set row_k_neg [scale -1.0 [getrow \$M \$k]]
for {set ij 0} {\$ij < [llength \$M]} {incr ij} {
if {\$ij == \$k} {continue}
set cur_row [getrow \$M \$ij]
lset M \$ij [add_vect \$cur_row [scale [lindex \$cur_row 0] \$row_k_neg]]
}

trace "after pivot:\nM:\n[show \$M]"

## new BinvN, Binvb from M-pivoted
set BinvN {}
for {set ij 1} {\$ij <= \$BinvN_cols} {incr ij} {
lappend BinvN [getcol \$M \$ij]
}
set BinvN [transpose \$BinvN]
set Binvb [getcol \$M end]

set cb {}
set cn {}
set Bt {}
set Nt {}
foreach ib \$iB {
lappend Bt [getcol \$cMatrix \$ib]
lappend cb [lindex \$cVect \$ib]
}
foreach in \$iN {
lappend Nt [getcol \$cMatrix \$in]
lappend cn [lindex \$cVect \$in]
}
set B [transpose \$Bt]
set N [transpose \$Nt]

trace "new BinvN:\n[show \$BinvN]"
trace "new Binvb:\n[show \$Binvb]"
}
}

set f {-4 -1 -1}
set c {
{2 1 2}
{3 3 1}
}
set b {4 3}
set l [list \$f \$c \$b]

set result [simplex_phase2 \$l]

if {"\$result" == "-1"} {
puts "unfeasible or unbound problem"
} else {
puts "\$result"
}```