min f(x) Ax = b x >= 0where 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
}
trace "problem ready for startup"
## 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"
}See also:http://en.wikipedia.org/wiki/Simplex_algorithm
Code to be replaced using math, part of tcllib
