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 018 package org.apache.commons.math3.optimization.linear; 019 020 import java.util.ArrayList; 021 import java.util.List; 022 023 import org.apache.commons.math3.exception.MaxCountExceededException; 024 import org.apache.commons.math3.optimization.PointValuePair; 025 import org.apache.commons.math3.util.Precision; 026 027 028 /** 029 * Solves a linear problem using the Two-Phase Simplex Method. 030 * 031 * @version $Id: SimplexSolver.java 1422230 2012-12-15 12:11:13Z erans $ 032 * @deprecated As of 3.1 (to be removed in 4.0). 033 * @since 2.0 034 */ 035 @Deprecated 036 public class SimplexSolver extends AbstractLinearOptimizer { 037 038 /** Default amount of error to accept for algorithm convergence. */ 039 private static final double DEFAULT_EPSILON = 1.0e-6; 040 041 /** Default amount of error to accept in floating point comparisons (as ulps). */ 042 private static final int DEFAULT_ULPS = 10; 043 044 /** Amount of error to accept for algorithm convergence. */ 045 private final double epsilon; 046 047 /** Amount of error to accept in floating point comparisons (as ulps). */ 048 private final int maxUlps; 049 050 /** 051 * Build a simplex solver with default settings. 052 */ 053 public SimplexSolver() { 054 this(DEFAULT_EPSILON, DEFAULT_ULPS); 055 } 056 057 /** 058 * Build a simplex solver with a specified accepted amount of error 059 * @param epsilon the amount of error to accept for algorithm convergence 060 * @param maxUlps amount of error to accept in floating point comparisons 061 */ 062 public SimplexSolver(final double epsilon, final int maxUlps) { 063 this.epsilon = epsilon; 064 this.maxUlps = maxUlps; 065 } 066 067 /** 068 * Returns the column with the most negative coefficient in the objective function row. 069 * @param tableau simple tableau for the problem 070 * @return column with the most negative coefficient 071 */ 072 private Integer getPivotColumn(SimplexTableau tableau) { 073 double minValue = 0; 074 Integer minPos = null; 075 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { 076 final double entry = tableau.getEntry(0, i); 077 // check if the entry is strictly smaller than the current minimum 078 // do not use a ulp/epsilon check 079 if (entry < minValue) { 080 minValue = entry; 081 minPos = i; 082 } 083 } 084 return minPos; 085 } 086 087 /** 088 * Returns the row with the minimum ratio as given by the minimum ratio test (MRT). 089 * @param tableau simple tableau for the problem 090 * @param col the column to test the ratio of. See {@link #getPivotColumn(SimplexTableau)} 091 * @return row with the minimum ratio 092 */ 093 private Integer getPivotRow(SimplexTableau tableau, final int col) { 094 // create a list of all the rows that tie for the lowest score in the minimum ratio test 095 List<Integer> minRatioPositions = new ArrayList<Integer>(); 096 double minRatio = Double.MAX_VALUE; 097 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) { 098 final double rhs = tableau.getEntry(i, tableau.getWidth() - 1); 099 final double entry = tableau.getEntry(i, col); 100 101 if (Precision.compareTo(entry, 0d, maxUlps) > 0) { 102 final double ratio = rhs / entry; 103 // check if the entry is strictly equal to the current min ratio 104 // do not use a ulp/epsilon check 105 final int cmp = Double.compare(ratio, minRatio); 106 if (cmp == 0) { 107 minRatioPositions.add(i); 108 } else if (cmp < 0) { 109 minRatio = ratio; 110 minRatioPositions = new ArrayList<Integer>(); 111 minRatioPositions.add(i); 112 } 113 } 114 } 115 116 if (minRatioPositions.size() == 0) { 117 return null; 118 } else if (minRatioPositions.size() > 1) { 119 // there's a degeneracy as indicated by a tie in the minimum ratio test 120 121 // 1. check if there's an artificial variable that can be forced out of the basis 122 if (tableau.getNumArtificialVariables() > 0) { 123 for (Integer row : minRatioPositions) { 124 for (int i = 0; i < tableau.getNumArtificialVariables(); i++) { 125 int column = i + tableau.getArtificialVariableOffset(); 126 final double entry = tableau.getEntry(row, column); 127 if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) { 128 return row; 129 } 130 } 131 } 132 } 133 134 // 2. apply Bland's rule to prevent cycling: 135 // take the row for which the corresponding basic variable has the smallest index 136 // 137 // see http://www.stanford.edu/class/msande310/blandrule.pdf 138 // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper) 139 // 140 // Additional heuristic: if we did not get a solution after half of maxIterations 141 // revert to the simple case of just returning the top-most row 142 // This heuristic is based on empirical data gathered while investigating MATH-828. 143 if (getIterations() < getMaxIterations() / 2) { 144 Integer minRow = null; 145 int minIndex = tableau.getWidth(); 146 final int varStart = tableau.getNumObjectiveFunctions(); 147 final int varEnd = tableau.getWidth() - 1; 148 for (Integer row : minRatioPositions) { 149 for (int i = varStart; i < varEnd && !row.equals(minRow); i++) { 150 final Integer basicRow = tableau.getBasicRow(i); 151 if (basicRow != null && basicRow.equals(row)) { 152 if (i < minIndex) { 153 minIndex = i; 154 minRow = row; 155 } 156 } 157 } 158 } 159 return minRow; 160 } 161 } 162 return minRatioPositions.get(0); 163 } 164 165 /** 166 * Runs one iteration of the Simplex method on the given model. 167 * @param tableau simple tableau for the problem 168 * @throws MaxCountExceededException if the maximal iteration count has been exceeded 169 * @throws UnboundedSolutionException if the model is found not to have a bounded solution 170 */ 171 protected void doIteration(final SimplexTableau tableau) 172 throws MaxCountExceededException, UnboundedSolutionException { 173 174 incrementIterationsCounter(); 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 * @param tableau simple tableau for the problem 198 * @throws MaxCountExceededException if the maximal iteration count has been exceeded 199 * @throws UnboundedSolutionException if the model is found not to have a bounded solution 200 * @throws NoFeasibleSolutionException if there is no feasible solution 201 */ 202 protected void solvePhase1(final SimplexTableau tableau) 203 throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException { 204 205 // make sure we're in Phase 1 206 if (tableau.getNumArtificialVariables() == 0) { 207 return; 208 } 209 210 while (!tableau.isOptimal()) { 211 doIteration(tableau); 212 } 213 214 // if W is not zero then we have no feasible solution 215 if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) { 216 throw new NoFeasibleSolutionException(); 217 } 218 } 219 220 /** {@inheritDoc} */ 221 @Override 222 public PointValuePair doOptimize() 223 throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException { 224 final SimplexTableau tableau = 225 new SimplexTableau(getFunction(), 226 getConstraints(), 227 getGoalType(), 228 restrictToNonNegative(), 229 epsilon, 230 maxUlps); 231 232 solvePhase1(tableau); 233 tableau.dropPhase1Objective(); 234 235 while (!tableau.isOptimal()) { 236 doIteration(tableau); 237 } 238 return tableau.getSolution(); 239 } 240 241 }