Add redundant presolve for redundant constraints

src/simplex.lisp

@@ -1,4 +1,3 @@
-
 (uiop:define-package :linear-programming/simplex
   (:use :cl
         :iterate
@@ -414,26 +413,47 @@ unchanged."
              (num-art-vars (tableau-var-count solved-art-tab))
              (num-constraints (tableau-constraint-count main-tab)))
 
-         ;; Check that all artificial variables are out of the basis
-         ;; Degeneracy can allow an artificial variable to be zero, but still in the basis
+         ;; Remove artificial variables from basis
          (iter (for basis-col in-vector (tableau-basis-columns solved-art-tab))
                (for i from 0)
            (when (>= basis-col num-vars)
-             ;; Need to find non-basis variable with a nonzero in this row
-             (when (/= 0 (aref art-matrix i num-art-vars))
+             (unless (fp= 0 (aref art-matrix i num-art-vars)
+                         (tableau-fp-tolerance-factor solved-art-tab))
                (error (format nil "Artificial variable ~S still non-zero" basis-col)))
+
              (let ((new-col -1))
+               ;; Find non-basis variable with nonzero coefficient
                (iter (for j from 0 below num-vars)
-                 (when (and (/= 0 (aref art-matrix i j))
-                            (iter (for new-col in-vector art-basis)
-                              (always (/= new-col j))))
+                 (when (and (not (fp= 0 (aref art-matrix i j)
+                                     (/ (tableau-fp-tolerance-factor solved-art-tab) 2)))
+                            (not (find j art-basis)))
                    (setf new-col j)
                    (return)))
-               (when (= new-col -1)
-                 (error "Artificial variable still in basis and cannot be replaced"))
-               (n-pivot-row solved-art-tab new-col i))))
 
-         ;; Copy tableau coefficients/RHS
+               (if (= new-col -1)
+                   ;; Handle redundant constraint case
+                   (let ((row-is-zero t))
+                     (iter (for j from 0 below num-vars)
+                       (unless (fp= 0 (aref art-matrix i j)
+                                   (/ (tableau-fp-tolerance-factor solved-art-tab) 2))
+                         (setf row-is-zero nil)
+                         (return)))
+                     (cond
+                       (row-is-zero
+                        ;; Redundant constraint - use any non-basis variable
+                        (iter (for j from 0 below num-vars)
+                          (unless (find j art-basis)
+                            (setf new-col j)
+                            (return)))
+                        (when (= new-col -1)
+                          (setf new-col 0))
+                        (when (not (fp= 0 (aref art-matrix i new-col)
+                                        (/ (tableau-fp-tolerance-factor solved-art-tab) 2)))
+                          (n-pivot-row solved-art-tab new-col i)))
+                       (t (error "Artificial variable still in basis and cannot be replaced"))))
+                   (n-pivot-row solved-art-tab new-col i)))))
+
+         ;; Copy coefficients and RHS to main tableau
          (iter (for row from 0 below num-constraints)
            (iter (for col from 0 below num-vars)
              (setf (aref main-matrix row col) (aref art-matrix row col)))
@@ -443,12 +463,13 @@ unchanged."
          ;; Update basis columns and objective row to match
          (iter (for basis-col in-vector art-basis)
                (for i from 0)
-           (setf (aref (tableau-basis-columns main-tab) i) basis-col)
-           (let ((scale (aref main-matrix num-constraints basis-col)))
-             (when (/= 0 scale)
-               (iter (for col from 0 to num-vars)
-                 (decf (aref main-matrix num-constraints col)
-                       (* scale (aref main-matrix i col))))))))
+           (when (< basis-col num-vars)
+             (setf (aref (tableau-basis-columns main-tab) i) basis-col)
+             (let ((scale (aref main-matrix num-constraints basis-col)))
+               (when (not (fp= 0 scale (/ (tableau-fp-tolerance-factor main-tab) 2)))
+                 (iter (for col from 0 to num-vars)
+                   (decf (aref main-matrix num-constraints col)
+                         (* scale (aref main-matrix i col)))))))))
        (n-solve-tableau main-tab)))
     (t
      (check-type tableau tableau)

t/simplex.lisp

@@ -403,3 +403,20 @@
              (is (= 1/2 x))
              (is (= 7 y))
              (is (= 0 z))))))
+
+(test tableau-redundant-constraints
+  (declare (notinline tableau-reduced-cost))
+  (let* ((problem (make-linear-problem (min (= n (+ b1 b2 b3 b4 b5)))
+                                       (= 7 (+ b1 b3 b4))
+                                       (= 5 (+ b4 b5))
+                                       (= 12 (+ b1 b2 b4 b5))
+                                       (= 7 (+ b1 b2 b5))
+                                       (= 2 (+ b1 b3 b5))))
+         (tableau (n-solve-tableau (build-tableau problem problem))))
+    (is (= 0 (tableau-reduced-cost tableau 'b1)))
+    (is (= 0 (tableau-reduced-cost tableau 'b2)))
+    (is (= -1 (tableau-reduced-cost tableau 'b3)))
+    (is (= 0 (tableau-reduced-cost tableau 'b4)))
+    (is (= 0 (tableau-reduced-cost tableau 'b5)))
+    (is (= 12 (tableau-variable tableau 'n)))
+    (signals error (tableau-reduced-cost tableau 'foo))))