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.io.IOException;
21  import java.io.ObjectInputStream;
22  import java.io.ObjectOutputStream;
23  import java.io.Serializable;
24  import java.util.ArrayList;
25  import java.util.Collection;
26  import java.util.HashSet;
27  import java.util.List;
28  import java.util.Set;
29  import java.util.TreeSet;
30  
31  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
32  import org.apache.commons.math3.linear.MatrixUtils;
33  import org.apache.commons.math3.linear.RealMatrix;
34  import org.apache.commons.math3.linear.RealVector;
35  import org.apache.commons.math3.optimization.GoalType;
36  import org.apache.commons.math3.optimization.PointValuePair;
37  import org.apache.commons.math3.util.FastMath;
38  import org.apache.commons.math3.util.Precision;
39  
40  /**
41   * A tableau for use in the Simplex method.
42   *
43   * <p>
44   * Example:
45   * <pre>
46   *   W |  Z |  x1 |  x2 |  x- | s1 |  s2 |  a1 |  RHS
47   * ---------------------------------------------------
48   *  -1    0    0     0     0     0     0     1     0   &lt;= phase 1 objective
49   *   0    1   -15   -10    0     0     0     0     0   &lt;= phase 2 objective
50   *   0    0    1     0     0     1     0     0     2   &lt;= constraint 1
51   *   0    0    0     1     0     0     1     0     3   &lt;= constraint 2
52   *   0    0    1     1     0     0     0     1     4   &lt;= constraint 3
53   * </pre>
54   * W: Phase 1 objective function</br>
55   * Z: Phase 2 objective function</br>
56   * x1 &amp; x2: Decision variables</br>
57   * x-: Extra decision variable to allow for negative values</br>
58   * s1 &amp; s2: Slack/Surplus variables</br>
59   * a1: Artificial variable</br>
60   * RHS: Right hand side</br>
61   * </p>
62   * @version $Id: SimplexTableau.java 1524213 2013-09-17 20:32:50Z tn $
63   * @deprecated As of 3.1 (to be removed in 4.0).
64   * @since 2.0
65   */
66  @Deprecated
67  class SimplexTableau implements Serializable {
68  
69      /** Column label for negative vars. */
70      private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
71  
72      /** Default amount of error to accept in floating point comparisons (as ulps). */
73      private static final int DEFAULT_ULPS = 10;
74  
75      /** The cut-off threshold to zero-out entries. */
76      private static final double CUTOFF_THRESHOLD = 1e-12;
77  
78      /** Serializable version identifier. */
79      private static final long serialVersionUID = -1369660067587938365L;
80  
81      /** Linear objective function. */
82      private final LinearObjectiveFunction f;
83  
84      /** Linear constraints. */
85      private final List<LinearConstraint> constraints;
86  
87      /** Whether to restrict the variables to non-negative values. */
88      private final boolean restrictToNonNegative;
89  
90      /** The variables each column represents */
91      private final List<String> columnLabels = new ArrayList<String>();
92  
93      /** Simple tableau. */
94      private transient RealMatrix tableau;
95  
96      /** Number of decision variables. */
97      private final int numDecisionVariables;
98  
99      /** Number of slack variables. */
100     private final int numSlackVariables;
101 
102     /** Number of artificial variables. */
103     private int numArtificialVariables;
104 
105     /** Amount of error to accept when checking for optimality. */
106     private final double epsilon;
107 
108     /** Amount of error to accept in floating point comparisons. */
109     private final int maxUlps;
110 
111     /**
112      * Build a tableau for a linear problem.
113      * @param f linear objective function
114      * @param constraints linear constraints
115      * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE}
116      * @param restrictToNonNegative whether to restrict the variables to non-negative values
117      * @param epsilon amount of error to accept when checking for optimality
118      */
119     SimplexTableau(final LinearObjectiveFunction f,
120                    final Collection<LinearConstraint> constraints,
121                    final GoalType goalType, final boolean restrictToNonNegative,
122                    final double epsilon) {
123         this(f, constraints, goalType, restrictToNonNegative, epsilon, DEFAULT_ULPS);
124     }
125 
126     /**
127      * Build a tableau for a linear problem.
128      * @param f linear objective function
129      * @param constraints linear constraints
130      * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE}
131      * @param restrictToNonNegative whether to restrict the variables to non-negative values
132      * @param epsilon amount of error to accept when checking for optimality
133      * @param maxUlps amount of error to accept in floating point comparisons
134      */
135     SimplexTableau(final LinearObjectiveFunction f,
136                    final Collection<LinearConstraint> constraints,
137                    final GoalType goalType, final boolean restrictToNonNegative,
138                    final double epsilon,
139                    final int maxUlps) {
140         this.f                      = f;
141         this.constraints            = normalizeConstraints(constraints);
142         this.restrictToNonNegative  = restrictToNonNegative;
143         this.epsilon                = epsilon;
144         this.maxUlps                = maxUlps;
145         this.numDecisionVariables   = f.getCoefficients().getDimension() +
146                                       (restrictToNonNegative ? 0 : 1);
147         this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
148                                       getConstraintTypeCounts(Relationship.GEQ);
149         this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
150                                       getConstraintTypeCounts(Relationship.GEQ);
151         this.tableau = createTableau(goalType == GoalType.MAXIMIZE);
152         initializeColumnLabels();
153     }
154 
155     /**
156      * Initialize the labels for the columns.
157      */
158     protected void initializeColumnLabels() {
159       if (getNumObjectiveFunctions() == 2) {
160         columnLabels.add("W");
161       }
162       columnLabels.add("Z");
163       for (int i = 0; i < getOriginalNumDecisionVariables(); i++) {
164         columnLabels.add("x" + i);
165       }
166       if (!restrictToNonNegative) {
167         columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL);
168       }
169       for (int i = 0; i < getNumSlackVariables(); i++) {
170         columnLabels.add("s" + i);
171       }
172       for (int i = 0; i < getNumArtificialVariables(); i++) {
173         columnLabels.add("a" + i);
174       }
175       columnLabels.add("RHS");
176     }
177 
178     /**
179      * Create the tableau by itself.
180      * @param maximize if true, goal is to maximize the objective function
181      * @return created tableau
182      */
183     protected RealMatrix createTableau(final boolean maximize) {
184 
185         // create a matrix of the correct size
186         int width = numDecisionVariables + numSlackVariables +
187         numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
188         int height = constraints.size() + getNumObjectiveFunctions();
189         Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);
190 
191         // initialize the objective function rows
192         if (getNumObjectiveFunctions() == 2) {
193             matrix.setEntry(0, 0, -1);
194         }
195         int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
196         matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);
197         RealVector objectiveCoefficients =
198             maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
199         copyArray(objectiveCoefficients.toArray(), matrix.getDataRef()[zIndex]);
200         matrix.setEntry(zIndex, width - 1,
201             maximize ? f.getConstantTerm() : -1 * f.getConstantTerm());
202 
203         if (!restrictToNonNegative) {
204             matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
205                 getInvertedCoefficientSum(objectiveCoefficients));
206         }
207 
208         // initialize the constraint rows
209         int slackVar = 0;
210         int artificialVar = 0;
211         for (int i = 0; i < constraints.size(); i++) {
212             LinearConstraint constraint = constraints.get(i);
213             int row = getNumObjectiveFunctions() + i;
214 
215             // decision variable coefficients
216             copyArray(constraint.getCoefficients().toArray(), matrix.getDataRef()[row]);
217 
218             // x-
219             if (!restrictToNonNegative) {
220                 matrix.setEntry(row, getSlackVariableOffset() - 1,
221                     getInvertedCoefficientSum(constraint.getCoefficients()));
222             }
223 
224             // RHS
225             matrix.setEntry(row, width - 1, constraint.getValue());
226 
227             // slack variables
228             if (constraint.getRelationship() == Relationship.LEQ) {
229                 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);  // slack
230             } else if (constraint.getRelationship() == Relationship.GEQ) {
231                 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess
232             }
233 
234             // artificial variables
235             if ((constraint.getRelationship() == Relationship.EQ) ||
236                     (constraint.getRelationship() == Relationship.GEQ)) {
237                 matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
238                 matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
239                 matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
240             }
241         }
242 
243         return matrix;
244     }
245 
246     /**
247      * Get new versions of the constraints which have positive right hand sides.
248      * @param originalConstraints original (not normalized) constraints
249      * @return new versions of the constraints
250      */
251     public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) {
252         List<LinearConstraint> normalized = new ArrayList<LinearConstraint>(originalConstraints.size());
253         for (LinearConstraint constraint : originalConstraints) {
254             normalized.add(normalize(constraint));
255         }
256         return normalized;
257     }
258 
259     /**
260      * Get a new equation equivalent to this one with a positive right hand side.
261      * @param constraint reference constraint
262      * @return new equation
263      */
264     private LinearConstraint normalize(final LinearConstraint constraint) {
265         if (constraint.getValue() < 0) {
266             return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
267                                         constraint.getRelationship().oppositeRelationship(),
268                                         -1 * constraint.getValue());
269         }
270         return new LinearConstraint(constraint.getCoefficients(),
271                                     constraint.getRelationship(), constraint.getValue());
272     }
273 
274     /**
275      * Get the number of objective functions in this tableau.
276      * @return 2 for Phase 1.  1 for Phase 2.
277      */
278     protected final int getNumObjectiveFunctions() {
279         return this.numArtificialVariables > 0 ? 2 : 1;
280     }
281 
282     /**
283      * Get a count of constraints corresponding to a specified relationship.
284      * @param relationship relationship to count
285      * @return number of constraint with the specified relationship
286      */
287     private int getConstraintTypeCounts(final Relationship relationship) {
288         int count = 0;
289         for (final LinearConstraint constraint : constraints) {
290             if (constraint.getRelationship() == relationship) {
291                 ++count;
292             }
293         }
294         return count;
295     }
296 
297     /**
298      * Get the -1 times the sum of all coefficients in the given array.
299      * @param coefficients coefficients to sum
300      * @return the -1 times the sum of all coefficients in the given array.
301      */
302     protected static double getInvertedCoefficientSum(final RealVector coefficients) {
303         double sum = 0;
304         for (double coefficient : coefficients.toArray()) {
305             sum -= coefficient;
306         }
307         return sum;
308     }
309 
310     /**
311      * Checks whether the given column is basic.
312      * @param col index of the column to check
313      * @return the row that the variable is basic in.  null if the column is not basic
314      */
315     protected Integer getBasicRow(final int col) {
316         Integer row = null;
317         for (int i = 0; i < getHeight(); i++) {
318             final double entry = getEntry(i, col);
319             if (Precision.equals(entry, 1d, maxUlps) && (row == null)) {
320                 row = i;
321             } else if (!Precision.equals(entry, 0d, maxUlps)) {
322                 return null;
323             }
324         }
325         return row;
326     }
327 
328     /**
329      * Removes the phase 1 objective function, positive cost non-artificial variables,
330      * and the non-basic artificial variables from this tableau.
331      */
332     protected void dropPhase1Objective() {
333         if (getNumObjectiveFunctions() == 1) {
334             return;
335         }
336 
337         Set<Integer> columnsToDrop = new TreeSet<Integer>();
338         columnsToDrop.add(0);
339 
340         // positive cost non-artificial variables
341         for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) {
342             final double entry = tableau.getEntry(0, i);
343             if (Precision.compareTo(entry, 0d, epsilon) > 0) {
344                 columnsToDrop.add(i);
345             }
346         }
347 
348         // non-basic artificial variables
349         for (int i = 0; i < getNumArtificialVariables(); i++) {
350             int col = i + getArtificialVariableOffset();
351             if (getBasicRow(col) == null) {
352                 columnsToDrop.add(col);
353             }
354         }
355 
356         double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()];
357         for (int i = 1; i < getHeight(); i++) {
358             int col = 0;
359             for (int j = 0; j < getWidth(); j++) {
360                 if (!columnsToDrop.contains(j)) {
361                     matrix[i - 1][col++] = tableau.getEntry(i, j);
362                 }
363             }
364         }
365 
366         // remove the columns in reverse order so the indices are correct
367         Integer[] drop = columnsToDrop.toArray(new Integer[columnsToDrop.size()]);
368         for (int i = drop.length - 1; i >= 0; i--) {
369             columnLabels.remove((int) drop[i]);
370         }
371 
372         this.tableau = new Array2DRowRealMatrix(matrix);
373         this.numArtificialVariables = 0;
374     }
375 
376     /**
377      * @param src the source array
378      * @param dest the destination array
379      */
380     private void copyArray(final double[] src, final double[] dest) {
381         System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length);
382     }
383 
384     /**
385      * Returns whether the problem is at an optimal state.
386      * @return whether the model has been solved
387      */
388     boolean isOptimal() {
389         for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) {
390             final double entry = tableau.getEntry(0, i);
391             if (Precision.compareTo(entry, 0d, epsilon) < 0) {
392                 return false;
393             }
394         }
395         return true;
396     }
397 
398     /**
399      * Get the current solution.
400      * @return current solution
401      */
402     protected PointValuePair getSolution() {
403       int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL);
404       Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null;
405       double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset());
406 
407       Set<Integer> basicRows = new HashSet<Integer>();
408       double[] coefficients = new double[getOriginalNumDecisionVariables()];
409       for (int i = 0; i < coefficients.length; i++) {
410           int colIndex = columnLabels.indexOf("x" + i);
411           if (colIndex < 0) {
412             coefficients[i] = 0;
413             continue;
414           }
415           Integer basicRow = getBasicRow(colIndex);
416           if (basicRow != null && basicRow == 0) {
417               // if the basic row is found to be the objective function row
418               // set the coefficient to 0 -> this case handles unconstrained
419               // variables that are still part of the objective function
420               coefficients[i] = 0;
421           } else if (basicRows.contains(basicRow)) {
422               // if multiple variables can take a given value
423               // then we choose the first and set the rest equal to 0
424               coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative);
425           } else {
426               basicRows.add(basicRow);
427               coefficients[i] =
428                   (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
429                   (restrictToNonNegative ? 0 : mostNegative);
430           }
431       }
432       return new PointValuePair(coefficients, f.getValue(coefficients));
433     }
434 
435     /**
436      * Subtracts a multiple of one row from another.
437      * <p>
438      * After application of this operation, the following will hold:
439      * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre>
440      *
441      * @param dividendRow index of the row
442      * @param divisor value of the divisor
443      */
444     protected void divideRow(final int dividendRow, final double divisor) {
445         for (int j = 0; j < getWidth(); j++) {
446             tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor);
447         }
448     }
449 
450     /**
451      * Subtracts a multiple of one row from another.
452      * <p>
453      * After application of this operation, the following will hold:
454      * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre>
455      *
456      * @param minuendRow row index
457      * @param subtrahendRow row index
458      * @param multiple multiplication factor
459      */
460     protected void subtractRow(final int minuendRow, final int subtrahendRow,
461                                final double multiple) {
462         for (int i = 0; i < getWidth(); i++) {
463             double result = tableau.getEntry(minuendRow, i) - tableau.getEntry(subtrahendRow, i) * multiple;
464             // cut-off values smaller than the CUTOFF_THRESHOLD, otherwise may lead to numerical instabilities
465             if (FastMath.abs(result) < CUTOFF_THRESHOLD) {
466                 result = 0.0;
467             }
468             tableau.setEntry(minuendRow, i, result);
469         }
470     }
471 
472     /**
473      * Get the width of the tableau.
474      * @return width of the tableau
475      */
476     protected final int getWidth() {
477         return tableau.getColumnDimension();
478     }
479 
480     /**
481      * Get the height of the tableau.
482      * @return height of the tableau
483      */
484     protected final int getHeight() {
485         return tableau.getRowDimension();
486     }
487 
488     /**
489      * Get an entry of the tableau.
490      * @param row row index
491      * @param column column index
492      * @return entry at (row, column)
493      */
494     protected final double getEntry(final int row, final int column) {
495         return tableau.getEntry(row, column);
496     }
497 
498     /**
499      * Set an entry of the tableau.
500      * @param row row index
501      * @param column column index
502      * @param value for the entry
503      */
504     protected final void setEntry(final int row, final int column,
505                                   final double value) {
506         tableau.setEntry(row, column, value);
507     }
508 
509     /**
510      * Get the offset of the first slack variable.
511      * @return offset of the first slack variable
512      */
513     protected final int getSlackVariableOffset() {
514         return getNumObjectiveFunctions() + numDecisionVariables;
515     }
516 
517     /**
518      * Get the offset of the first artificial variable.
519      * @return offset of the first artificial variable
520      */
521     protected final int getArtificialVariableOffset() {
522         return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
523     }
524 
525     /**
526      * Get the offset of the right hand side.
527      * @return offset of the right hand side
528      */
529     protected final int getRhsOffset() {
530         return getWidth() - 1;
531     }
532 
533     /**
534      * Get the number of decision variables.
535      * <p>
536      * If variables are not restricted to positive values, this will include 1 extra decision variable to represent
537      * the absolute value of the most negative variable.
538      *
539      * @return number of decision variables
540      * @see #getOriginalNumDecisionVariables()
541      */
542     protected final int getNumDecisionVariables() {
543         return numDecisionVariables;
544     }
545 
546     /**
547      * Get the original number of decision variables.
548      * @return original number of decision variables
549      * @see #getNumDecisionVariables()
550      */
551     protected final int getOriginalNumDecisionVariables() {
552         return f.getCoefficients().getDimension();
553     }
554 
555     /**
556      * Get the number of slack variables.
557      * @return number of slack variables
558      */
559     protected final int getNumSlackVariables() {
560         return numSlackVariables;
561     }
562 
563     /**
564      * Get the number of artificial variables.
565      * @return number of artificial variables
566      */
567     protected final int getNumArtificialVariables() {
568         return numArtificialVariables;
569     }
570 
571     /**
572      * Get the tableau data.
573      * @return tableau data
574      */
575     protected final double[][] getData() {
576         return tableau.getData();
577     }
578 
579     @Override
580     public boolean equals(Object other) {
581 
582       if (this == other) {
583         return true;
584       }
585 
586       if (other instanceof SimplexTableau) {
587           SimplexTableau rhs = (SimplexTableau) other;
588           return (restrictToNonNegative  == rhs.restrictToNonNegative) &&
589                  (numDecisionVariables   == rhs.numDecisionVariables) &&
590                  (numSlackVariables      == rhs.numSlackVariables) &&
591                  (numArtificialVariables == rhs.numArtificialVariables) &&
592                  (epsilon                == rhs.epsilon) &&
593                  (maxUlps                == rhs.maxUlps) &&
594                  f.equals(rhs.f) &&
595                  constraints.equals(rhs.constraints) &&
596                  tableau.equals(rhs.tableau);
597       }
598       return false;
599     }
600 
601     @Override
602     public int hashCode() {
603         return Boolean.valueOf(restrictToNonNegative).hashCode() ^
604                numDecisionVariables ^
605                numSlackVariables ^
606                numArtificialVariables ^
607                Double.valueOf(epsilon).hashCode() ^
608                maxUlps ^
609                f.hashCode() ^
610                constraints.hashCode() ^
611                tableau.hashCode();
612     }
613 
614     /**
615      * Serialize the instance.
616      * @param oos stream where object should be written
617      * @throws IOException if object cannot be written to stream
618      */
619     private void writeObject(ObjectOutputStream oos)
620         throws IOException {
621         oos.defaultWriteObject();
622         MatrixUtils.serializeRealMatrix(tableau, oos);
623     }
624 
625     /**
626      * Deserialize the instance.
627      * @param ois stream from which the object should be read
628      * @throws ClassNotFoundException if a class in the stream cannot be found
629      * @throws IOException if object cannot be read from the stream
630      */
631     private void readObject(ObjectInputStream ois)
632       throws ClassNotFoundException, IOException {
633         ois.defaultReadObject();
634         MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
635     }
636 }