Fixed a OutOfBoundException in simplex solver when some constraints are tight

JIRA: MATH-293

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@813301 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java

@@ -74,6 +74,9 @@ class SimplexTableau implements Serializable {
     /** Whether to restrict the variables to non-negative values. */
     private final boolean restrictToNonNegative;
 
+    /** The variables each column represents */
+    private final List<String> columnLabels = new ArrayList<String>();
+
     /** Simple tableau. */
     private transient RealMatrix tableau;
 
@@ -113,6 +116,27 @@ class SimplexTableau implements Serializable {
         this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                       getConstraintTypeCounts(Relationship.GEQ);
         this.tableau = createTableau(goalType == GoalType.MAXIMIZE);
+        initializeColumnLabels();
+    }
+
+    protected void initializeColumnLabels() {
+      if (getNumObjectiveFunctions() == 2) {
+        columnLabels.add("W");
+      }
+      columnLabels.add("Z");
+      for (int i = 0; i < getOriginalNumDecisionVariables(); i++) {
+        columnLabels.add("x" + i);
+      }
+      if (!restrictToNonNegative) {
+        columnLabels.add("x-");
+      }
+      for (int i = 0; i < getNumSlackVariables(); i++) {
+        columnLabels.add("s" + i);
+      }
+      for (int i = 0; i < getNumArtificialVariables(); i++) {
+        columnLabels.add("a" + i);
+      }
+      columnLabels.add("RHS");
     }
 
     /**
@@ -301,6 +325,10 @@ class SimplexTableau implements Serializable {
           }
         }
 
+        for (int i = columnsToDrop.size() - 1; i >= 0; i--) {
+          columnLabels.remove((int) columnsToDrop.get(i));
+        }
+
         this.tableau = new Array2DRowRealMatrix(matrix);
         this.numArtificialVariables = 0;
     }
@@ -332,12 +360,19 @@ class SimplexTableau implements Serializable {
      * @return current solution
      */
     protected RealPointValuePair getSolution() {
-      double[] coefficients = new double[getOriginalNumDecisionVariables()];
+      int negativeVarColumn = columnLabels.indexOf("x-");
-      Integer negativeVarBasicRow = getBasicRow(getNegativeDecisionVariableOffset());
+      Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null;
       double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset());
+
       Set<Integer> basicRows = new HashSet<Integer>();
+      double[] coefficients = new double[getOriginalNumDecisionVariables()];
       for (int i = 0; i < coefficients.length; i++) {
-          Integer basicRow = getBasicRow(getNumObjectiveFunctions() + i);
+          int colIndex = columnLabels.indexOf("x" + i);
+          if (colIndex < 0) {
+            coefficients[i] = 0;
+            continue;
+          }
+          Integer basicRow = getBasicRow(colIndex);
           if (basicRows.contains(basicRow)) {
               // if multiple variables can take a given value
               // then we choose the first and set the rest equal to 0
@@ -442,15 +477,6 @@ class SimplexTableau implements Serializable {
         return getWidth() - 1;
     }
 
-    /**
-     * Returns the offset of the extra decision variable added when there is a
-     * negative decision variable in the original problem.
-     * @return the offset of x-
-     */
-    protected final int getNegativeDecisionVariableOffset() {
-      return getNumObjectiveFunctions() + getOriginalNumDecisionVariables();
-    }
-
     /**
      * Get the number of decision variables.
      * <p>
@@ -471,7 +497,7 @@ class SimplexTableau implements Serializable {
      * @see #getNumDecisionVariables()
      */
     protected final int getOriginalNumDecisionVariables() {
-        return restrictToNonNegative ? numDecisionVariables : numDecisionVariables - 1;
+        return f.getCoefficients().getDimension();
     }
 
     /**
@@ -562,4 +588,5 @@ class SimplexTableau implements Serializable {
         ois.defaultReadObject();
         MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
     }
+
 }

src/site/xdoc/changes.xml

@@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
   </properties>
   <body>
     <release version="2.1" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-293" due-to="Benjamin McCann">
+        Fixed a OutOfBoundException in simplex solver when some constraints are tight.
+      </action>
       <action dev="luc" type="fix" issue="MATH-291" due-to="Sebb">
         Fixed misleading number formats in error messages for adaptive
         stepsize integrators.

src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java

@@ -116,6 +116,43 @@ public class SimplexSolverTest {
         solver.optimize(f, constraints, GoalType.MINIMIZE, true);
     }
 
+    @Test
+    public void testMath293() throws OptimizationException {
+      LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.8, 0.2, 0.7, 0.3, 0.4, 0.6}, 0 );
+      Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+      constraints.add(new LinearConstraint(new double[] { 1, 0, 1, 0, 1, 0 }, Relationship.EQ, 30.0));
+      constraints.add(new LinearConstraint(new double[] { 0, 1, 0, 1, 0, 1 }, Relationship.EQ, 30.0));
+      constraints.add(new LinearConstraint(new double[] { 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 }, Relationship.GEQ, 10.0));
+      constraints.add(new LinearConstraint(new double[] { 0.0, 0.0, 0.7, 0.3, 0.0, 0.0 }, Relationship.GEQ, 10.0));
+      constraints.add(new LinearConstraint(new double[] { 0.0, 0.0, 0.0, 0.0, 0.4, 0.6 }, Relationship.GEQ, 10.0));
+
+      SimplexSolver solver = new SimplexSolver();
+      RealPointValuePair solution1 = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
+
+      Assert.assertEquals(15.7143, solution1.getPoint()[0], .0001);
+      Assert.assertEquals(0.0, solution1.getPoint()[1], .0001);
+      Assert.assertEquals(14.2857, solution1.getPoint()[2], .0001);
+      Assert.assertEquals(0.0, solution1.getPoint()[3], .0001);
+      Assert.assertEquals(0.0, solution1.getPoint()[4], .0001);
+      Assert.assertEquals(30.0, solution1.getPoint()[5], .0001);
+      Assert.assertEquals(40.57143, solution1.getValue(), .0001);
+
+      double valA = 0.8 * solution1.getPoint()[0] + 0.2 * solution1.getPoint()[1];
+      double valB = 0.7 * solution1.getPoint()[2] + 0.3 * solution1.getPoint()[3];
+      double valC = 0.4 * solution1.getPoint()[4] + 0.6 * solution1.getPoint()[5];
+
+      f = new LinearObjectiveFunction(new double[] { 0.8, 0.2, 0.7, 0.3, 0.4, 0.6}, 0 );
+      constraints = new ArrayList<LinearConstraint>();
+      constraints.add(new LinearConstraint(new double[] { 1, 0, 1, 0, 1, 0 }, Relationship.EQ, 30.0));
+      constraints.add(new LinearConstraint(new double[] { 0, 1, 0, 1, 0, 1 }, Relationship.EQ, 30.0));
+      constraints.add(new LinearConstraint(new double[] { 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 }, Relationship.GEQ, valA));
+      constraints.add(new LinearConstraint(new double[] { 0.0, 0.0, 0.7, 0.3, 0.0, 0.0 }, Relationship.GEQ, valB));
+      constraints.add(new LinearConstraint(new double[] { 0.0, 0.0, 0.0, 0.0, 0.4, 0.6 }, Relationship.GEQ, valC));
+
+      RealPointValuePair solution2 = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
+      Assert.assertEquals(40.57143, solution2.getValue(), .0001);
+    }
+
     @Test
     public void testSimplexSolver() throws OptimizationException {
         LinearObjectiveFunction f =