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