001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math3.optim.linear; 018 019 import java.util.ArrayList; 020 import java.util.List; 021 import org.apache.commons.math3.exception.TooManyIterationsException; 022 import org.apache.commons.math3.optim.PointValuePair; 023 import org.apache.commons.math3.util.Precision; 024 025 /** 026 * Solves a linear problem using the "Two-Phase Simplex" method. 027 * 028 * @version $Id: SimplexSolver.java 1416643 2012-12-03 19:37:14Z tn $ 029 * @since 2.0 030 */ 031 public class SimplexSolver extends LinearOptimizer { 032 /** Default amount of error to accept for algorithm convergence. */ 033 private static final double DEFAULT_EPSILON = 1.0e-6; 034 035 /** Default amount of error to accept in floating point comparisons (as ulps). */ 036 private static final int DEFAULT_ULPS = 10; 037 038 /** Amount of error to accept for algorithm convergence. */ 039 private final double epsilon; 040 041 /** Amount of error to accept in floating point comparisons (as ulps). */ 042 private final int maxUlps; 043 044 /** 045 * Builds a simplex solver with default settings. 046 */ 047 public SimplexSolver() { 048 this(DEFAULT_EPSILON, DEFAULT_ULPS); 049 } 050 051 /** 052 * Builds a simplex solver with a specified accepted amount of error. 053 * 054 * @param epsilon Amount of error to accept for algorithm convergence. 055 * @param maxUlps Amount of error to accept in floating point comparisons. 056 */ 057 public SimplexSolver(final double epsilon, 058 final int maxUlps) { 059 this.epsilon = epsilon; 060 this.maxUlps = maxUlps; 061 } 062 063 /** 064 * Returns the column with the most negative coefficient in the objective function row. 065 * 066 * @param tableau Simple tableau for the problem. 067 * @return the column with the most negative coefficient. 068 */ 069 private Integer getPivotColumn(SimplexTableau tableau) { 070 double minValue = 0; 071 Integer minPos = null; 072 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { 073 final double entry = tableau.getEntry(0, i); 074 // check if the entry is strictly smaller than the current minimum 075 // do not use a ulp/epsilon check 076 if (entry < minValue) { 077 minValue = entry; 078 minPos = i; 079 } 080 } 081 return minPos; 082 } 083 084 /** 085 * Returns the row with the minimum ratio as given by the minimum ratio test (MRT). 086 * 087 * @param tableau Simple tableau for the problem. 088 * @param col Column to test the ratio of (see {@link #getPivotColumn(SimplexTableau)}). 089 * @return the row with the minimum ratio. 090 */ 091 private Integer getPivotRow(SimplexTableau tableau, final int col) { 092 // create a list of all the rows that tie for the lowest score in the minimum ratio test 093 List<Integer> minRatioPositions = new ArrayList<Integer>(); 094 double minRatio = Double.MAX_VALUE; 095 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) { 096 final double rhs = tableau.getEntry(i, tableau.getWidth() - 1); 097 final double entry = tableau.getEntry(i, col); 098 099 if (Precision.compareTo(entry, 0d, maxUlps) > 0) { 100 final double ratio = rhs / entry; 101 // check if the entry is strictly equal to the current min ratio 102 // do not use a ulp/epsilon check 103 final int cmp = Double.compare(ratio, minRatio); 104 if (cmp == 0) { 105 minRatioPositions.add(i); 106 } else if (cmp < 0) { 107 minRatio = ratio; 108 minRatioPositions = new ArrayList<Integer>(); 109 minRatioPositions.add(i); 110 } 111 } 112 } 113 114 if (minRatioPositions.size() == 0) { 115 return null; 116 } else if (minRatioPositions.size() > 1) { 117 // there's a degeneracy as indicated by a tie in the minimum ratio test 118 119 // 1. check if there's an artificial variable that can be forced out of the basis 120 if (tableau.getNumArtificialVariables() > 0) { 121 for (Integer row : minRatioPositions) { 122 for (int i = 0; i < tableau.getNumArtificialVariables(); i++) { 123 int column = i + tableau.getArtificialVariableOffset(); 124 final double entry = tableau.getEntry(row, column); 125 if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) { 126 return row; 127 } 128 } 129 } 130 } 131 132 // 2. apply Bland's rule to prevent cycling: 133 // take the row for which the corresponding basic variable has the smallest index 134 // 135 // see http://www.stanford.edu/class/msande310/blandrule.pdf 136 // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper) 137 // 138 // Additional heuristic: if we did not get a solution after half of maxIterations 139 // revert to the simple case of just returning the top-most row 140 // This heuristic is based on empirical data gathered while investigating MATH-828. 141 if (getEvaluations() < getMaxEvaluations() / 2) { 142 Integer minRow = null; 143 int minIndex = tableau.getWidth(); 144 final int varStart = tableau.getNumObjectiveFunctions(); 145 final int varEnd = tableau.getWidth() - 1; 146 for (Integer row : minRatioPositions) { 147 for (int i = varStart; i < varEnd && !row.equals(minRow); i++) { 148 final Integer basicRow = tableau.getBasicRow(i); 149 if (basicRow != null && basicRow.equals(row)) { 150 if (i < minIndex) { 151 minIndex = i; 152 minRow = row; 153 } 154 } 155 } 156 } 157 return minRow; 158 } 159 } 160 return minRatioPositions.get(0); 161 } 162 163 /** 164 * Runs one iteration of the Simplex method on the given model. 165 * 166 * @param tableau Simple tableau for the problem. 167 * @throws TooManyIterationsException if the allowed number of iterations has been exhausted. 168 * @throws UnboundedSolutionException if the model is found not to have a bounded solution. 169 */ 170 protected void doIteration(final SimplexTableau tableau) 171 throws TooManyIterationsException, 172 UnboundedSolutionException { 173 174 incrementIterationCount(); 175 176 Integer pivotCol = getPivotColumn(tableau); 177 Integer pivotRow = getPivotRow(tableau, pivotCol); 178 if (pivotRow == null) { 179 throw new UnboundedSolutionException(); 180 } 181 182 // set the pivot element to 1 183 double pivotVal = tableau.getEntry(pivotRow, pivotCol); 184 tableau.divideRow(pivotRow, pivotVal); 185 186 // set the rest of the pivot column to 0 187 for (int i = 0; i < tableau.getHeight(); i++) { 188 if (i != pivotRow) { 189 final double multiplier = tableau.getEntry(i, pivotCol); 190 tableau.subtractRow(i, pivotRow, multiplier); 191 } 192 } 193 } 194 195 /** 196 * Solves Phase 1 of the Simplex method. 197 * 198 * @param tableau Simple tableau for the problem. 199 * @throws TooManyIterationsException if the allowed number of iterations has been exhausted. 200 * @throws UnboundedSolutionException if the model is found not to have a bounded solution. 201 * @throws NoFeasibleSolutionException if there is no feasible solution? 202 */ 203 protected void solvePhase1(final SimplexTableau tableau) 204 throws TooManyIterationsException, 205 UnboundedSolutionException, 206 NoFeasibleSolutionException { 207 208 // make sure we're in Phase 1 209 if (tableau.getNumArtificialVariables() == 0) { 210 return; 211 } 212 213 while (!tableau.isOptimal()) { 214 doIteration(tableau); 215 } 216 217 // if W is not zero then we have no feasible solution 218 if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) { 219 throw new NoFeasibleSolutionException(); 220 } 221 } 222 223 /** {@inheritDoc} */ 224 @Override 225 public PointValuePair doOptimize() 226 throws TooManyIterationsException, 227 UnboundedSolutionException, 228 NoFeasibleSolutionException { 229 final SimplexTableau tableau = 230 new SimplexTableau(getFunction(), 231 getConstraints(), 232 getGoalType(), 233 isRestrictedToNonNegative(), 234 epsilon, 235 maxUlps); 236 237 solvePhase1(tableau); 238 tableau.dropPhase1Objective(); 239 240 while (!tableau.isOptimal()) { 241 doIteration(tableau); 242 } 243 return tableau.getSolution(); 244 } 245 }