001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.optimization.linear;
019    
020    import org.apache.commons.math.optimization.OptimizationException;
021    import org.apache.commons.math.optimization.RealPointValuePair;
022    import org.apache.commons.math.util.MathUtils;
023    
024    
025    /**
026     * Solves a linear problem using the Two-Phase Simplex Method.
027     * @version $Revision: 797801 $ $Date: 2009-07-25 13:22:06 -0400 (Sat, 25 Jul 2009) $
028     * @since 2.0
029     */
030    public class SimplexSolver extends AbstractLinearOptimizer {
031    
032        /** Default amount of error to accept in floating point comparisons. */ 
033        private static final double DEFAULT_EPSILON = 1.0e-6;
034    
035        /** Amount of error to accept in floating point comparisons. */ 
036        protected final double epsilon;  
037    
038        /**
039         * Build a simplex solver with default settings.
040         */
041        public SimplexSolver() {
042            this(DEFAULT_EPSILON);
043        }
044    
045        /**
046         * Build a simplex solver with a specified accepted amount of error
047         * @param epsilon the amount of error to accept in floating point comparisons
048         */
049        public SimplexSolver(final double epsilon) {
050            this.epsilon = epsilon;
051        }
052    
053        /**
054         * Returns the column with the most negative coefficient in the objective function row.
055         * @param tableau simple tableau for the problem
056         * @return column with the most negative coefficient
057         */
058        private Integer getPivotColumn(SimplexTableau tableau) {
059            double minValue = 0;
060            Integer minPos = null;
061            for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
062                if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) {
063                    minValue = tableau.getEntry(0, i);
064                    minPos = i;
065                }
066            }
067            return minPos;
068        }
069    
070        /**
071         * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
072         * @param tableau simple tableau for the problem
073         * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
074         * @return row with the minimum ratio
075         */
076        private Integer getPivotRow(final int col, final SimplexTableau tableau) {
077            double minRatio = Double.MAX_VALUE;
078            Integer minRatioPos = null;
079            for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
080                double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
081                if (MathUtils.compareTo(tableau.getEntry(i, col), 0, epsilon) >= 0) {
082                    double ratio = rhs / tableau.getEntry(i, col);
083                    if (ratio < minRatio) {
084                        minRatio = ratio;
085                        minRatioPos = i; 
086                    }
087                }
088            }
089            return minRatioPos;
090        }
091    
092    
093        /**
094         * Runs one iteration of the Simplex method on the given model.
095         * @param tableau simple tableau for the problem
096         * @throws OptimizationException if the maximal iteration count has been
097         * exceeded or if the model is found not to have a bounded solution
098         */
099        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    }