SimplexSolver.java

  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.  *      https://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. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */
  21. package org.hipparchus.optim.linear;

  22. import java.util.ArrayList;
  23. import java.util.List;

  24. import org.hipparchus.exception.MathIllegalStateException;
  25. import org.hipparchus.optim.LocalizedOptimFormats;
  26. import org.hipparchus.optim.OptimizationData;
  27. import org.hipparchus.optim.PointValuePair;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.Precision;

  30. /**
  31.  * Solves a linear problem using the "Two-Phase Simplex" method.
  32.  * <p>
  33.  * The {@link SimplexSolver} supports the following {@link OptimizationData} data provided
  34.  * as arguments to {@link #optimize(OptimizationData...)}:
  35.  * <ul>
  36.  *   <li>objective function: {@link LinearObjectiveFunction} - mandatory</li>
  37.  *   <li>linear constraints {@link LinearConstraintSet} - mandatory</li>
  38.  *   <li>type of optimization: {@link org.hipparchus.optim.nonlinear.scalar.GoalType GoalType}
  39.  *    - optional, default: {@link org.hipparchus.optim.nonlinear.scalar.GoalType#MINIMIZE MINIMIZE}</li>
  40.  *   <li>whether to allow negative values as solution: {@link NonNegativeConstraint} - optional, default: true</li>
  41.  *   <li>pivot selection rule: {@link PivotSelectionRule} - optional, default {@link PivotSelectionRule#DANTZIG}</li>
  42.  *   <li>callback for the best solution: {@link SolutionCallback} - optional</li>
  43.  *   <li>maximum number of iterations: {@link org.hipparchus.optim.MaxIter} - optional, default: {@link Integer#MAX_VALUE}</li>
  44.  * </ul>
  45.  * <p>
  46.  * <b>Note:</b> Depending on the problem definition, the default convergence criteria
  47.  * may be too strict, resulting in {@link MathIllegalStateException} or
  48.  * {@link MathIllegalStateException}. In such a case it is advised to adjust these
  49.  * criteria with more appropriate values, e.g. relaxing the epsilon value.
  50.  * <p>
  51.  * Default convergence criteria:
  52.  * <ul>
  53.  *   <li>Algorithm convergence: 1e-6</li>
  54.  *   <li>Floating-point comparisons: 10 ulp</li>
  55.  *   <li>Cut-Off value: 1e-10</li>
  56.   * </ul>
  57.  * <p>
  58.  * The cut-off value has been introduced to handle the case of very small pivot elements
  59.  * in the Simplex tableau, as these may lead to numerical instabilities and degeneracy.
  60.  * Potential pivot elements smaller than this value will be treated as if they were zero
  61.  * and are thus not considered by the pivot selection mechanism. The default value is safe
  62.  * for many problems, but may need to be adjusted in case of very small coefficients
  63.  * used in either the {@link LinearConstraint} or {@link LinearObjectiveFunction}.
  64.  *
  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.     /** Default cut-off value. */
  70.     static final double DEFAULT_CUT_OFF = 1e-10;

  71.     /** Default amount of error to accept for algorithm convergence. */
  72.     private static final double DEFAULT_EPSILON = 1.0e-6;

  73.     /** Amount of error to accept for algorithm convergence. */
  74.     private final double epsilon;

  75.     /** Amount of error to accept in floating point comparisons (as ulps). */
  76.     private final int maxUlps;

  77.     /**
  78.      * Cut-off value for entries in the tableau: values smaller than the cut-off
  79.      * are treated as zero to improve numerical stability.
  80.      */
  81.     private final double cutOff;

  82.     /** The pivot selection method to use. */
  83.     private PivotSelectionRule pivotSelection;

  84.     /**
  85.      * The solution callback to access the best solution found so far in case
  86.      * the optimizer fails to find an optimal solution within the iteration limits.
  87.      */
  88.     private SolutionCallback solutionCallback;

  89.     /**
  90.      * Builds a simplex solver with default settings.
  91.      */
  92.     public SimplexSolver() {
  93.         this(DEFAULT_EPSILON, DEFAULT_ULPS, DEFAULT_CUT_OFF);
  94.     }

  95.     /**
  96.      * Builds a simplex solver with a specified accepted amount of error.
  97.      *
  98.      * @param epsilon Amount of error to accept for algorithm convergence.
  99.      */
  100.     public SimplexSolver(final double epsilon) {
  101.         this(epsilon, DEFAULT_ULPS, DEFAULT_CUT_OFF);
  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.      * @param maxUlps Amount of error to accept in floating point comparisons.
  108.      */
  109.     public SimplexSolver(final double epsilon, final int maxUlps) {
  110.         this(epsilon, maxUlps, DEFAULT_CUT_OFF);
  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.      * @param cutOff Values smaller than the cutOff are treated as zero.
  118.      */
  119.     public SimplexSolver(final double epsilon, final int maxUlps, final double cutOff) {
  120.         this.epsilon = epsilon;
  121.         this.maxUlps = maxUlps;
  122.         this.cutOff = cutOff;
  123.         this.pivotSelection = PivotSelectionRule.DANTZIG;
  124.     }

  125.     /**
  126.      * {@inheritDoc}
  127.      *
  128.      * @param optData Optimization data. In addition to those documented in
  129.      * {@link LinearOptimizer#optimize(OptimizationData...)
  130.      * LinearOptimizer}, this method will register the following data:
  131.      * <ul>
  132.      *  <li>{@link SolutionCallback}</li>
  133.      *  <li>{@link PivotSelectionRule}</li>
  134.      * </ul>
  135.      *
  136.      * @return {@inheritDoc}
  137.      * @throws MathIllegalStateException if the maximal number of iterations is exceeded.
  138.      * @throws org.hipparchus.exception.MathIllegalArgumentException if the dimension
  139.      * of the constraints does not match the dimension of the objective function
  140.      */
  141.     @Override
  142.     public PointValuePair optimize(OptimizationData... optData)
  143.         throws MathIllegalStateException {
  144.         // Set up base class and perform computation.
  145.         return super.optimize(optData);
  146.     }

  147.     /**
  148.      * {@inheritDoc}
  149.      *
  150.      * @param optData Optimization data.
  151.      * In addition to those documented in
  152.      * {@link LinearOptimizer#parseOptimizationData(OptimizationData[])
  153.      * LinearOptimizer}, this method will register the following data:
  154.      * <ul>
  155.      *  <li>{@link SolutionCallback}</li>
  156.      *  <li>{@link PivotSelectionRule}</li>
  157.      * </ul>
  158.      */
  159.     @Override
  160.     protected void parseOptimizationData(OptimizationData... optData) {
  161.         // Allow base class to register its own data.
  162.         super.parseOptimizationData(optData);

  163.         // reset the callback before parsing
  164.         solutionCallback = null;

  165.         for (OptimizationData data : optData) {
  166.             if (data instanceof SolutionCallback) {
  167.                 solutionCallback = (SolutionCallback) data;
  168.                 continue;
  169.             }
  170.             if (data instanceof PivotSelectionRule) {
  171.                 pivotSelection = (PivotSelectionRule) data;
  172.                 continue;
  173.             }
  174.         }
  175.     }

  176.     /**
  177.      * Returns the column with the most negative coefficient in the objective function row.
  178.      *
  179.      * @param tableau Simple tableau for the problem.
  180.      * @return the column with the most negative coefficient.
  181.      */
  182.     private Integer getPivotColumn(SimplexTableau tableau) {
  183.         double minValue = 0;
  184.         Integer minPos = null;
  185.         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
  186.             final double entry = tableau.getEntry(0, i);
  187.             // check if the entry is strictly smaller than the current minimum
  188.             // do not use a ulp/epsilon check
  189.             if (entry < minValue) {
  190.                 minValue = entry;
  191.                 minPos = i;

  192.                 // Bland's rule: chose the entering column with the lowest index
  193.                 if (pivotSelection == PivotSelectionRule.BLAND && isValidPivotColumn(tableau, i)) {
  194.                     break;
  195.                 }
  196.             }
  197.         }
  198.         return minPos;
  199.     }

  200.     /**
  201.      * Checks whether the given column is valid pivot column, i.e. will result
  202.      * in a valid pivot row.
  203.      * <p>
  204.      * When applying Bland's rule to select the pivot column, it may happen that
  205.      * there is no corresponding pivot row. This method will check if the selected
  206.      * pivot column will return a valid pivot row.
  207.      *
  208.      * @param tableau simplex tableau for the problem
  209.      * @param col the column to test
  210.      * @return {@code true} if the pivot column is valid, {@code false} otherwise
  211.      */
  212.     private boolean isValidPivotColumn(SimplexTableau tableau, int col) {
  213.         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
  214.             final double entry = tableau.getEntry(i, col);

  215.             // do the same check as in getPivotRow
  216.             if (Precision.compareTo(entry, 0d, cutOff) > 0) {
  217.                 return true;
  218.             }
  219.         }
  220.         return false;
  221.     }

  222.     /**
  223.      * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
  224.      *
  225.      * @param tableau Simplex tableau for the problem.
  226.      * @param col Column to test the ratio of (see {@link #getPivotColumn(SimplexTableau)}).
  227.      * @return the row with the minimum ratio.
  228.      */
  229.     private Integer getPivotRow(SimplexTableau tableau, final int col) {
  230.         // create a list of all the rows that tie for the lowest score in the minimum ratio test
  231.         List<Integer> minRatioPositions = new ArrayList<>();
  232.         double minRatio = Double.MAX_VALUE;
  233.         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
  234.             final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
  235.             final double entry = tableau.getEntry(i, col);

  236.             // only consider pivot elements larger than the cutOff threshold
  237.             // selecting others may lead to degeneracy or numerical instabilities
  238.             if (Precision.compareTo(entry, 0d, cutOff) > 0) {
  239.                 final double ratio = FastMath.abs(rhs / entry);
  240.                 // check if the entry is strictly equal to the current min ratio
  241.                 // do not use a ulp/epsilon check
  242.                 final int cmp = Double.compare(ratio, minRatio);
  243.                 if (cmp == 0) {
  244.                     minRatioPositions.add(i);
  245.                 } else if (cmp < 0) {
  246.                     minRatio = ratio;
  247.                     minRatioPositions.clear();
  248.                     minRatioPositions.add(i);
  249.                 }
  250.             }
  251.         }

  252.         if (minRatioPositions.isEmpty()) {
  253.             return null;
  254.         } else if (minRatioPositions.size() > 1) {
  255.             // there's a degeneracy as indicated by a tie in the minimum ratio test

  256.             // 1. check if there's an artificial variable that can be forced out of the basis
  257.             if (tableau.getNumArtificialVariables() > 0) {
  258.                 for (Integer row : minRatioPositions) {
  259.                     for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
  260.                         int column = i + tableau.getArtificialVariableOffset();
  261.                         final double entry = tableau.getEntry(row, column);
  262.                         if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) {
  263.                             return row;
  264.                         }
  265.                     }
  266.                 }
  267.             }

  268.             // 2. apply Bland's rule to prevent cycling:
  269.             //    take the row for which the corresponding basic variable has the smallest index
  270.             //
  271.             // see http://www.stanford.edu/class/msande310/blandrule.pdf
  272.             // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper)

  273.             Integer minRow = null;
  274.             int minIndex = tableau.getWidth();
  275.             for (Integer row : minRatioPositions) {
  276.                 final int basicVar = tableau.getBasicVariable(row);
  277.                 if (basicVar < minIndex) {
  278.                     minIndex = basicVar;
  279.                     minRow = row;
  280.                 }
  281.             }
  282.             return minRow;
  283.         }
  284.         return minRatioPositions.get(0);
  285.     }

  286.     /**
  287.      * Runs one iteration of the Simplex method on the given model.
  288.      *
  289.      * @param tableau Simple tableau for the problem.
  290.      * @throws MathIllegalStateException if the allowed number of iterations has been exhausted.
  291.      * @throws MathIllegalStateException if the model is found not to have a bounded solution.
  292.      */
  293.     protected void doIteration(final SimplexTableau tableau)
  294.         throws MathIllegalStateException {

  295.         incrementIterationCount();

  296.         Integer pivotCol = getPivotColumn(tableau);
  297.         Integer pivotRow = getPivotRow(tableau, pivotCol);
  298.         if (pivotRow == null) {
  299.             throw new MathIllegalStateException(LocalizedOptimFormats.UNBOUNDED_SOLUTION);
  300.         }

  301.         tableau.performRowOperations(pivotCol, pivotRow);
  302.     }

  303.     /**
  304.      * Solves Phase 1 of the Simplex method.
  305.      *
  306.      * @param tableau Simple tableau for the problem.
  307.      * @throws MathIllegalStateException if the allowed number of iterations has been exhausted,
  308.      * or if the model is found not to have a bounded solution, or if there is no feasible solution
  309.      */
  310.     protected void solvePhase1(final SimplexTableau tableau)
  311.         throws MathIllegalStateException {

  312.         // make sure we're in Phase 1
  313.         if (tableau.getNumArtificialVariables() == 0) {
  314.             return;
  315.         }

  316.         while (!tableau.isOptimal()) {
  317.             doIteration(tableau);
  318.         }

  319.         // if W is not zero then we have no feasible solution
  320.         if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) {
  321.             throw new MathIllegalStateException(LocalizedOptimFormats.NO_FEASIBLE_SOLUTION);
  322.         }
  323.     }

  324.     /** {@inheritDoc} */
  325.     @Override
  326.     public PointValuePair doOptimize()
  327.         throws MathIllegalStateException {

  328.         // reset the tableau to indicate a non-feasible solution in case
  329.         // we do not pass phase 1 successfully
  330.         if (solutionCallback != null) {
  331.             solutionCallback.setTableau(null);
  332.         }

  333.         final SimplexTableau tableau =
  334.             new SimplexTableau(getFunction(),
  335.                                getConstraints(),
  336.                                getGoalType(),
  337.                                isRestrictedToNonNegative(),
  338.                                epsilon,
  339.                                maxUlps);

  340.         solvePhase1(tableau);
  341.         tableau.dropPhase1Objective();

  342.         // after phase 1, we are sure to have a feasible solution
  343.         if (solutionCallback != null) {
  344.             solutionCallback.setTableau(tableau);
  345.         }

  346.         while (!tableau.isOptimal()) {
  347.             doIteration(tableau);
  348.         }

  349.         // check that the solution respects the nonNegative restriction in case
  350.         // the epsilon/cutOff values are too large for the actual linear problem
  351.         // (e.g. with very small constraint coefficients), the solver might actually
  352.         // find a non-valid solution (with negative coefficients).
  353.         final PointValuePair solution = tableau.getSolution();
  354.         if (isRestrictedToNonNegative()) {
  355.             final double[] coeff = solution.getPoint();
  356.             for (int i = 0; i < coeff.length; i++) {
  357.                 if (Precision.compareTo(coeff[i], 0, epsilon) < 0) {
  358.                     throw new MathIllegalStateException(LocalizedOptimFormats.NO_FEASIBLE_SOLUTION);
  359.                 }
  360.             }
  361.         }
  362.         return solution;
  363.     }
  364. }