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  
18  package org.apache.commons.math.optimization.linear;
19  
20  import org.apache.commons.math.optimization.OptimizationException;
21  import org.apache.commons.math.optimization.RealPointValuePair;
22  import org.apache.commons.math.util.MathUtils;
23  
24  
25  /**
26   * Solves a linear problem using the Two-Phase Simplex Method.
27   * @version $Revision: 797801 $ $Date: 2009-07-25 13:22:06 -0400 (Sat, 25 Jul 2009) $
28   * @since 2.0
29   */
30  public class SimplexSolver extends AbstractLinearOptimizer {
31  
32      /** Default amount of error to accept in floating point comparisons. */ 
33      private static final double DEFAULT_EPSILON = 1.0e-6;
34  
35      /** Amount of error to accept in floating point comparisons. */ 
36      protected final double epsilon;  
37  
38      /**
39       * Build a simplex solver with default settings.
40       */
41      public SimplexSolver() {
42          this(DEFAULT_EPSILON);
43      }
44  
45      /**
46       * Build a simplex solver with a specified accepted amount of error
47       * @param epsilon the amount of error to accept in floating point comparisons
48       */
49      public SimplexSolver(final double epsilon) {
50          this.epsilon = epsilon;
51      }
52  
53      /**
54       * Returns the column with the most negative coefficient in the objective function row.
55       * @param tableau simple tableau for the problem
56       * @return column with the most negative coefficient
57       */
58      private Integer getPivotColumn(SimplexTableau tableau) {
59          double minValue = 0;
60          Integer minPos = null;
61          for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
62              if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) {
63                  minValue = tableau.getEntry(0, i);
64                  minPos = i;
65              }
66          }
67          return minPos;
68      }
69  
70      /**
71       * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
72       * @param tableau simple tableau for the problem
73       * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
74       * @return row with the minimum ratio
75       */
76      private Integer getPivotRow(final int col, final SimplexTableau tableau) {
77          double minRatio = Double.MAX_VALUE;
78          Integer minRatioPos = null;
79          for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
80              double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
81              if (MathUtils.compareTo(tableau.getEntry(i, col), 0, epsilon) >= 0) {
82                  double ratio = rhs / tableau.getEntry(i, col);
83                  if (ratio < minRatio) {
84                      minRatio = ratio;
85                      minRatioPos = i; 
86                  }
87              }
88          }
89          return minRatioPos;
90      }
91  
92  
93      /**
94       * Runs one iteration of the Simplex method on the given model.
95       * @param tableau simple tableau for the problem
96       * @throws OptimizationException if the maximal iteration count has been
97       * exceeded or if the model is found not to have a bounded solution
98       */
99      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 }