Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

[MATH-828] Fixed numerical instabilities in SimplexSolver leading to unbounded solutions. 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1366707 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Thomas Neidhart 2012-07-28 16:37:40 +00:00
parent 6b537b20f9
commit 574c9abbf2
3 changed files with 107 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
src/main/java/org/apache/commons/math3/optimization/linear/SimplexSolver.java
View File

@ -71,7 +71,9 @@ public class SimplexSolver extends AbstractLinearOptimizer {
Integer minPos = null;
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
final double entry = tableau.getEntry(0, i);
if (Precision.compareTo(entry, minValue, maxUlps) < 0) {
// check if the entry is strictly smaller than the current minimum
// do not use a ulp/epsilon check
if (entry < minValue) {
minValue = entry;
minPos = i;
}
@ -95,7 +97,9 @@ public class SimplexSolver extends AbstractLinearOptimizer {
if (Precision.compareTo(entry, 0d, maxUlps) > 0) {
final double ratio = rhs / entry;
final int cmp = Precision.compareTo(ratio, minRatio, maxUlps);
// check if the entry is strictly equal to the current min ratio
// do not use a ulp/epsilon check
final int cmp = Double.compare(ratio, minRatio);
if (cmp == 0) {
minRatioPositions.add(i);
} else if (cmp < 0) {
@ -107,20 +111,40 @@ public class SimplexSolver extends AbstractLinearOptimizer {
}
if (minRatioPositions.size() == 0) {
return null;
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();
final double entry = tableau.getEntry(row, column);
if (Precision.equals(entry, 1d, maxUlps) &&
row.equals(tableau.getBasicRow(column))) {
return row;
}
// there's a degeneracy as indicated by a tie in the minimum ratio test
// 1. 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();
final double entry = tableau.getEntry(row, column);
if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) {
return row;
}
}
}
}
// 2. apply Bland's rule to prevent cycling:
// take the row for which the corresponding basic variable has the smallest index
//
// see http://www.stanford.edu/class/msande310/blandrule.pdf
// see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper)
Integer minRow = null;
int minIndex = tableau.getWidth();
for (Integer row : minRatioPositions) {
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1 && minRow != row; i++) {
if (row == tableau.getBasicRow(i)) {
if (i < minIndex) {
minIndex = i;
minRow = row;
}
}
}
}
return minRow;
}
return minRatioPositions.get(0);
}
@ -149,7 +173,7 @@ public class SimplexSolver extends AbstractLinearOptimizer {
// set the rest of the pivot column to 0
for (int i = 0; i < tableau.getHeight(); i++) {
if (i != pivotRow) {
double multiplier = tableau.getEntry(i, pivotCol);
final double multiplier = tableau.getEntry(i, pivotCol);
tableau.subtractRow(i, pivotRow, multiplier);
}
}

15
src/main/java/org/apache/commons/math3/optimization/linear/SimplexTableau.java
View File

@ -27,6 +27,8 @@ import java.util.HashSet;
import java.util.List;
import java.util.Set;
import odk.lang.FastMath;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
@ -68,6 +70,9 @@ class SimplexTableau implements Serializable {
/** Default amount of error to accept in floating point comparisons (as ulps). */
private static final int DEFAULT_ULPS = 10;
/** The cut-off threshold to zero-out entries. */
private static final double CUTOFF_THRESHOLD = 1e-12;
/** Serializable version identifier. */
private static final long serialVersionUID = -1369660067587938365L;
@ -453,8 +458,14 @@ class SimplexTableau implements Serializable {
*/
protected void subtractRow(final int minuendRow, final int subtrahendRow,
final double multiple) {
tableau.setRowVector(minuendRow, tableau.getRowVector(minuendRow)
.subtract(tableau.getRowVector(subtrahendRow).mapMultiply(multiple)));
for (int i = 0; i < getWidth(); i++) {
double result = tableau.getEntry(minuendRow, i) - tableau.getEntry(subtrahendRow, i) * multiple;
// cut-off values smaller than the CUTOFF_THRESHOLD, otherwise may lead to numerical instabilities
if (FastMath.abs(result) < CUTOFF_THRESHOLD) {
result = 0.0;
}
tableau.setEntry(minuendRow, i, result);
}
}
/**

55
src/test/java/org/apache/commons/math3/optimization/linear/SimplexSolverTest.java
View File

@ -21,6 +21,7 @@ import org.junit.Assert;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
import org.apache.commons.math3.optimization.GoalType;
import org.apache.commons.math3.optimization.PointValuePair;
@ -29,6 +30,27 @@ import org.junit.Test;
public class SimplexSolverTest {
@Test
public void testMath828() {
LinearObjectiveFunction f = new LinearObjectiveFunction(
new double[] { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 0.0);
ArrayList <LinearConstraint>constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] {0.0, 39.0, 23.0, 96.0, 15.0, 48.0, 9.0, 21.0, 48.0, 36.0, 76.0, 19.0, 88.0, 17.0, 16.0, 36.0,}, Relationship.GEQ, 15.0));
constraints.add(new LinearConstraint(new double[] {0.0, 59.0, 93.0, 12.0, 29.0, 78.0, 73.0, 87.0, 32.0, 70.0, 68.0, 24.0, 11.0, 26.0, 65.0, 25.0,}, Relationship.GEQ, 29.0));
constraints.add(new LinearConstraint(new double[] {0.0, 74.0, 5.0, 82.0, 6.0, 97.0, 55.0, 44.0, 52.0, 54.0, 5.0, 93.0, 91.0, 8.0, 20.0, 97.0,}, Relationship.GEQ, 6.0));
constraints.add(new LinearConstraint(new double[] {8.0, -3.0, -28.0, -72.0, -8.0, -31.0, -31.0, -74.0, -47.0, -59.0, -24.0, -57.0, -56.0, -16.0, -92.0, -59.0,}, Relationship.GEQ, 0.0));
constraints.add(new LinearConstraint(new double[] {25.0, -7.0, -99.0, -78.0, -25.0, -14.0, -16.0, -89.0, -39.0, -56.0, -53.0, -9.0, -18.0, -26.0, -11.0, -61.0,}, Relationship.GEQ, 0.0));
constraints.add(new LinearConstraint(new double[] {33.0, -95.0, -15.0, -4.0, -33.0, -3.0, -20.0, -96.0, -27.0, -13.0, -80.0, -24.0, -3.0, -13.0, -57.0, -76.0,}, Relationship.GEQ, 0.0));
constraints.add(new LinearConstraint(new double[] {7.0, -95.0, -39.0, -93.0, -7.0, -94.0, -94.0, -62.0, -76.0, -26.0, -53.0, -57.0, -31.0, -76.0, -53.0, -52.0,}, Relationship.GEQ, 0.0));
double epsilon = 1e-6;
PointValuePair solution = new SimplexSolver().optimize(f, constraints, GoalType.MINIMIZE, true);
Assert.assertEquals(1.0d, solution.getValue(), epsilon);
Assert.assertTrue(validSolution(solution, constraints, epsilon));
}
@Test
public void testMath781() {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 2, 6, 7 }, 0);
@ -560,4 +582,37 @@ public class SimplexSolverTest {
return new LinearConstraint(lhs, relationship, rhs);
}
private static boolean validSolution(PointValuePair solution, List<LinearConstraint> constraints, double epsilon) {
double[] vals = solution.getPoint();
for (LinearConstraint c : constraints) {
double[] coeffs = c.getCoefficients().toArray();
double result = 0.0d;
for (int i = 0; i < vals.length; i++) {
result += vals[i] * coeffs[i];
}
switch (c.getRelationship()) {
case EQ:
if (!Precision.equals(result, c.getValue(), epsilon)) {
return false;
}
break;
case GEQ:
if (Precision.compareTo(result, c.getValue(), epsilon) < 0) {
return false;
}
break;
case LEQ:
if (Precision.compareTo(result, c.getValue(), epsilon) > 0) {
return false;
}
break;
}
}
return true;
}
}