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