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.math.optimization.linear; 019 020 import org.apache.commons.math.optimization.OptimizationException; 021 import org.apache.commons.math.optimization.RealPointValuePair; 022 import org.apache.commons.math.util.MathUtils; 023 024 025 /** 026 * Solves a linear problem using the Two-Phase Simplex Method. 027 * @version $Revision: 797801 $ $Date: 2009-07-25 13:22:06 -0400 (Sat, 25 Jul 2009) $ 028 * @since 2.0 029 */ 030 public class SimplexSolver extends AbstractLinearOptimizer { 031 032 /** Default amount of error to accept in floating point comparisons. */ 033 private static final double DEFAULT_EPSILON = 1.0e-6; 034 035 /** Amount of error to accept in floating point comparisons. */ 036 protected final double epsilon; 037 038 /** 039 * Build a simplex solver with default settings. 040 */ 041 public SimplexSolver() { 042 this(DEFAULT_EPSILON); 043 } 044 045 /** 046 * Build a simplex solver with a specified accepted amount of error 047 * @param epsilon the amount of error to accept in floating point comparisons 048 */ 049 public SimplexSolver(final double epsilon) { 050 this.epsilon = epsilon; 051 } 052 053 /** 054 * Returns the column with the most negative coefficient in the objective function row. 055 * @param tableau simple tableau for the problem 056 * @return column with the most negative coefficient 057 */ 058 private Integer getPivotColumn(SimplexTableau tableau) { 059 double minValue = 0; 060 Integer minPos = null; 061 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { 062 if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) { 063 minValue = tableau.getEntry(0, i); 064 minPos = i; 065 } 066 } 067 return minPos; 068 } 069 070 /** 071 * Returns the row with the minimum ratio as given by the minimum ratio test (MRT). 072 * @param tableau simple tableau for the problem 073 * @param col the column to test the ratio of. See {@link #getPivotColumn(SimplexTableau)} 074 * @return row with the minimum ratio 075 */ 076 private Integer getPivotRow(final int col, final SimplexTableau tableau) { 077 double minRatio = Double.MAX_VALUE; 078 Integer minRatioPos = null; 079 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) { 080 double rhs = tableau.getEntry(i, tableau.getWidth() - 1); 081 if (MathUtils.compareTo(tableau.getEntry(i, col), 0, epsilon) >= 0) { 082 double ratio = rhs / tableau.getEntry(i, col); 083 if (ratio < minRatio) { 084 minRatio = ratio; 085 minRatioPos = i; 086 } 087 } 088 } 089 return minRatioPos; 090 } 091 092 093 /** 094 * Runs one iteration of the Simplex method on the given model. 095 * @param tableau simple tableau for the problem 096 * @throws OptimizationException if the maximal iteration count has been 097 * exceeded or if the model is found not to have a bounded solution 098 */ 099 protected void doIteration(final SimplexTableau tableau) 100 throws OptimizationException { 101 102 incrementIterationsCounter(); 103 104 Integer pivotCol = getPivotColumn(tableau); 105 Integer pivotRow = getPivotRow(pivotCol, tableau); 106 if (pivotRow == null) { 107 throw new UnboundedSolutionException(); 108 } 109 110 // set the pivot element to 1 111 double pivotVal = tableau.getEntry(pivotRow, pivotCol); 112 tableau.divideRow(pivotRow, pivotVal); 113 114 // set the rest of the pivot column to 0 115 for (int i = 0; i < tableau.getHeight(); i++) { 116 if (i != pivotRow) { 117 double multiplier = tableau.getEntry(i, pivotCol); 118 tableau.subtractRow(i, pivotRow, multiplier); 119 } 120 } 121 } 122 123 /** 124 * Checks whether Phase 1 is solved. 125 * @param tableau simple tableau for the problem 126 * @return whether Phase 1 is solved 127 */ 128 private boolean isPhase1Solved(final SimplexTableau tableau) { 129 if (tableau.getNumArtificialVariables() == 0) { 130 return true; 131 } 132 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { 133 if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { 134 return false; 135 } 136 } 137 return true; 138 } 139 140 /** 141 * Returns whether the problem is at an optimal state. 142 * @param tableau simple tableau for the problem 143 * @return whether the model has been solved 144 */ 145 public boolean isOptimal(final SimplexTableau tableau) { 146 if (tableau.getNumArtificialVariables() > 0) { 147 return false; 148 } 149 for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) { 150 if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { 151 return false; 152 } 153 } 154 return true; 155 } 156 157 /** 158 * Solves Phase 1 of the Simplex method. 159 * @param tableau simple tableau for the problem 160 * @exception OptimizationException if the maximal number of iterations is 161 * exceeded, or if the problem is found not to have a bounded solution, or 162 * if there is no feasible solution 163 */ 164 protected void solvePhase1(final SimplexTableau tableau) 165 throws OptimizationException { 166 // make sure we're in Phase 1 167 if (tableau.getNumArtificialVariables() == 0) { 168 return; 169 } 170 171 while (!isPhase1Solved(tableau)) { 172 doIteration(tableau); 173 } 174 175 // if W is not zero then we have no feasible solution 176 if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0, epsilon)) { 177 throw new NoFeasibleSolutionException(); 178 } 179 } 180 181 /** {@inheritDoc} */ 182 @Override 183 public RealPointValuePair doOptimize() 184 throws OptimizationException { 185 final SimplexTableau tableau = 186 new SimplexTableau(f, constraints, goalType, restrictToNonNegative, epsilon); 187 solvePhase1(tableau); 188 tableau.discardArtificialVariables(); 189 while (!isOptimal(tableau)) { 190 doIteration(tableau); 191 } 192 return tableau.getSolution(); 193 } 194 195 }