diff --git a/src/simplex.lisp b/src/simplex.lisp index 81d038c..13d076e 100644 --- a/src/simplex.lisp +++ b/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) diff --git a/t/simplex.lisp b/t/simplex.lisp index aa042bd..9a46cc8 100644 --- a/t/simplex.lisp +++ b/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))))