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