Home | History | Annotate | Download | only in linear 
      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 
     18 package org.apache.commons.math.optimization.linear;
     19 
     20 import java.util.ArrayList;
     21 import java.util.List;
     22 
     23 import org.apache.commons.math.optimization.OptimizationException;
     24 import org.apache.commons.math.optimization.RealPointValuePair;
     25 import org.apache.commons.math.util.MathUtils;
     26 
     27 
     28 /**
     29  * Solves a linear problem using the Two-Phase Simplex Method.
     30  * @version $Revision: 812831 $ $Date: 2009-09-09 10:48:03 +0200 (mer. 09 sept. 2009) $
     31  * @since 2.0
     32  */
     33 public class SimplexSolver extends AbstractLinearOptimizer {
     34 
     35     /** Default amount of error to accept in floating point comparisons. */
     36     private static final double DEFAULT_EPSILON = 1.0e-6;
     37 
     38     /** Amount of error to accept in floating point comparisons. */
     39     protected final double epsilon;
     40 
     41     /**
     42      * Build a simplex solver with default settings.
     43      */
     44     public SimplexSolver() {
     45         this(DEFAULT_EPSILON);
     46     }
     47 
     48     /**
     49      * Build a simplex solver with a specified accepted amount of error
     50      * @param epsilon the amount of error to accept in floating point comparisons
     51      */
     52     public SimplexSolver(final double epsilon) {
     53         this.epsilon = epsilon;
     54     }
     55 
     56     /**
     57      * Returns the column with the most negative coefficient in the objective function row.
     58      * @param tableau simple tableau for the problem
     59      * @return column with the most negative coefficient
     60      */
     61     private Integer getPivotColumn(SimplexTableau tableau) {
     62         double minValue = 0;
     63         Integer minPos = null;
     64         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
     65             if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) {
     66                 minValue = tableau.getEntry(0, i);
     67                 minPos = i;
     68             }
     69         }
     70         return minPos;
     71     }
     72 
     73     /**
     74      * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
     75      * @param tableau simple tableau for the problem
     76      * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
     77      * @return row with the minimum ratio
     78      */
     79     private Integer getPivotRow(SimplexTableau tableau, final int col) {
     80         // create a list of all the rows that tie for the lowest score in the minimum ratio test
     81         List<Integer> minRatioPositions = new ArrayList<Integer>();
     82         double minRatio = Double.MAX_VALUE;
     83         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
     84             final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
     85             final double entry = tableau.getEntry(i, col);
     86             if (MathUtils.compareTo(entry, 0, epsilon) > 0) {
     87                 final double ratio = rhs / entry;
     88                 if (MathUtils.equals(ratio, minRatio, epsilon)) {
     89                     minRatioPositions.add(i);
     90                 } else if (ratio < minRatio) {
     91                     minRatio = ratio;
     92                     minRatioPositions = new ArrayList<Integer>();
     93                     minRatioPositions.add(i);
     94                 }
     95             }
     96         }
     97 
     98         if (minRatioPositions.size() == 0) {
     99           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 }
    186