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 }