FIX: linprog_simplex: Fix bug in Phase 1

Fix #725

Author: oyamad

quantecon/optimize/linprog_simplex.py

@@ -452,13 +452,31 @@ def solve_phase_1(tableau, basis, max_iter=10**6, piv_options=PivOptions()):
     See `solve_tableau`.
 
     """
+    L = tableau.shape[0] - 1
+    nm = tableau.shape[1] - (L+1)  # n + m
+
     success, status, num_iter_1 = \
         solve_tableau(tableau, basis, max_iter, skip_aux=False,
                       piv_options=piv_options)
 
-    if success and tableau[-1, -1] > piv_options.fea_tol:  # Infeasible
+    if not success:  # max_iter exceeded
+        return success, status, num_iter_1
+    if tableau[-1, -1] > piv_options.fea_tol:  # Infeasible
         success = False
         status = 2
+        return success, status, num_iter_1
+
+    # Check artifial variables have been eliminated
+    tol_piv = piv_options.tol_piv
+    for i in range(L):
+        if basis[i] >= nm:  # Artifial variable not eliminated
+            for j in range(nm):
+                if tableau[i, j] < -tol_piv or \
+                   tableau[i, j] > tol_piv:  # Treated nonzero
+                    _pivoting(tableau, j, i)
+                    basis[i] = j
+                    num_iter_1 += 1
+                    break
 
     return success, status, num_iter_1
 

quantecon/optimize/tests/test_linprog_simplex.py

@@ -239,3 +239,15 @@ def test_bug_8662(self):
         desired_fun = -36.0000000000
         res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
         _assert_success(res, c, b_ub=b_ub, desired_fun=desired_fun)
+
+
+class TestLinprogSimplex:
+    def test_phase_1_bug_725(self):
+        # Identified a bug in Phase 1
+        # https://github.com/QuantEcon/QuantEcon.py/issues/725
+        c = np.array([-4.09555556, 4.59044444])
+        A_ub = np.array([[1, 0.1], [-1, -0.1], [1, 1]])
+        b_ub = np.array([9.1, -0.1, 0.1])
+        desired_x = [0.1, 0]
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_success(res, c, b_ub=b_ub, desired_x=desired_x)