working

Author: scheinerman

development/examples.jl

@@ -183,6 +183,17 @@ function phase_one_example()
     return Tableau(A, b, c, false)
 end
 
+function phase_one_infeasible()
+    A = [
+        2 0 5 2
+        2 1 3 1
+        -1 1 4 3
+    ]
+    b = [5, -1, 2]
+    c = [1, -1, 4, 4]
+    return Tableau(A, b, c, false)
+end
+
 function phase_one_trouble()
     A = [5 1 3; 1 1 3]
     b = [5, 5]

src/SimplexTableaux.jl

@@ -20,6 +20,7 @@ include("Exact.jl")
 export Tableau,
     basic_vector,
     big_M_solve,
+    big_M_tableau,
     check_basis,
     set_basis!,
     find_a_basis,

src/big_M.jl

@@ -17,7 +17,7 @@ function big_M_tableau(T::Tableau, M::Int=1000)
     AA = hcat(AA, eye(Int, m))
 
     # create artifical c 
-    cc = vcat(0*c, M*ones(Int, m))
+    cc = vcat(c, M*ones(Int, m))
 
     TT = Tableau(AA, bb, cc, false)
 
@@ -34,11 +34,25 @@ Solve the LP `T` using the big-M method.
 """
 function big_M_solve(T::Tableau, M::Int=1000)
     TT = big_M_tableau(T, M)
+    println("Solving augmented tableau")
     x = simplex_solve!(TT)
+
+    if isnothing(x)
+        return nothing
+    end
+    # check that artificial vars are all zero 
+
+    arts = x[(T.n_vars + 1):end]
+
+    if any(arts .≠ 0)
+        @info "LP is possibly infeasible. Try a larger value than M = $M?"
+        return nothing
+    end
+
     x = x[1:T.n_vars]
     B = infer_basis!(T, x)
     set_basis!(T, B)
-    println("Final tableau\n")
+    println("\nFinal tableau\n")
     println(T)
     v = value(T)
     println("Minimial value = $v = $(Float64(v))")