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.math3.optimization.linear;
19  
20  import java.util.ArrayList;
21  import java.util.List;
22  
23  import org.apache.commons.math3.exception.MaxCountExceededException;
24  import org.apache.commons.math3.optimization.PointValuePair;
25  import org.apache.commons.math3.util.Precision;
26  
27  
28  /**
29   * Solves a linear problem using the Two-Phase Simplex Method.
30   *
31   * @deprecated As of 3.1 (to be removed in 4.0).
32   * @since 2.0
33   */
34  @Deprecated
35  public class SimplexSolver extends AbstractLinearOptimizer {
36  
37      /** Default amount of error to accept for algorithm convergence. */
38      private static final double DEFAULT_EPSILON = 1.0e-6;
39  
40      /** Default amount of error to accept in floating point comparisons (as ulps). */
41      private static final int DEFAULT_ULPS = 10;
42  
43      /** Amount of error to accept for algorithm convergence. */
44      private final double epsilon;
45  
46      /** Amount of error to accept in floating point comparisons (as ulps). */
47      private final int maxUlps;
48  
49      /**
50       * Build a simplex solver with default settings.
51       */
52      public SimplexSolver() {
53          this(DEFAULT_EPSILON, DEFAULT_ULPS);
54      }
55  
56      /**
57       * Build a simplex solver with a specified accepted amount of error
58       * @param epsilon the amount of error to accept for algorithm convergence
59       * @param maxUlps amount of error to accept in floating point comparisons
60       */
61      public SimplexSolver(final double epsilon, final int maxUlps) {
62          this.epsilon = epsilon;
63          this.maxUlps = maxUlps;
64      }
65  
66      /**
67       * Returns the column with the most negative coefficient in the objective function row.
68       * @param tableau simple tableau for the problem
69       * @return column with the most negative coefficient
70       */
71      private Integer getPivotColumn(SimplexTableau tableau) {
72          double minValue = 0;
73          Integer minPos = null;
74          for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
75              final double entry = tableau.getEntry(0, i);
76              // check if the entry is strictly smaller than the current minimum
77              // do not use a ulp/epsilon check
78              if (entry < minValue) {
79                  minValue = entry;
80                  minPos = i;
81              }
82          }
83          return minPos;
84      }
85  
86      /**
87       * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
88       * @param tableau simple tableau for the problem
89       * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
90       * @return row with the minimum ratio
91       */
92      private Integer getPivotRow(SimplexTableau tableau, final int col) {
93          // create a list of all the rows that tie for the lowest score in the minimum ratio test
94          List<Integer> minRatioPositions = new ArrayList<Integer>();
95          double minRatio = Double.MAX_VALUE;
96          for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
97              final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
98              final double entry = tableau.getEntry(i, col);
99  
100             if (Precision.compareTo(entry, 0d, maxUlps) > 0) {
101                 final double ratio = rhs / entry;
102                 // check if the entry is strictly equal to the current min ratio
103                 // do not use a ulp/epsilon check
104                 final int cmp = Double.compare(ratio, minRatio);
105                 if (cmp == 0) {
106                     minRatioPositions.add(i);
107                 } else if (cmp < 0) {
108                     minRatio = ratio;
109                     minRatioPositions = new ArrayList<Integer>();
110                     minRatioPositions.add(i);
111                 }
112             }
113         }
114 
115         if (minRatioPositions.size() == 0) {
116             return null;
117         } else if (minRatioPositions.size() > 1) {
118             // there's a degeneracy as indicated by a tie in the minimum ratio test
119 
120             // 1. check if there's an artificial variable that can be forced out of the basis
121             if (tableau.getNumArtificialVariables() > 0) {
122                 for (Integer row : minRatioPositions) {
123                     for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
124                         int column = i + tableau.getArtificialVariableOffset();
125                         final double entry = tableau.getEntry(row, column);
126                         if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) {
127                             return row;
128                         }
129                     }
130                 }
131             }
132 
133             // 2. apply Bland's rule to prevent cycling:
134             //    take the row for which the corresponding basic variable has the smallest index
135             //
136             // see http://www.stanford.edu/class/msande310/blandrule.pdf
137             // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper)
138             //
139             // Additional heuristic: if we did not get a solution after half of maxIterations
140             //                       revert to the simple case of just returning the top-most row
141             // This heuristic is based on empirical data gathered while investigating MATH-828.
142             if (getIterations() < getMaxIterations() / 2) {
143                 Integer minRow = null;
144                 int minIndex = tableau.getWidth();
145                 final int varStart = tableau.getNumObjectiveFunctions();
146                 final int varEnd = tableau.getWidth() - 1;
147                 for (Integer row : minRatioPositions) {
148                     for (int i = varStart; i < varEnd && !row.equals(minRow); i++) {
149                         final Integer basicRow = tableau.getBasicRow(i);
150                         if (basicRow != null && basicRow.equals(row) && i < minIndex) {
151                             minIndex = i;
152                             minRow = row;
153                         }
154                     }
155                 }
156                 return minRow;
157             }
158         }
159         return minRatioPositions.get(0);
160     }
161 
162     /**
163      * Runs one iteration of the Simplex method on the given model.
164      * @param tableau simple tableau for the problem
165      * @throws MaxCountExceededException if the maximal iteration count has been exceeded
166      * @throws UnboundedSolutionException if the model is found not to have a bounded solution
167      */
168     protected void doIteration(final SimplexTableau tableau)
169         throws MaxCountExceededException, UnboundedSolutionException {
170 
171         incrementIterationsCounter();
172 
173         Integer pivotCol = getPivotColumn(tableau);
174         Integer pivotRow = getPivotRow(tableau, pivotCol);
175         if (pivotRow == null) {
176             throw new UnboundedSolutionException();
177         }
178 
179         // set the pivot element to 1
180         double pivotVal = tableau.getEntry(pivotRow, pivotCol);
181         tableau.divideRow(pivotRow, pivotVal);
182 
183         // set the rest of the pivot column to 0
184         for (int i = 0; i < tableau.getHeight(); i++) {
185             if (i != pivotRow) {
186                 final double multiplier = tableau.getEntry(i, pivotCol);
187                 tableau.subtractRow(i, pivotRow, multiplier);
188             }
189         }
190     }
191 
192     /**
193      * Solves Phase 1 of the Simplex method.
194      * @param tableau simple tableau for the problem
195      * @throws MaxCountExceededException if the maximal iteration count has been exceeded
196      * @throws UnboundedSolutionException if the model is found not to have a bounded solution
197      * @throws NoFeasibleSolutionException if there is no feasible solution
198      */
199     protected void solvePhase1(final SimplexTableau tableau)
200         throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException {
201 
202         // make sure we're in Phase 1
203         if (tableau.getNumArtificialVariables() == 0) {
204             return;
205         }
206 
207         while (!tableau.isOptimal()) {
208             doIteration(tableau);
209         }
210 
211         // if W is not zero then we have no feasible solution
212         if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) {
213             throw new NoFeasibleSolutionException();
214         }
215     }
216 
217     /** {@inheritDoc} */
218     @Override
219     public PointValuePair doOptimize()
220         throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException {
221         final SimplexTableau tableau =
222             new SimplexTableau(getFunction(),
223                                getConstraints(),
224                                getGoalType(),
225                                restrictToNonNegative(),
226                                epsilon,
227                                maxUlps);
228 
229         solvePhase1(tableau);
230         tableau.dropPhase1Objective();
231 
232         while (!tableau.isOptimal()) {
233             doIteration(tableau);
234         }
235         return tableau.getSolution();
236     }
237 
238 }