View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.optim.linear;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.apache.commons.math3.exception.TooManyIterationsException;
23  import org.apache.commons.math3.optim.OptimizationData;
24  import org.apache.commons.math3.optim.PointValuePair;
25  import org.apache.commons.math3.util.FastMath;
26  import org.apache.commons.math3.util.Precision;
27  
28  /**
29   * Solves a linear problem using the "Two-Phase Simplex" method.
30   * <p>
31   * The {@link SimplexSolver} supports the following {@link OptimizationData} data provided
32   * as arguments to {@link #optimize(OptimizationData...)}:
33   * <ul>
34   *   <li>objective function: {@link LinearObjectiveFunction} - mandatory</li>
35   *   <li>linear constraints {@link LinearConstraintSet} - mandatory</li>
36   *   <li>type of optimization: {@link org.apache.commons.math3.optim.nonlinear.scalar.GoalType GoalType}
37   *    - optional, default: {@link org.apache.commons.math3.optim.nonlinear.scalar.GoalType#MINIMIZE MINIMIZE}</li>
38   *   <li>whether to allow negative values as solution: {@link NonNegativeConstraint} - optional, default: true</li>
39   *   <li>pivot selection rule: {@link PivotSelectionRule} - optional, default {@link PivotSelectionRule#DANTZIG}</li>
40   *   <li>callback for the best solution: {@link SolutionCallback} - optional</li>
41   *   <li>maximum number of iterations: {@link org.apache.commons.math3.optim.MaxIter} - optional, default: {@link Integer#MAX_VALUE}</li>
42   * </ul>
43   * <p>
44   * <b>Note:</b> Depending on the problem definition, the default convergence criteria
45   * may be too strict, resulting in {@link NoFeasibleSolutionException} or
46   * {@link TooManyIterationsException}. In such a case it is advised to adjust these
47   * criteria with more appropriate values, e.g. relaxing the epsilon value.
48   * <p>
49   * Default convergence criteria:
50   * <ul>
51   *   <li>Algorithm convergence: 1e-6</li>
52   *   <li>Floating-point comparisons: 10 ulp</li>
53   *   <li>Cut-Off value: 1e-10</li>
54    * </ul>
55   * <p>
56   * The cut-off value has been introduced to handle the case of very small pivot elements
57   * in the Simplex tableau, as these may lead to numerical instabilities and degeneracy.
58   * Potential pivot elements smaller than this value will be treated as if they were zero
59   * and are thus not considered by the pivot selection mechanism. The default value is safe
60   * for many problems, but may need to be adjusted in case of very small coefficients
61   * used in either the {@link LinearConstraint} or {@link LinearObjectiveFunction}.
62   *
63   * @version $Id: SimplexSolver.java 1553855 2013-12-28 15:55:24Z erans $
64   * @since 2.0
65   */
66  public class SimplexSolver extends LinearOptimizer {
67      /** Default amount of error to accept in floating point comparisons (as ulps). */
68      static final int DEFAULT_ULPS = 10;
69  
70      /** Default cut-off value. */
71      static final double DEFAULT_CUT_OFF = 1e-10;
72  
73      /** Default amount of error to accept for algorithm convergence. */
74      private static final double DEFAULT_EPSILON = 1.0e-6;
75  
76      /** Amount of error to accept for algorithm convergence. */
77      private final double epsilon;
78  
79      /** Amount of error to accept in floating point comparisons (as ulps). */
80      private final int maxUlps;
81  
82      /**
83       * Cut-off value for entries in the tableau: values smaller than the cut-off
84       * are treated as zero to improve numerical stability.
85       */
86      private final double cutOff;
87  
88      /** The pivot selection method to use. */
89      private PivotSelectionRule pivotSelection;
90  
91      /**
92       * The solution callback to access the best solution found so far in case
93       * the optimizer fails to find an optimal solution within the iteration limits.
94       */
95      private SolutionCallback solutionCallback;
96  
97      /**
98       * Builds a simplex solver with default settings.
99       */
100     public SimplexSolver() {
101         this(DEFAULT_EPSILON, DEFAULT_ULPS, DEFAULT_CUT_OFF);
102     }
103 
104     /**
105      * Builds a simplex solver with a specified accepted amount of error.
106      *
107      * @param epsilon Amount of error to accept for algorithm convergence.
108      */
109     public SimplexSolver(final double epsilon) {
110         this(epsilon, DEFAULT_ULPS, DEFAULT_CUT_OFF);
111     }
112 
113     /**
114      * Builds a simplex solver with a specified accepted amount of error.
115      *
116      * @param epsilon Amount of error to accept for algorithm convergence.
117      * @param maxUlps Amount of error to accept in floating point comparisons.
118      */
119     public SimplexSolver(final double epsilon, final int maxUlps) {
120         this(epsilon, maxUlps, DEFAULT_CUT_OFF);
121     }
122 
123     /**
124      * Builds a simplex solver with a specified accepted amount of error.
125      *
126      * @param epsilon Amount of error to accept for algorithm convergence.
127      * @param maxUlps Amount of error to accept in floating point comparisons.
128      * @param cutOff Values smaller than the cutOff are treated as zero.
129      */
130     public SimplexSolver(final double epsilon, final int maxUlps, final double cutOff) {
131         this.epsilon = epsilon;
132         this.maxUlps = maxUlps;
133         this.cutOff = cutOff;
134         this.pivotSelection = PivotSelectionRule.DANTZIG;
135     }
136 
137     /**
138      * {@inheritDoc}
139      *
140      * @param optData Optimization data. In addition to those documented in
141      * {@link LinearOptimizer#optimize(OptimizationData...)
142      * LinearOptimizer}, this method will register the following data:
143      * <ul>
144      *  <li>{@link SolutionCallback}</li>
145      *  <li>{@link PivotSelectionRule}</li>
146      * </ul>
147      *
148      * @return {@inheritDoc}
149      * @throws TooManyIterationsException if the maximal number of iterations is exceeded.
150      */
151     @Override
152     public PointValuePair optimize(OptimizationData... optData)
153         throws TooManyIterationsException {
154         // Set up base class and perform computation.
155         return super.optimize(optData);
156     }
157 
158     /**
159      * {@inheritDoc}
160      *
161      * @param optData Optimization data.
162      * In addition to those documented in
163      * {@link LinearOptimizer#parseOptimizationData(OptimizationData[])
164      * LinearOptimizer}, this method will register the following data:
165      * <ul>
166      *  <li>{@link SolutionCallback}</li>
167      *  <li>{@link PivotSelectionRule}</li>
168      * </ul>
169      */
170     @Override
171     protected void parseOptimizationData(OptimizationData... optData) {
172         // Allow base class to register its own data.
173         super.parseOptimizationData(optData);
174 
175         // reset the callback before parsing
176         solutionCallback = null;
177 
178         for (OptimizationData data : optData) {
179             if (data instanceof SolutionCallback) {
180                 solutionCallback = (SolutionCallback) data;
181                 continue;
182             }
183             if (data instanceof PivotSelectionRule) {
184                 pivotSelection = (PivotSelectionRule) data;
185                 continue;
186             }
187         }
188     }
189 
190     /**
191      * Returns the column with the most negative coefficient in the objective function row.
192      *
193      * @param tableau Simple tableau for the problem.
194      * @return the column with the most negative coefficient.
195      */
196     private Integer getPivotColumn(SimplexTableau tableau) {
197         double minValue = 0;
198         Integer minPos = null;
199         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
200             final double entry = tableau.getEntry(0, i);
201             // check if the entry is strictly smaller than the current minimum
202             // do not use a ulp/epsilon check
203             if (entry < minValue) {
204                 minValue = entry;
205                 minPos = i;
206 
207                 // Bland's rule: chose the entering column with the lowest index
208                 if (pivotSelection == PivotSelectionRule.BLAND && isValidPivotColumn(tableau, i)) {
209                     break;
210                 }
211             }
212         }
213         return minPos;
214     }
215 
216     /**
217      * Checks whether the given column is valid pivot column, i.e. will result
218      * in a valid pivot row.
219      * <p>
220      * When applying Bland's rule to select the pivot column, it may happen that
221      * there is no corresponding pivot row. This method will check if the selected
222      * pivot column will return a valid pivot row.
223      *
224      * @param tableau simplex tableau for the problem
225      * @param col the column to test
226      * @return {@code true} if the pivot column is valid, {@code false} otherwise
227      */
228     private boolean isValidPivotColumn(SimplexTableau tableau, int col) {
229         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
230             final double entry = tableau.getEntry(i, col);
231 
232             // do the same check as in getPivotRow
233             if (Precision.compareTo(entry, 0d, cutOff) > 0) {
234                 return true;
235             }
236         }
237         return false;
238     }
239 
240     /**
241      * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
242      *
243      * @param tableau Simplex tableau for the problem.
244      * @param col Column to test the ratio of (see {@link #getPivotColumn(SimplexTableau)}).
245      * @return the row with the minimum ratio.
246      */
247     private Integer getPivotRow(SimplexTableau tableau, final int col) {
248         // create a list of all the rows that tie for the lowest score in the minimum ratio test
249         List<Integer> minRatioPositions = new ArrayList<Integer>();
250         double minRatio = Double.MAX_VALUE;
251         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
252             final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
253             final double entry = tableau.getEntry(i, col);
254 
255             // only consider pivot elements larger than the cutOff threshold
256             // selecting others may lead to degeneracy or numerical instabilities
257             if (Precision.compareTo(entry, 0d, cutOff) > 0) {
258                 final double ratio = FastMath.abs(rhs / entry);
259                 // check if the entry is strictly equal to the current min ratio
260                 // do not use a ulp/epsilon check
261                 final int cmp = Double.compare(ratio, minRatio);
262                 if (cmp == 0) {
263                     minRatioPositions.add(i);
264                 } else if (cmp < 0) {
265                     minRatio = ratio;
266                     minRatioPositions.clear();
267                     minRatioPositions.add(i);
268                 }
269             }
270         }
271 
272         if (minRatioPositions.size() == 0) {
273             return null;
274         } else if (minRatioPositions.size() > 1) {
275             // there's a degeneracy as indicated by a tie in the minimum ratio test
276 
277             // 1. check if there's an artificial variable that can be forced out of the basis
278             if (tableau.getNumArtificialVariables() > 0) {
279                 for (Integer row : minRatioPositions) {
280                     for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
281                         int column = i + tableau.getArtificialVariableOffset();
282                         final double entry = tableau.getEntry(row, column);
283                         if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) {
284                             return row;
285                         }
286                     }
287                 }
288             }
289 
290             // 2. apply Bland's rule to prevent cycling:
291             //    take the row for which the corresponding basic variable has the smallest index
292             //
293             // see http://www.stanford.edu/class/msande310/blandrule.pdf
294             // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper)
295 
296             Integer minRow = null;
297             int minIndex = tableau.getWidth();
298             for (Integer row : minRatioPositions) {
299                 final int basicVar = tableau.getBasicVariable(row);
300                 if (basicVar < minIndex) {
301                     minIndex = basicVar;
302                     minRow = row;
303                 }
304             }
305             return minRow;
306         }
307         return minRatioPositions.get(0);
308     }
309 
310     /**
311      * Runs one iteration of the Simplex method on the given model.
312      *
313      * @param tableau Simple tableau for the problem.
314      * @throws TooManyIterationsException if the allowed number of iterations has been exhausted.
315      * @throws UnboundedSolutionException if the model is found not to have a bounded solution.
316      */
317     protected void doIteration(final SimplexTableau tableau)
318         throws TooManyIterationsException,
319                UnboundedSolutionException {
320 
321         incrementIterationCount();
322 
323         Integer pivotCol = getPivotColumn(tableau);
324         Integer pivotRow = getPivotRow(tableau, pivotCol);
325         if (pivotRow == null) {
326             throw new UnboundedSolutionException();
327         }
328 
329         tableau.performRowOperations(pivotCol, pivotRow);
330     }
331 
332     /**
333      * Solves Phase 1 of the Simplex method.
334      *
335      * @param tableau Simple tableau for the problem.
336      * @throws TooManyIterationsException if the allowed number of iterations has been exhausted.
337      * @throws UnboundedSolutionException if the model is found not to have a bounded solution.
338      * @throws NoFeasibleSolutionException if there is no feasible solution?
339      */
340     protected void solvePhase1(final SimplexTableau tableau)
341         throws TooManyIterationsException,
342                UnboundedSolutionException,
343                NoFeasibleSolutionException {
344 
345         // make sure we're in Phase 1
346         if (tableau.getNumArtificialVariables() == 0) {
347             return;
348         }
349 
350         while (!tableau.isOptimal()) {
351             doIteration(tableau);
352         }
353 
354         // if W is not zero then we have no feasible solution
355         if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) {
356             throw new NoFeasibleSolutionException();
357         }
358     }
359 
360     /** {@inheritDoc} */
361     @Override
362     public PointValuePair doOptimize()
363         throws TooManyIterationsException,
364                UnboundedSolutionException,
365                NoFeasibleSolutionException {
366 
367         // reset the tableau to indicate a non-feasible solution in case
368         // we do not pass phase 1 successfully
369         if (solutionCallback != null) {
370             solutionCallback.setTableau(null);
371         }
372 
373         final SimplexTableau tableau =
374             new SimplexTableau(getFunction(),
375                                getConstraints(),
376                                getGoalType(),
377                                isRestrictedToNonNegative(),
378                                epsilon,
379                                maxUlps);
380 
381         solvePhase1(tableau);
382         tableau.dropPhase1Objective();
383 
384         // after phase 1, we are sure to have a feasible solution
385         if (solutionCallback != null) {
386             solutionCallback.setTableau(tableau);
387         }
388 
389         while (!tableau.isOptimal()) {
390             doIteration(tableau);
391         }
392 
393         // check that the solution respects the nonNegative restriction in case
394         // the epsilon/cutOff values are too large for the actual linear problem
395         // (e.g. with very small constraint coefficients), the solver might actually
396         // find a non-valid solution (with negative coefficients).
397         final PointValuePair solution = tableau.getSolution();
398         if (isRestrictedToNonNegative()) {
399             final double[] coeff = solution.getPoint();
400             for (int i = 0; i < coeff.length; i++) {
401                 if (Precision.compareTo(coeff[i], 0, epsilon) < 0) {
402                     throw new NoFeasibleSolutionException();
403                 }
404             }
405         }
406         return solution;
407     }
408 }