From e5a77a7651775ecaa0efe9b479698c0b6ec3d73a Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Tue, 8 Sep 2009 08:40:12 +0000 Subject: [PATCH] added Benjamin's patch from 2009-09-07 JIRA: MATH-286 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@812390 13f79535-47bb-0310-9956-ffa450edef68 --- .../optimization/linear/SimplexSolver.java | 83 +++++++-------- .../optimization/linear/SimplexTableau.java | 28 +++-- .../linear/SimplexSolverTest.java | 100 ++++++++++-------- .../linear/SimplexTableauTest.java | 19 ++-- 4 files changed, 121 insertions(+), 109 deletions(-) diff --git a/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java b/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java index 9a1a24b0e..b7c01e921 100644 --- a/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java +++ b/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java @@ -17,6 +17,9 @@ package org.apache.commons.math.optimization.linear; +import java.util.ArrayList; +import java.util.List; + import org.apache.commons.math.optimization.OptimizationException; import org.apache.commons.math.optimization.RealPointValuePair; import org.apache.commons.math.util.MathUtils; @@ -73,23 +76,42 @@ public class SimplexSolver extends AbstractLinearOptimizer { * @param col the column to test the ratio of. See {@link #getPivotColumn(SimplexTableau)} * @return row with the minimum ratio */ - private Integer getPivotRow(final int col, final SimplexTableau tableau) { + private Integer getPivotRow(SimplexTableau tableau, final int col) { + // create a list of all the rows that tie for the lowest score in the minimum ratio test + List minRatioPositions = new ArrayList(); double minRatio = Double.MAX_VALUE; - Integer minRatioPos = null; for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) { final double rhs = tableau.getEntry(i, tableau.getWidth() - 1); final double entry = tableau.getEntry(i, col); if (MathUtils.compareTo(entry, 0, epsilon) > 0) { final double ratio = rhs / entry; - if (ratio < minRatio) { + if (MathUtils.equals(ratio, minRatio, epsilon)) { + minRatioPositions.add(i); + } else if (ratio < minRatio) { minRatio = ratio; - minRatioPos = i; + minRatioPositions = new ArrayList(); + minRatioPositions.add(i); } } } - return minRatioPos; - } + if (minRatioPositions.size() == 0) { + return null; + } else if (minRatioPositions.size() > 1) { + // there's a degeneracy as indicated by a tie in the minimum ratio test + // check if there's an artificial variable that can be forced out of the basis + for (Integer row : minRatioPositions) { + for (int i = 0; i < tableau.getNumArtificialVariables(); i++) { + int column = i + tableau.getArtificialVariableOffset(); + if (MathUtils.equals(tableau.getEntry(row, column), 1, epsilon) && + row.equals(tableau.getBasicRow(column))) { + return row; + } + } + } + } + return minRatioPositions.get(0); + } /** * Runs one iteration of the Simplex method on the given model. @@ -103,7 +125,7 @@ public class SimplexSolver extends AbstractLinearOptimizer { incrementIterationsCounter(); Integer pivotCol = getPivotColumn(tableau); - Integer pivotRow = getPivotRow(pivotCol, tableau); + Integer pivotRow = getPivotRow(tableau, pivotCol); if (pivotRow == null) { throw new UnboundedSolutionException(); } @@ -121,40 +143,6 @@ public class SimplexSolver extends AbstractLinearOptimizer { } } - /** - * Checks whether Phase 1 is solved. - * @param tableau simple tableau for the problem - * @return whether Phase 1 is solved - */ - private boolean isPhase1Solved(final SimplexTableau tableau) { - if (tableau.getNumArtificialVariables() == 0) { - return true; - } - for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { - if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { - return false; - } - } - return true; - } - - /** - * Returns whether the problem is at an optimal state. - * @param tableau simple tableau for the problem - * @return whether the model has been solved - */ - public boolean isOptimal(final SimplexTableau tableau) { - if (tableau.getNumArtificialVariables() > 0) { - return false; - } - for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { - if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { - return false; - } - } - return true; - } - /** * Solves Phase 1 of the Simplex method. * @param tableau simple tableau for the problem @@ -162,14 +150,14 @@ public class SimplexSolver extends AbstractLinearOptimizer { * exceeded, or if the problem is found not to have a bounded solution, or * if there is no feasible solution */ - protected void solvePhase1(final SimplexTableau tableau) - throws OptimizationException { + protected void solvePhase1(final SimplexTableau tableau) throws OptimizationException { + // make sure we're in Phase 1 if (tableau.getNumArtificialVariables() == 0) { return; } - while (!isPhase1Solved(tableau)) { + while (!tableau.isOptimal()) { doIteration(tableau); } @@ -181,13 +169,14 @@ public class SimplexSolver extends AbstractLinearOptimizer { /** {@inheritDoc} */ @Override - public RealPointValuePair doOptimize() - throws OptimizationException { + public RealPointValuePair doOptimize() throws OptimizationException { final SimplexTableau tableau = new SimplexTableau(function, linearConstraints, goal, nonNegative, epsilon); + solvePhase1(tableau); tableau.discardArtificialVariables(); - while (!isOptimal(tableau)) { + + while (!tableau.isOptimal()) { doIteration(tableau); } return tableau.getSolution(); diff --git a/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java b/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java index 3a5edfdaa..5f419c129 100644 --- a/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java +++ b/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java @@ -27,9 +27,9 @@ import java.util.HashSet; import java.util.List; import java.util.Set; +import org.apache.commons.math.linear.Array2DRowRealMatrix; import org.apache.commons.math.linear.MatrixUtils; import org.apache.commons.math.linear.RealMatrix; -import org.apache.commons.math.linear.Array2DRowRealMatrix; import org.apache.commons.math.linear.RealVector; import org.apache.commons.math.optimization.GoalType; import org.apache.commons.math.optimization.RealPointValuePair; @@ -106,7 +106,8 @@ class SimplexTableau implements Serializable { this.constraints = normalizeConstraints(constraints); this.restrictToNonNegative = restrictToNonNegative; this.epsilon = epsilon; - this.numDecisionVariables = getNumVariables() + (restrictToNonNegative ? 0 : 1); + this.numDecisionVariables = f.getCoefficients().getDimension() + + (restrictToNonNegative ? 0 : 1); this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) + getConstraintTypeCounts(Relationship.GEQ); this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) + @@ -183,13 +184,6 @@ class SimplexTableau implements Serializable { } - /** Get the number of variables. - * @return number of variables - */ - public int getNumVariables() { - return f.getCoefficients().getDimension(); - } - /** * Get new versions of the constraints which have positive right hand sides. * @param originalConstraints original (not normalized) constraints @@ -270,7 +264,7 @@ class SimplexTableau implements Serializable { * @param col index of the column to check * @return the row that the variable is basic in. null if the column is not basic */ - private Integer getBasicRow(final int col) { + Integer getBasicRow(final int col) { return getBasicRow(col, true); } @@ -323,7 +317,6 @@ class SimplexTableau implements Serializable { this.numArtificialVariables = 0; } - /** * @param src the source array * @param dest the destination array @@ -332,6 +325,19 @@ class SimplexTableau implements Serializable { System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length); } + /** + * Returns whether the problem is at an optimal state. + * @return whether the model has been solved + */ + boolean isOptimal() { + for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) { + if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { + return false; + } + } + return true; + } + /** * Get the current solution. * diff --git a/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java b/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java index 584d42d50..3e0ce2bfe 100644 --- a/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java +++ b/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java @@ -17,13 +17,11 @@ package org.apache.commons.math.optimization.linear; -import static org.junit.Assert.assertEquals; +import org.junit.Assert; import java.util.ArrayList; import java.util.Collection; -import org.apache.commons.math.linear.RealVector; -import org.apache.commons.math.linear.ArrayRealVector; import org.apache.commons.math.optimization.GoalType; import org.apache.commons.math.optimization.OptimizationException; import org.apache.commons.math.optimization.RealPointValuePair; @@ -42,20 +40,34 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true); - assertEquals(0.0, solution.getPoint()[0], .0000001); - assertEquals(1.0, solution.getPoint()[1], .0000001); - assertEquals(1.0, solution.getPoint()[2], .0000001); - assertEquals(3.0, solution.getValue(), .0000001); + Assert.assertEquals(0.0, solution.getPoint()[0], .0000001); + Assert.assertEquals(1.0, solution.getPoint()[1], .0000001); + Assert.assertEquals(1.0, solution.getPoint()[2], .0000001); + Assert.assertEquals(3.0, solution.getValue(), .0000001); } @Test public void testMath286() throws OptimizationException { - LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.2, 0.3 }, 0 ); - Collection constraints = new ArrayList(); - constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 23.0)); + LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.2, 0.3 }, 0 ); + Collection constraints = new ArrayList(); + constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 23.0)); - RealPointValuePair solution = new SimplexSolver().optimize(f, constraints, GoalType.MAXIMIZE, true); - assertEquals(6.9, solution.getValue(), .0000001); + SimplexSolver solver = new SimplexSolver(); + RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true); + Assert.assertEquals(6.9, solution.getValue(), .0000001); + } + + @Test + public void testDegeneracy() throws OptimizationException { + LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.8, 0.7 }, 0 ); + Collection constraints = new ArrayList(); + constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.LEQ, 18.0)); + constraints.add(new LinearConstraint(new double[] { 1, 0 }, Relationship.GEQ, 10.0)); + constraints.add(new LinearConstraint(new double[] { 0, 1 }, Relationship.GEQ, 8.0)); + + SimplexSolver solver = new SimplexSolver(); + RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true); + Assert.assertEquals(13.6, solution.getValue(), .0000001); } @Test @@ -70,7 +82,7 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true); - assertEquals(10.0, solution.getValue(), .0000001); + Assert.assertEquals(10.0, solution.getValue(), .0000001); } @Test @@ -80,9 +92,9 @@ public class SimplexSolverTest { constraints.add(new LinearConstraint(new double[] { 2, 0 }, Relationship.GEQ, -1.0)); SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true); - assertEquals(0, solution.getValue(), .0000001); - assertEquals(0, solution.getPoint()[0], .0000001); - assertEquals(0, solution.getPoint()[1], .0000001); + Assert.assertEquals(0, solution.getValue(), .0000001); + Assert.assertEquals(0, solution.getPoint()[0], .0000001); + Assert.assertEquals(0, solution.getPoint()[1], .0000001); } @Test(expected=NoFeasibleSolutionException.class) @@ -105,9 +117,9 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false); - assertEquals(2.0, solution.getPoint()[0], 0.0); - assertEquals(2.0, solution.getPoint()[1], 0.0); - assertEquals(57.0, solution.getValue(), 0.0); + Assert.assertEquals(2.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(2.0, solution.getPoint()[1], 0.0); + Assert.assertEquals(57.0, solution.getValue(), 0.0); } @Test @@ -118,8 +130,8 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false); - assertEquals(10.0, solution.getPoint()[0], 0.0); - assertEquals(30.0, solution.getValue(), 0.0); + Assert.assertEquals(10.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(30.0, solution.getValue(), 0.0); } /** @@ -136,9 +148,9 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false); - assertEquals(2.0, solution.getPoint()[0], 0.0); - assertEquals(2.0, solution.getPoint()[1], 0.0); - assertEquals(50.0, solution.getValue(), 0.0); + Assert.assertEquals(2.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(2.0, solution.getPoint()[1], 0.0); + Assert.assertEquals(50.0, solution.getValue(), 0.0); } @Test @@ -151,9 +163,9 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, false); - assertEquals(4.0, solution.getPoint()[0], 0.0); - assertEquals(0.0, solution.getPoint()[1], 0.0); - assertEquals(-13.0, solution.getValue(), 0.0); + Assert.assertEquals(4.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(0.0, solution.getPoint()[1], 0.0); + Assert.assertEquals(-13.0, solution.getValue(), 0.0); } @Test @@ -165,9 +177,9 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false); - assertEquals(-2.0, solution.getPoint()[0], 0.0); - assertEquals(8.0, solution.getPoint()[1], 0.0); - assertEquals(12.0, solution.getValue(), 0.0); + Assert.assertEquals(-2.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(8.0, solution.getPoint()[1], 0.0); + Assert.assertEquals(12.0, solution.getValue(), 0.0); } @Test(expected = NoFeasibleSolutionException.class) @@ -203,12 +215,12 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true); - assertEquals(2902.92783505155, solution.getPoint()[0], .0000001); - assertEquals(480.419243986254, solution.getPoint()[1], .0000001); - assertEquals(0.0, solution.getPoint()[2], .0000001); - assertEquals(0.0, solution.getPoint()[3], .0000001); - assertEquals(0.0, solution.getPoint()[4], .0000001); - assertEquals(1438556.7491409, solution.getValue(), .0000001); + Assert.assertEquals(2902.92783505155, solution.getPoint()[0], .0000001); + Assert.assertEquals(480.419243986254, solution.getPoint()[1], .0000001); + Assert.assertEquals(0.0, solution.getPoint()[2], .0000001); + Assert.assertEquals(0.0, solution.getPoint()[3], .0000001); + Assert.assertEquals(0.0, solution.getPoint()[4], .0000001); + Assert.assertEquals(1438556.7491409, solution.getValue(), .0000001); } @Test @@ -222,10 +234,10 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false); - assertEquals(1.0, solution.getPoint()[0], 0.0); - assertEquals(1.0, solution.getPoint()[1], 0.0); - assertEquals(0.0, solution.getPoint()[2], 0.0); - assertEquals(15.0, solution.getValue(), 0.0); + Assert.assertEquals(1.0, solution.getPoint()[0], 0.0); + Assert.assertEquals(1.0, solution.getPoint()[1], 0.0); + Assert.assertEquals(0.0, solution.getPoint()[2], 0.0); + Assert.assertEquals(15.0, solution.getValue(), 0.0); } @Test @@ -236,7 +248,7 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true); - assertEquals(0, solution.getValue(), .0000001); + Assert.assertEquals(0, solution.getValue(), .0000001); } @Test @@ -363,7 +375,7 @@ public class SimplexSolverTest { SimplexSolver solver = new SimplexSolver(); RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true); - assertEquals(7518.0, solution.getValue(), .0000001); + Assert.assertEquals(7518.0, solution.getValue(), .0000001); } /** @@ -385,13 +397,13 @@ public class SimplexSolverTest { String[] equationParts = s.split("[>|<]?="); double rhs = Double.parseDouble(equationParts[1].trim()); - RealVector lhs = new ArrayRealVector(numCoefficients); + double[] lhs = new double[numCoefficients]; String left = equationParts[0].replaceAll(" ?x", ""); String[] coefficients = left.split(" "); for (String coefficient : coefficients) { double value = coefficient.charAt(0) == '-' ? -1 : 1; int index = Integer.parseInt(coefficient.replaceFirst("[+|-]", "").trim()); - lhs.setEntry(index, value); + lhs[index] = value; } return new LinearConstraint(lhs, relationship, rhs); } diff --git a/src/test/java/org/apache/commons/math/optimization/linear/SimplexTableauTest.java b/src/test/java/org/apache/commons/math/optimization/linear/SimplexTableauTest.java index f45fcc72a..0c8edd283 100644 --- a/src/test/java/org/apache/commons/math/optimization/linear/SimplexTableauTest.java +++ b/src/test/java/org/apache/commons/math/optimization/linear/SimplexTableauTest.java @@ -22,11 +22,12 @@ import java.util.Collection; import org.apache.commons.math.TestUtils; import org.apache.commons.math.optimization.GoalType; +import org.junit.Assert; +import org.junit.Test; -import junit.framework.TestCase; - -public class SimplexTableauTest extends TestCase { +public class SimplexTableauTest { + @Test public void testInitialization() { LinearObjectiveFunction f = createFunction(); Collection constraints = createConstraints(); @@ -42,6 +43,7 @@ public class SimplexTableauTest extends TestCase { assertMatrixEquals(expectedInitialTableau, tableau.getData()); } + @Test public void testdiscardArtificialVariables() { LinearObjectiveFunction f = createFunction(); Collection constraints = createConstraints(); @@ -57,6 +59,7 @@ public class SimplexTableauTest extends TestCase { assertMatrixEquals(expectedTableau, tableau.getData()); } + @Test public void testTableauWithNoArtificialVars() { LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {15, 10}, 0); Collection constraints = new ArrayList(); @@ -74,13 +77,15 @@ public class SimplexTableauTest extends TestCase { assertMatrixEquals(initialTableau, tableau.getData()); } + @Test public void testSerial() { LinearObjectiveFunction f = createFunction(); Collection constraints = createConstraints(); SimplexTableau tableau = new SimplexTableau(f, constraints, GoalType.MAXIMIZE, false, 1.0e-6); - assertEquals(tableau, TestUtils.serializeAndRecover(tableau)); + Assert.assertEquals(tableau, TestUtils.serializeAndRecover(tableau)); } + private LinearObjectiveFunction createFunction() { return new LinearObjectiveFunction(new double[] {15, 10}, 0); } @@ -94,11 +99,11 @@ public class SimplexTableauTest extends TestCase { } private void assertMatrixEquals(double[][] expected, double[][] result) { - assertEquals("Wrong number of rows.", expected.length, result.length); + Assert.assertEquals("Wrong number of rows.", expected.length, result.length); for (int i = 0; i < expected.length; i++) { - assertEquals("Wrong number of columns.", expected[i].length, result[i].length); + Assert.assertEquals("Wrong number of columns.", expected[i].length, result[i].length); for (int j = 0; j < expected[i].length; j++) { - assertEquals("Wrong value at position [" + i + "," + j + "]", expected[i][j], result[i][j]); + Assert.assertEquals("Wrong value at position [" + i + "," + j + "]", expected[i][j], result[i][j], 1.0e-15); } } }