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