added an implementation of Dantzig's simplex algorithm

to solve constrained linear optimization problems
JIRA: MATH-246

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@758920 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-03-26 23:25:30 +00:00
parent b2db570b31
commit 61e775706d
19 changed files with 1766 additions and 3 deletions

View File

@ -4,6 +4,9 @@ Copyright 2001-2008 The Apache Software Foundation
This product includes software developed by
The Apache Software Foundation (http://www.apache.org/).
This product includes software developed by
Benjamin McCann (http://www.benmccann.com) and distributed with with
the following copyright: Copyright 2009 Google Inc.
This product includes software translated from the lmder, lmpar
and qrsolv Fortran routines from the Minpack package and

View File

@ -143,6 +143,7 @@
</contributor>
<contributor>
<name>Benjamin McCann</name>
<url>http://www.benmccann.com"</url>
</contributor>
<contributor>
<name>Fredrik Norin</name>

View File

@ -146,6 +146,14 @@ public class MessagesResources_fr
{ "unable to bracket optimum in line search",
"impossible d''encadrer l''optimum lors de la recherche lin\u00e9aire" },
// org.apache.commons.math.optimization.linear2.NoFeasibleSolutionException
{ "no feasible solution",
"aucune solution r\u00e9alisable" },
// org.apache.commons.math.optimization.linear2.UnboundedSolutionException
{ "unbounded solution",
"solution non born\u00e9e" },
// org.apache.commons.math.geometry.CardanEulerSingularityException
{ "Cardan angles singularity",
"singularit\u00e9 d''angles de Cardan" },

View File

@ -20,7 +20,7 @@ package org.apache.commons.math.optimization;
import org.apache.commons.math.ConvergenceException;
/**
* This class represents exceptions thrown by the estimation solvers.
* This class represents exceptions thrown by optimizers.
*
* @version $Revision$ $Date$
* @since 1.2

View File

@ -53,7 +53,7 @@ public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
/** Simple constructor with default settings.
* <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
* and the maximal number of evaluation is set to
* {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}.
* {@link AbstractLinearOptimizer#DEFAULT_MAX_ITERATIONS}.
* @param useLU if true, the normal equations will be solved using LU
* decomposition, otherwise they will be solved using QR decomposition
*/

View File

@ -63,7 +63,7 @@ public class NonLinearConjugateGradientOptimizer
/** Simple constructor with default settings.
* <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
* and the maximal number of evaluation is set to
* {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_EVALUATIONS}.
* {@link AbstractLinearOptimizer#DEFAULT_MAX_EVALUATIONS}.
* @param updateFormula formula to use for updating the &beta; parameter,
* must be one of {@link UpdateFormula#FLETCHER_REEVES} or {@link
* UpdateFormula#POLAK_RIBIERE}

View File

@ -0,0 +1,124 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.util.Collection;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealPointValuePair;
/**
* Base class for implementing linear optimizers.
* <p>This base class handles the boilerplate methods associated to thresholds
* settings and iterations counters.</p>
* @version $Revision$ $Date$
* @since 2.0
*
*/
public abstract class AbstractLinearOptimizer implements LinearOptimizer {
/** Serializable version identifier */
private static final long serialVersionUID = 8581325080951819398L;
/** Default maximal number of iterations allowed. */
public static final int DEFAULT_MAX_ITERATIONS = 100;
/** Maximal number of iterations allowed. */
private int maxIterations;
/** Number of iterations already performed. */
private int iterations;
/** Linear objective function. */
protected LinearObjectiveFunction f;
/** Linear constraints. */
protected Collection<LinearConstraint> constraints;
/** Type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE}. */
protected GoalType goalType;
/** Whether to restrict the variables to non-negative values. */
protected boolean restrictToNonNegative;
/** Simple constructor with default settings.
* <p>The maximal number of evaluation is set to its default value.</p>
*/
protected AbstractLinearOptimizer() {
setMaxIterations(DEFAULT_MAX_ITERATIONS);
}
/** {@inheritDoc} */
public void setMaxIterations(int maxIterations) {
this.maxIterations = maxIterations;
}
/** {@inheritDoc} */
public int getMaxIterations() {
return maxIterations;
}
/** {@inheritDoc} */
public int getIterations() {
return iterations;
}
/** Increment the iterations counter by 1.
* @exception OptimizationException if the maximal number
* of iterations is exceeded
*/
protected void incrementIterationsCounter()
throws OptimizationException {
if (++iterations > maxIterations) {
if (++iterations > maxIterations) {
throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
}
}
}
/** {@inheritDoc} */
public RealPointValuePair optimize(final LinearObjectiveFunction f,
final Collection<LinearConstraint> constraints,
final GoalType goalType, final boolean restrictToNonNegative)
throws OptimizationException {
// store linear problem characteristics
this.f = f;
this.constraints = constraints;
this.goalType = goalType;
this.restrictToNonNegative = restrictToNonNegative;
iterations = 0;
// solve the problem
return doOptimize();
}
/** Perform the bulk of optimization algorithm.
* @return the point/value pair giving the optimal value for objective function
* @exception OptimizationException if no solution fulfilling the constraints
* can be found in the allowed number of iterations
*/
abstract protected RealPointValuePair doOptimize()
throws OptimizationException;
}

View File

@ -0,0 +1,183 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.io.Serializable;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.RealVectorImpl;
/**
* A linear constraint for a linear optimization problem.
* <p>
* A linear constraint has one of the forms:
* <ul>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> = v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> &lt;= v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> >= v</li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> =
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> &lt;=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> >=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* </ul>
* The c<sub>i</sub>, l<sub>i</sub> or r<sub>i</sub> are the coefficients of the constraints, the x<sub>i</sub>
* are the coordinates of the current point and v is the value of the constraint.
* </p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class LinearConstraint implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = -764632794033034092L;
/** Coefficients of the constraint (left hand side). */
private final RealVector coefficients;
/** Relationship between left and right hand sides (=, &lt;=, >=). */
private final Relationship relationship;
/** Value of the constraint (right hand side). */
private final double value;
/**
* Build a constraint involving a single linear equation.
* <p>
* A linear constraint with a single linear equation has one of the forms:
* <ul>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> = v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> &lt;= v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> >= v</li>
* </ul>
* </p>
* @param coefficients The coefficients of the constraint (left hand side)
* @param relationship The type of (in)equality used in the constraint
* @param value The value of the constraint (right hand side)
*/
public LinearConstraint(final double[] coefficients, final Relationship relationship,
final double value) {
this(new RealVectorImpl(coefficients), relationship, value);
}
/**
* Build a constraint involving a single linear equation.
* <p>
* A linear constraint with a single linear equation has one of the forms:
* <ul>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> = v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> &lt;= v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> >= v</li>
* </ul>
* </p>
* @param coefficients The coefficients of the constraint (left hand side)
* @param relationship The type of (in)equality used in the constraint
* @param value The value of the constraint (right hand side)
*/
public LinearConstraint(final RealVector coefficients, final Relationship relationship,
final double value) {
this.coefficients = coefficients;
this.relationship = relationship;
this.value = value;
}
/**
* Build a constraint involving two linear equations.
* <p>
* A linear constraint with two linear equation has one of the forms:
* <ul>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> =
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> &lt;=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> >=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* </ul>
* </p>
* @param lhsCoefficients The coefficients of the linear expression on the left hand side of the constraint
* @param lhsConstant The constant term of the linear expression on the left hand side of the constraint
* @param relationship The type of (in)equality used in the constraint
* @param rhsCoefficients The coefficients of the linear expression on the right hand side of the constraint
* @param rhsConstant The constant term of the linear expression on the right hand side of the constraint
*/
public LinearConstraint(final double[] lhsCoefficients, final double lhsConstant,
final Relationship relationship,
final double[] rhsCoefficients, final double rhsConstant) {
double[] sub = new double[lhsCoefficients.length];
for (int i = 0; i < sub.length; ++i) {
sub[i] = lhsCoefficients[i] - rhsCoefficients[i];
}
this.coefficients = new RealVectorImpl(sub, false);
this.relationship = relationship;
this.value = rhsConstant - lhsConstant;
}
/**
* Build a constraint involving two linear equations.
* <p>
* A linear constraint with two linear equation has one of the forms:
* <ul>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> =
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> &lt;=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> >=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* </ul>
* </p>
* @param lhsCoefficients The coefficients of the linear expression on the left hand side of the constraint
* @param lhsConstant The constant term of the linear expression on the left hand side of the constraint
* @param relationship The type of (in)equality used in the constraint
* @param rhsCoefficients The coefficients of the linear expression on the right hand side of the constraint
* @param rhsConstant The constant term of the linear expression on the right hand side of the constraint
*/
public LinearConstraint(final RealVector lhsCoefficients, final double lhsConstant,
final Relationship relationship,
final RealVector rhsCoefficients, final double rhsConstant) {
this.coefficients = lhsCoefficients.subtract(rhsCoefficients);
this.relationship = relationship;
this.value = rhsConstant - lhsConstant;
}
/**
* Get the coefficients of the constraint (left hand side).
* @return coefficients of the constraint (left hand side)
*/
public RealVector getCoefficients() {
return coefficients;
}
/**
* Get the relationship between left and right hand sides.
* @return relationship between left and right hand sides
*/
public Relationship getRelationship() {
return relationship;
}
/**
* Get the value of the constraint (right hand side).
* @return value of the constraint (right hand side)
*/
public double getValue() {
return value;
}
}

View File

@ -0,0 +1,100 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.io.Serializable;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.RealVectorImpl;
/**
* An objective function for a linear optimization problem.
* <p>
* A linear objective function has one the form:
* <pre>
* c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> + d
* </pre>
* The c<sub>i</sub> and d are the coefficients of the equation,
* the x<sub>i</sub> are the coordinates of the current point.
* </p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class LinearObjectiveFunction implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = -4531815507568396090L;
/** Coefficients of the constraint (c<sub>i</sub>). */
private final RealVector coefficients;
/** Constant term of the linear equation. */
private final double constantTerm;
/**
* @param coefficients The coefficients for the linear equation being optimized
* @param constantTerm The constant term of the linear equation
*/
public LinearObjectiveFunction(double[] coefficients, double constantTerm) {
this(new RealVectorImpl(coefficients), constantTerm);
}
/**
* @param coefficients The coefficients for the linear equation being optimized
* @param constantTerm The constant term of the linear equation
*/
public LinearObjectiveFunction(RealVector coefficients, double constantTerm) {
this.coefficients = coefficients;
this.constantTerm = constantTerm;
}
/**
* Get the coefficients of the linear equation being optimized.
* @return coefficients of the linear equation being optimized
*/
public RealVector getCoefficients() {
return coefficients;
}
/**
* Get the constant of the linear equation being optimized.
* @return constant of the linear equation being optimized
*/
public double getConstantTerm() {
return constantTerm;
}
/**
* Compute the value of the linear equation at the current point
* @param point point at which linear equation must be evaluated
* @return value of the linear equation at the current point
*/
public double getValue(final double[] point) {
return coefficients.dotProduct(point) + constantTerm;
}
/**
* Compute the value of the linear equation at the current point
* @param point point at which linear equation must be evaluated
* @return value of the linear equation at the current point
*/
public double getValue(final RealVector point) {
return coefficients.dotProduct(point) + constantTerm;
}
}

View File

@ -0,0 +1,91 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.io.Serializable;
import java.util.Collection;
import org.apache.commons.math.analysis.DifferentiableMultivariateRealFunction;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealPointValuePair;
/**
* This interface represents an optimization algorithm for linear problems.
* <p>Optimization algorithms find the input point set that either {@link GoalType
* maximize or minimize} an objective function. In the linear case the form of
* the function is restricted to
* <pre>
* c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> = v
* </pre>
* and there may be linear constraints too, of one of the forms:
* <ul>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> = v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> &lt;= v</li>
* <li>c<sub>1</sub>x<sub>1</sub> + ... c<sub>n</sub>x<sub>n</sub> >= v</li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> =
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> &lt;=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* <li>l<sub>1</sub>x<sub>1</sub> + ... l<sub>n</sub>x<sub>n</sub> + l<sub>cst</sub> >=
* r<sub>1</sub>x<sub>1</sub> + ... r<sub>n</sub>x<sub>n</sub> + r<sub>cst</sub></li>
* </ul>
* where the c<sub>i</sub>, l<sub>i</sub> or r<sub>i</sub> are the coefficients of
* the constraints, the x<sub>i</sub> are the coordinates of the current point and
* v is the value of the constraint.
* </p>
* @version $Revision$ $Date$
* @since 2.0
*/
public interface LinearOptimizer extends Serializable {
/** Set the maximal number of iterations of the algorithm.
* @param maxIterations maximal number of function calls
*/
void setMaxIterations(int maxIterations);
/** Get the maximal number of iterations of the algorithm.
* @return maximal number of iterations
*/
int getMaxIterations();
/** Get the number of iterations realized by the algorithm.
* <p>
* The number of evaluations corresponds to the last call to the
* {@link #optimize(DifferentiableMultivariateRealFunction, GoalType, double[]) optimize}
* method. It is 0 if the method has not been called yet.
* </p>
* @return number of iterations
*/
int getIterations();
/** Optimizes an objective function.
* @param f linear objective function
* @param constraints linear constraints
* @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
* or {@link GoalType#MINIMIZE}
* @param restrictToNonNegative whether to restrict the variables to non-negative values
* @return point/value pair giving the optimal value for objective function
* @exception OptimizationException if no solution fulfilling the constraints
* can be found in the allowed number of iterations
*/
RealPointValuePair optimize(LinearObjectiveFunction f, Collection<LinearConstraint> constraints,
GoalType goalType, boolean restrictToNonNegative)
throws OptimizationException;
}

View File

@ -0,0 +1,40 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import org.apache.commons.math.optimization.OptimizationException;
/**
* This class represents exceptions thrown by optimizers when no solution
* fulfills the constraints.
* @version $Revision$ $Date$
* @since 2.0
*/
public class NoFeasibleSolutionException extends OptimizationException {
/** Serializable version identifier. */
private static final long serialVersionUID = -3044253632189082760L;
/**
* Simple constructor using a default message.
*/
public NoFeasibleSolutionException() {
super("no feasible solution");
}
}

View File

@ -0,0 +1,67 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
/**
* Types of relationships between two cells in a Solver {@link LinearConstraint}.
* @version $Revision$ $Date$
* @since 2.0
*/
public enum Relationship {
/** Equality relationship. */
EQ("="),
/** Lesser than or equal relationship. */
LEQ("<="),
/** Greater than or equal relationship. */
GEQ(">=");
/** Display string for the relationship. */
private String stringValue;
/** Simple constructor.
* @param stringValue display string for the relationship
*/
private Relationship(String stringValue) {
this.stringValue = stringValue;
}
/** {@inheritDoc} */
@Override
public String toString() {
return stringValue;
}
/**
* Get the relationship obtained when multiplying all coefficients by -1.
* @return relationship obtained when multiplying all coefficients by -1
*/
public Relationship oppositeRelationship() {
switch (this) {
case LEQ :
return GEQ;
case GEQ :
return LEQ;
default :
return EQ;
}
}
}

View File

@ -0,0 +1,197 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealPointValuePair;
import org.apache.commons.math.util.MathUtils;
/**
* Solves a linear problem using the Two-Phase Simplex Method.
* @version $Revision$ $Date$
* @since 2.0
*/
public class SimplexSolver extends AbstractLinearOptimizer {
/** Serializable version identifier. */
private static final long serialVersionUID = -4886937648715323786L;
/** Default amount of error to accept in floating point comparisons. */
private static final double DEFAULT_EPSILON = 1.0e-10;
/** Amount of error to accept in floating point comparisons. */
protected final double epsilon;
/**
* Build a simplex solver with default settings.
*/
public SimplexSolver() {
this(DEFAULT_EPSILON);
}
/**
* Build a simplex solver with a specified accepted amount of error
* @param epsilon the amount of error to accept in floating point comparisons
*/
public SimplexSolver(final double epsilon) {
this.epsilon = epsilon;
}
/**
* Returns the column with the most negative coefficient in the objective function row.
* @param tableau simple tableau for the problem
* @return column with the most negative coefficient
*/
private Integer getPivotColumn(SimplexTableau tableau) {
double minValue = 0;
Integer minPos = null;
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
if (tableau.getEntry(0, i) < minValue) {
minValue = tableau.getEntry(0, i);
minPos = i;
}
}
return minPos;
}
/**
* Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
* @param tableau simple tableau for the problem
* @param col the column to test the ratio of. See {@link #getPivotColumn()}
* @return row with the minimum ratio
*/
private Integer getPivotRow(final int col, final SimplexTableau tableau) {
double minRatio = Double.MAX_VALUE;
Integer minRatioPos = null;
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
if (tableau.getEntry(i, col) >= 0) {
double ratio = rhs / tableau.getEntry(i, col);
if (ratio < minRatio) {
minRatio = ratio;
minRatioPos = i;
}
}
}
return minRatioPos;
}
/**
* Runs one iteration of the Simplex method on the given model.
* @param tableau simple tableau for the problem
* @throws OptimizationException if the maximal iteration count has been
* exceeded or if the model is found not to have a bounded solution
*/
protected void doIteration(final SimplexTableau tableau)
throws OptimizationException {
incrementIterationsCounter();
Integer pivotCol = getPivotColumn(tableau);
Integer pivotRow = getPivotRow(pivotCol, tableau);
if (pivotRow == null) {
throw new UnboundedSolutionException();
}
// set the pivot element to 1
double pivotVal = tableau.getEntry(pivotRow, pivotCol);
tableau.divideRow(pivotRow, pivotVal);
// set the rest of the pivot column to 0
for (int i = 0; i < tableau.getHeight(); i++) {
if (i != pivotRow) {
double multiplier = tableau.getEntry(i, pivotCol);
tableau.subtractRow(i, pivotRow, multiplier);
}
}
}
/**
* Checks whether Phase 1 is solved.
* @param tableau simple tableau for the problem
* @return whether Phase 1 is solved
*/
private boolean isPhase1Solved(final SimplexTableau tableau) {
if (tableau.getNumArtificialVariables() == 0) {
return true;
}
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
if (tableau.getEntry(0, i) < 0) {
return false;
}
}
return true;
}
/**
* Returns whether the problem is at an optimal state.
* @param tableau simple tableau for the problem
* @return whether the model has been solved
*/
public boolean isOptimal(final SimplexTableau tableau) {
if (tableau.getNumArtificialVariables() > 0) {
return false;
}
for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
if (tableau.getEntry(0, i) < 0) {
return false;
}
}
return true;
}
/**
* Solves Phase 1 of the Simplex method.
* @param tableau simple tableau for the problem
* @exception OptimizationException if the maximal number of iterations is
* exceeded, or if the problem is found not to have a bounded solution, or
* if there is no feasible solution
*/
protected void solvePhase1(final SimplexTableau tableau)
throws OptimizationException {
// make sure we're in Phase 1
if (tableau.getNumArtificialVariables() == 0) {
return;
}
while (!isPhase1Solved(tableau)) {
doIteration(tableau);
}
// if W is not zero then we have no feasible solution
if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0, epsilon)) {
throw new NoFeasibleSolutionException();
}
}
/** {@inheritDoc} */
public RealPointValuePair doOptimize()
throws OptimizationException {
final SimplexTableau tableau =
new SimplexTableau(f, constraints, goalType, restrictToNonNegative);
solvePhase1(tableau);
tableau.discardArtificialVariables();
while (!isOptimal(tableau)) {
doIteration(tableau);
}
return tableau.getSolution();
}
}

View File

@ -0,0 +1,473 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixImpl;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.RealPointValuePair;
/**
* A tableau for use in the Simplex method.
*
* <p>
* Example:
* <pre>
* W | Z | x1 | x2 | x- | s1 | s2 | a1 | RHS
* ---------------------------------------------------
* -1 0 0 0 0 0 0 1 0 &lt;= phase 1 objective
* 0 1 -15 -10 0 0 0 0 0 &lt;= phase 2 objective
* 0 0 1 0 0 1 0 0 2 &lt;= constraint 1
* 0 0 0 1 0 0 1 0 3 &lt;= constraint 2
* 0 0 1 1 0 0 0 1 4 &lt;= constraint 3
* </pre>
* W: Phase 1 objective function</br>
* Z: Phase 2 objective function</br>
* x1 &amp; x2: Decision variables</br>
* x-: Extra decision variable to allow for negative values</br>
* s1 &amp; s2: Slack/Surplus variables</br>
* a1: Artificial variable</br>
* RHS: Right hand side</br>
* </p>
* @version $Revision$ $Date$
* @since 2.0
*/
class SimplexTableau implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = -1369660067587938365L;
/** Linear objective function. */
private final LinearObjectiveFunction f;
/** Linear constraints. */
private final Collection<LinearConstraint> constraints;
/** Whether to restrict the variables to non-negative values. */
private final boolean restrictToNonNegative;
/** Simple tableau. */
protected RealMatrix tableau;
/** Number of decision variables. */
protected final int numDecisionVariables;
/** Number of slack variables. */
protected final int numSlackVariables;
/** Number of artificial variables. */
protected int numArtificialVariables;
/**
* Build a tableau for a linear problem.
* @param f linear objective function
* @param constraints linear constraints
* @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
* or {@link GoalType#MINIMIZE}
* @param restrictToNonNegative whether to restrict the variables to non-negative values
*/
SimplexTableau(final LinearObjectiveFunction f,
final Collection<LinearConstraint> constraints,
final GoalType goalType, final boolean restrictToNonNegative) {
this.f = f;
this.constraints = constraints;
this.restrictToNonNegative = restrictToNonNegative;
this.numDecisionVariables = getNumVariables() + (restrictToNonNegative ? 0 : 1);
this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) +
getConstraintTypeCounts(Relationship.GEQ);
this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
getConstraintTypeCounts(Relationship.GEQ);
this.tableau = new RealMatrixImpl(createTableau(goalType == GoalType.MAXIMIZE));
initialize();
}
/**
* Create the tableau by itself.
* @param maximize if true, goal is to maximize the objective function
* @return created tableau
*/
protected double[][] createTableau(final boolean maximize) {
// create a matrix of the correct size
List<LinearConstraint> constraints = getNormalizedConstraints();
int width = numDecisionVariables + numSlackVariables +
numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
int height = constraints.size() + getNumObjectiveFunctions();
double[][] matrix = new double[height][width];
// initialize the objective function rows
if (getNumObjectiveFunctions() == 2) {
matrix[0][0] = -1;
}
int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
matrix[zIndex][zIndex] = maximize ? 1 : -1;
RealVector objectiveCoefficients =
maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
copyArray(objectiveCoefficients.getData(), matrix[zIndex], getNumObjectiveFunctions());
matrix[zIndex][width - 1] =
maximize ? f.getConstantTerm() : -1 * f.getConstantTerm();
if (!restrictToNonNegative) {
matrix[zIndex][getSlackVariableOffset() - 1] =
getInvertedCoeffiecientSum(objectiveCoefficients);
}
// initialize the constraint rows
int slackVar = 0;
int artificialVar = 0;
for (int i = 0; i < constraints.size(); i++) {
LinearConstraint constraint = constraints.get(i);
int row = getNumObjectiveFunctions() + i;
// decision variable coefficients
copyArray(constraint.getCoefficients().getData(), matrix[row], 1);
// x-
if (!restrictToNonNegative) {
matrix[row][getSlackVariableOffset() - 1] =
getInvertedCoeffiecientSum(constraint.getCoefficients());
}
// RHS
matrix[row][width - 1] = constraint.getValue();
// slack variables
if (constraint.getRelationship() == Relationship.LEQ) {
matrix[row][getSlackVariableOffset() + slackVar++] = 1; // slack
} else if (constraint.getRelationship() == Relationship.GEQ) {
matrix[row][getSlackVariableOffset() + slackVar++] = -1; // excess
}
// artificial variables
if ((constraint.getRelationship() == Relationship.EQ) ||
(constraint.getRelationship() == Relationship.GEQ)) {
matrix[0][getArtificialVariableOffset() + artificialVar] = 1;
matrix[row][getArtificialVariableOffset() + artificialVar++] = 1;
}
}
return matrix;
}
/** Get the number of variables.
* @return number of variables
*/
public int getNumVariables() {
return f.getCoefficients().getDimension();
}
/**
* Get new versions of the constraints which have positive right hand sides.
* @return new versions of the constraints
*/
public List<LinearConstraint> getNormalizedConstraints() {
List<LinearConstraint> normalized = new ArrayList<LinearConstraint>();
for (LinearConstraint constraint : constraints) {
normalized.add(normalize(constraint));
}
return normalized;
}
/**
* Get a new equation equivalent to this one with a positive right hand side.
* @param constraint reference constraint
* @return new equation
*/
private LinearConstraint normalize(final LinearConstraint constraint) {
if (constraint.getValue() < 0) {
return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
constraint.getRelationship().oppositeRelationship(),
-1 * constraint.getValue());
}
return new LinearConstraint(constraint.getCoefficients(),
constraint.getRelationship(), constraint.getValue());
}
/**
* Get the number of objective functions in this tableau.
* @return 2 for Phase 1. 1 for Phase 2.
*/
protected final int getNumObjectiveFunctions() {
return this.numArtificialVariables > 0 ? 2 : 1;
}
/**
* Get a count of constraints corresponding to a specified relationship.
* @param relationship relationship to count
* @return number of constraint with the specified relationship
*/
private int getConstraintTypeCounts(final Relationship relationship) {
int count = 0;
for (final LinearConstraint constraint : constraints) {
if (constraint.getRelationship() == relationship) {
++count;
}
}
return count;
}
/**
* Puts the tableau in proper form by zeroing out the artificial variables
* in the objective function via elementary row operations.
*/
private void initialize() {
for (int artificialVar = 0; artificialVar < numArtificialVariables; artificialVar++) {
int row = getBasicRow(getArtificialVariableOffset() + artificialVar);
subtractRow(0, row, 1.0);
}
}
/**
* Get the -1 times the sum of all coefficients in the given array.
* @param coefficients coefficients to sum
* @return the -1 times the sum of all coefficients in the given array.
*/
protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
double sum = 0;
for (double coefficient : coefficients.getData()) {
sum -= coefficient;
}
return sum;
}
/**
* Checks whether the given column is basic.
* @param col index of the column to check
* @return the row that the variable is basic in. null if the column is not basic
*/
private Integer getBasicRow(final int col) {
Integer row = null;
for (int i = getNumObjectiveFunctions(); i < getHeight(); i++) {
if (getEntry(i, col) != 0.0) {
if (row == null) {
row = i;
} else {
return null;
}
}
}
return row;
}
/**
* Removes the phase 1 objective function and artificial variables from this tableau.
*/
protected void discardArtificialVariables() {
if (numArtificialVariables == 0) {
return;
}
int width = getWidth() - numArtificialVariables - 1;
int height = getHeight() - 1;
double[][] matrix = new double[height][width];
for (int i = 0; i < height; i++) {
for (int j = 0; j < width - 1; j++) {
matrix[i][j] = getEntry(i + 1, j + 1);
}
matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
}
this.tableau = new RealMatrixImpl(matrix);
this.numArtificialVariables = 0;
}
/**
* @param src the source array
* @param dest the destination array
* @param destPos the destination position
*/
private void copyArray(final double[] src, final double[] dest,
final int destPos) {
System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length);
}
/**
* Get the current solution.
* <p>
* {@link #solve} should be called first for this to be the optimal solution.
* </p>
* @return current solution
*/
protected RealPointValuePair getSolution() {
double[] coefficients = new double[getOriginalNumDecisionVariables()];
double mostNegative = getDecisionVariableValue(getOriginalNumDecisionVariables());
for (int i = 0; i < coefficients.length; i++) {
coefficients[i] =
getDecisionVariableValue(i) - (restrictToNonNegative ? 0 : mostNegative);
}
return new RealPointValuePair(coefficients, f.getValue(coefficients));
}
/**
* Get the value of the given decision variable. This is not the actual
* value as it is guaranteed to be >= 0 and thus must be corrected before
* being returned to the user.
*
* @param decisionVariable The index of the decision variable
* @return The value of the given decision variable.
*/
protected double getDecisionVariableValue(final int decisionVariable) {
Integer basicRow = getBasicRow(getNumObjectiveFunctions() + decisionVariable);
return basicRow == null ? 0 : getEntry(basicRow, getRhsOffset());
}
/**
* Subtracts a multiple of one row from another.
* <p>
* After application of this operation, the following will hold:
* minuendRow = minuendRow - multiple * subtrahendRow
* </p>
* @param dividendRow index of the row
* @param divisor value of the divisor
*/
protected void divideRow(final int dividendRow, final double divisor) {
for (int j = 0; j < getWidth(); j++) {
tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor);
}
}
/**
* Subtracts a multiple of one row from another.
* <p>
* After application of this operation, the following will hold:
* minuendRow = minuendRow - multiple * subtrahendRow
* </p>
* @param minuendRow row index
* @param subtrahendRow row index
* @param multiple multiplication factor
*/
protected void subtractRow(final int minuendRow, final int subtrahendRow,
final double multiple) {
for (int j = 0; j < getWidth(); j++) {
tableau.setEntry(minuendRow, j, tableau.getEntry(minuendRow, j) -
multiple * tableau.getEntry(subtrahendRow, j));
}
}
/**
* Get the width of the tableau.
* @return width of the tableau
*/
protected final int getWidth() {
return tableau.getColumnDimension();
}
/**
* Get the height of the tableau.
* @return height of the tableau
*/
protected final int getHeight() {
return tableau.getRowDimension();
}
/** Get an entry of the tableau.
* @param row row index
* @param column column index
* @return entry at (row, column)
*/
protected final double getEntry(final int row, final int column) {
return tableau.getEntry(row, column);
}
/** Set an entry of the tableau.
* @param row row index
* @param column column index
* param value for the entry
*/
protected final void setEntry(final int row, final int column,
final double value) {
tableau.setEntry(row, column, value);
}
/**
* Get the offset of the first slack variable.
* @return offset of the first slack variable
*/
protected final int getSlackVariableOffset() {
return getNumObjectiveFunctions() + numDecisionVariables;
}
/**
* Get the offset of the first artificial variable.
* @return offset of the first artificial variable
*/
protected final int getArtificialVariableOffset() {
return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
}
/**
* Get the offset of the right hand side.
* @return offset of the right hand side
*/
protected final int getRhsOffset() {
return getWidth() - 1;
}
/**
* Get the number of decision variables.
* <p>
* If variables are not restricted to positive values, this will include 1
* extra decision variable to represent the absolute value of the most
* negative variable.
* </p>
* @return number of decision variables
* @see #getOriginalNumDecisionVariables()
*/
protected final int getNumDecisionVariables() {
return numDecisionVariables;
}
/**
* Get the original number of decision variables.
* @return original number of decision variables
* @see #getNumDecisionVariables()
*/
protected final int getOriginalNumDecisionVariables() {
return restrictToNonNegative ? numDecisionVariables : numDecisionVariables - 1;
}
/**
* Get the number of slack variables.
* @return number of slack variables
*/
protected final int getNumSlackVariables() {
return numSlackVariables;
}
/**
* Get the number of artificial variables.
* @return number of artificial variables
*/
protected final int getNumArtificialVariables() {
return numArtificialVariables;
}
/**
* Get the tableau data.
* @return tableau data
*/
protected final double[][] getData() {
return tableau.getData();
}
}

View File

@ -0,0 +1,40 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import org.apache.commons.math.optimization.OptimizationException;
/**
* This class represents exceptions thrown by optimizers when a solution
* escapes to infinity.
* @version $Revision$ $Date$
* @since 2.0
*/
public class UnboundedSolutionException extends OptimizationException {
/** Serializable version identifier. */
private static final long serialVersionUID = 940539497277290619L;
/**
* Simple constructor using a default message.
*/
public UnboundedSolutionException() {
super("unbounded solution");
}
}

View File

@ -0,0 +1,22 @@
<html>
<!--
Licensed to the Apache Software Foundation (ASF) under one or more
contributor license agreements. See the NOTICE file distributed with
this work for additional information regarding copyright ownership.
The ASF licenses this file to You under the Apache License, Version 2.0
(the "License"); you may not use this file except in compliance with
the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-->
<!-- $Revision$ -->
<body>
This package provides optimization algorithms for linear constrained problems.
</body>
</html>

View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add" issue="MATH-246" due-to="Benjamin McCann">
Added an optimizer for constrained linear problems based on 2-phases standard simplex.
</action>
<action dev="luc" type="fix" issue="MATH-177" >
Redesigned the optimization framework for a simpler yet more powerful API.
Added non-linear conjugate gradient optimizer.

View File

@ -0,0 +1,313 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.util.ArrayList;
import java.util.Collection;
import junit.framework.TestCase;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.RealVectorImpl;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealPointValuePair;
public class SimplexSolverTest extends TestCase {
public void testSimplexSolver() throws OptimizationException {
LinearObjectiveFunction f =
new LinearObjectiveFunction(new double[] { 15, 10 }, 7);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 0 }, Relationship.LEQ, 2));
constraints.add(new LinearConstraint(new double[] { 0, 1 }, Relationship.LEQ, 3));
constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 4));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
assertEquals(2.0, solution.getPoint()[0]);
assertEquals(2.0, solution.getPoint()[1]);
assertEquals(57.0, solution.getValue());
}
/**
* With no artificial variables needed (no equals and no greater than
* constraints) we can go straight to Phase 2.
*/
public void testModelWithNoArtificialVars() throws OptimizationException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 15, 10 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 0 }, Relationship.LEQ, 2));
constraints.add(new LinearConstraint(new double[] { 0, 1 }, Relationship.LEQ, 3));
constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.LEQ, 4));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
assertEquals(2.0, solution.getPoint()[0]);
assertEquals(2.0, solution.getPoint()[1]);
assertEquals(50.0, solution.getValue());
}
public void testMinimization() throws OptimizationException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { -2, 1 }, -5);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 2 }, Relationship.LEQ, 6));
constraints.add(new LinearConstraint(new double[] { 3, 2 }, Relationship.LEQ, 12));
constraints.add(new LinearConstraint(new double[] { 0, 1 }, Relationship.GEQ, 0));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
assertEquals(4.0, solution.getPoint()[0]);
assertEquals(0.0, solution.getPoint()[1]);
assertEquals(-13.0, solution.getValue());
}
public void testSolutionWithNegativeDecisionVariable() throws OptimizationException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { -2, 1 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.GEQ, 6));
constraints.add(new LinearConstraint(new double[] { 1, 2 }, Relationship.LEQ, 14));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
assertEquals(-2.0, solution.getPoint()[0]);
assertEquals(8.0, solution.getPoint()[1]);
assertEquals(12.0, solution.getValue());
}
public void testInfeasibleSolution() throws UnboundedSolutionException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 15 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1 }, Relationship.LEQ, 1));
constraints.add(new LinearConstraint(new double[] { 1 }, Relationship.GEQ, 3));
SimplexSolver solver = new SimplexSolver();
try {
solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
fail("An exception should have been thrown.");
} catch (NoFeasibleSolutionException e) {
// expected;
} catch (OptimizationException e) {
fail("wrong exception caught");
}
}
public void testUnboundedSolution() throws NoFeasibleSolutionException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 15, 10 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 0 }, Relationship.EQ, 2));
SimplexSolver solver = new SimplexSolver();
try {
solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
fail("An exception should have been thrown.");
} catch (UnboundedSolutionException e) {
// expected;
} catch (OptimizationException e) {
fail("wrong exception caught");
}
}
public void testRestrictVariablesToNonNegative() throws OptimizationException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 409, 523, 70, 204, 339 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 43, 56, 345, 56, 5 }, Relationship.LEQ, 4567456));
constraints.add(new LinearConstraint(new double[] { 12, 45, 7, 56, 23 }, Relationship.LEQ, 56454));
constraints.add(new LinearConstraint(new double[] { 8, 768, 0, 34, 7456 }, Relationship.LEQ, 1923421));
constraints.add(new LinearConstraint(new double[] { 12342, 2342, 34, 678, 2342 }, Relationship.GEQ, 4356));
constraints.add(new LinearConstraint(new double[] { 45, 678, 76, 52, 23 }, Relationship.EQ, 456356));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
assertEquals(2902.92783505155, solution.getPoint()[0], .0000001);
assertEquals(480.419243986254, solution.getPoint()[1], .0000001);
assertEquals(0.0, solution.getPoint()[2], .0000001);
assertEquals(0.0, solution.getPoint()[3], .0000001);
assertEquals(0.0, solution.getPoint()[4], .0000001);
assertEquals(1438556.7491409, solution.getValue(), .0000001);
}
public void testSomething() throws OptimizationException {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 1, 1 }, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 0));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
assertEquals(0, solution.getValue(), .0000001);
}
public void testLargeModel() throws OptimizationException {
double[] objective = new double[] {
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 12, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
12, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 12, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 12, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 12, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 12, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1};
LinearObjectiveFunction f = new LinearObjectiveFunction(objective, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(equationFromString(objective.length, "x0 + x1 + x2 + x3 - x12 = 0"));
constraints.add(equationFromString(objective.length, "x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 - x13 = 0"));
constraints.add(equationFromString(objective.length, "x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 >= 49"));
constraints.add(equationFromString(objective.length, "x0 + x1 + x2 + x3 >= 42"));
constraints.add(equationFromString(objective.length, "x14 + x15 + x16 + x17 - x26 = 0"));
constraints.add(equationFromString(objective.length, "x18 + x19 + x20 + x21 + x22 + x23 + x24 + x25 - x27 = 0"));
constraints.add(equationFromString(objective.length, "x14 + x15 + x16 + x17 - x12 = 0"));
constraints.add(equationFromString(objective.length, "x18 + x19 + x20 + x21 + x22 + x23 + x24 + x25 - x13 = 0"));
constraints.add(equationFromString(objective.length, "x28 + x29 + x30 + x31 - x40 = 0"));
constraints.add(equationFromString(objective.length, "x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 - x41 = 0"));
constraints.add(equationFromString(objective.length, "x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 >= 49"));
constraints.add(equationFromString(objective.length, "x28 + x29 + x30 + x31 >= 42"));
constraints.add(equationFromString(objective.length, "x42 + x43 + x44 + x45 - x54 = 0"));
constraints.add(equationFromString(objective.length, "x46 + x47 + x48 + x49 + x50 + x51 + x52 + x53 - x55 = 0"));
constraints.add(equationFromString(objective.length, "x42 + x43 + x44 + x45 - x40 = 0"));
constraints.add(equationFromString(objective.length, "x46 + x47 + x48 + x49 + x50 + x51 + x52 + x53 - x41 = 0"));
constraints.add(equationFromString(objective.length, "x56 + x57 + x58 + x59 - x68 = 0"));
constraints.add(equationFromString(objective.length, "x60 + x61 + x62 + x63 + x64 + x65 + x66 + x67 - x69 = 0"));
constraints.add(equationFromString(objective.length, "x60 + x61 + x62 + x63 + x64 + x65 + x66 + x67 >= 51"));
constraints.add(equationFromString(objective.length, "x56 + x57 + x58 + x59 >= 44"));
constraints.add(equationFromString(objective.length, "x70 + x71 + x72 + x73 - x82 = 0"));
constraints.add(equationFromString(objective.length, "x74 + x75 + x76 + x77 + x78 + x79 + x80 + x81 - x83 = 0"));
constraints.add(equationFromString(objective.length, "x70 + x71 + x72 + x73 - x68 = 0"));
constraints.add(equationFromString(objective.length, "x74 + x75 + x76 + x77 + x78 + x79 + x80 + x81 - x69 = 0"));
constraints.add(equationFromString(objective.length, "x84 + x85 + x86 + x87 - x96 = 0"));
constraints.add(equationFromString(objective.length, "x88 + x89 + x90 + x91 + x92 + x93 + x94 + x95 - x97 = 0"));
constraints.add(equationFromString(objective.length, "x88 + x89 + x90 + x91 + x92 + x93 + x94 + x95 >= 51"));
constraints.add(equationFromString(objective.length, "x84 + x85 + x86 + x87 >= 44"));
constraints.add(equationFromString(objective.length, "x98 + x99 + x100 + x101 - x110 = 0"));
constraints.add(equationFromString(objective.length, "x102 + x103 + x104 + x105 + x106 + x107 + x108 + x109 - x111 = 0"));
constraints.add(equationFromString(objective.length, "x98 + x99 + x100 + x101 - x96 = 0"));
constraints.add(equationFromString(objective.length, "x102 + x103 + x104 + x105 + x106 + x107 + x108 + x109 - x97 = 0"));
constraints.add(equationFromString(objective.length, "x112 + x113 + x114 + x115 - x124 = 0"));
constraints.add(equationFromString(objective.length, "x116 + x117 + x118 + x119 + x120 + x121 + x122 + x123 - x125 = 0"));
constraints.add(equationFromString(objective.length, "x116 + x117 + x118 + x119 + x120 + x121 + x122 + x123 >= 49"));
constraints.add(equationFromString(objective.length, "x112 + x113 + x114 + x115 >= 42"));
constraints.add(equationFromString(objective.length, "x126 + x127 + x128 + x129 - x138 = 0"));
constraints.add(equationFromString(objective.length, "x130 + x131 + x132 + x133 + x134 + x135 + x136 + x137 - x139 = 0"));
constraints.add(equationFromString(objective.length, "x126 + x127 + x128 + x129 - x124 = 0"));
constraints.add(equationFromString(objective.length, "x130 + x131 + x132 + x133 + x134 + x135 + x136 + x137 - x125 = 0"));
constraints.add(equationFromString(objective.length, "x140 + x141 + x142 + x143 - x152 = 0"));
constraints.add(equationFromString(objective.length, "x144 + x145 + x146 + x147 + x148 + x149 + x150 + x151 - x153 = 0"));
constraints.add(equationFromString(objective.length, "x144 + x145 + x146 + x147 + x148 + x149 + x150 + x151 >= 59"));
constraints.add(equationFromString(objective.length, "x140 + x141 + x142 + x143 >= 42"));
constraints.add(equationFromString(objective.length, "x154 + x155 + x156 + x157 - x166 = 0"));
constraints.add(equationFromString(objective.length, "x158 + x159 + x160 + x161 + x162 + x163 + x164 + x165 - x167 = 0"));
constraints.add(equationFromString(objective.length, "x154 + x155 + x156 + x157 - x152 = 0"));
constraints.add(equationFromString(objective.length, "x158 + x159 + x160 + x161 + x162 + x163 + x164 + x165 - x153 = 0"));
constraints.add(equationFromString(objective.length, "x83 + x82 - x168 = 0"));
constraints.add(equationFromString(objective.length, "x111 + x110 - x169 = 0"));
constraints.add(equationFromString(objective.length, "x170 - x182 = 0"));
constraints.add(equationFromString(objective.length, "x171 - x183 = 0"));
constraints.add(equationFromString(objective.length, "x172 - x184 = 0"));
constraints.add(equationFromString(objective.length, "x173 - x185 = 0"));
constraints.add(equationFromString(objective.length, "x174 - x186 = 0"));
constraints.add(equationFromString(objective.length, "x175 + x176 - x187 = 0"));
constraints.add(equationFromString(objective.length, "x177 - x188 = 0"));
constraints.add(equationFromString(objective.length, "x178 - x189 = 0"));
constraints.add(equationFromString(objective.length, "x179 - x190 = 0"));
constraints.add(equationFromString(objective.length, "x180 - x191 = 0"));
constraints.add(equationFromString(objective.length, "x181 - x192 = 0"));
constraints.add(equationFromString(objective.length, "x170 - x26 = 0"));
constraints.add(equationFromString(objective.length, "x171 - x27 = 0"));
constraints.add(equationFromString(objective.length, "x172 - x54 = 0"));
constraints.add(equationFromString(objective.length, "x173 - x55 = 0"));
constraints.add(equationFromString(objective.length, "x174 - x168 = 0"));
constraints.add(equationFromString(objective.length, "x177 - x169 = 0"));
constraints.add(equationFromString(objective.length, "x178 - x138 = 0"));
constraints.add(equationFromString(objective.length, "x179 - x139 = 0"));
constraints.add(equationFromString(objective.length, "x180 - x166 = 0"));
constraints.add(equationFromString(objective.length, "x181 - x167 = 0"));
constraints.add(equationFromString(objective.length, "x193 - x205 = 0"));
constraints.add(equationFromString(objective.length, "x194 - x206 = 0"));
constraints.add(equationFromString(objective.length, "x195 - x207 = 0"));
constraints.add(equationFromString(objective.length, "x196 - x208 = 0"));
constraints.add(equationFromString(objective.length, "x197 - x209 = 0"));
constraints.add(equationFromString(objective.length, "x198 + x199 - x210 = 0"));
constraints.add(equationFromString(objective.length, "x200 - x211 = 0"));
constraints.add(equationFromString(objective.length, "x201 - x212 = 0"));
constraints.add(equationFromString(objective.length, "x202 - x213 = 0"));
constraints.add(equationFromString(objective.length, "x203 - x214 = 0"));
constraints.add(equationFromString(objective.length, "x204 - x215 = 0"));
constraints.add(equationFromString(objective.length, "x193 - x182 = 0"));
constraints.add(equationFromString(objective.length, "x194 - x183 = 0"));
constraints.add(equationFromString(objective.length, "x195 - x184 = 0"));
constraints.add(equationFromString(objective.length, "x196 - x185 = 0"));
constraints.add(equationFromString(objective.length, "x197 - x186 = 0"));
constraints.add(equationFromString(objective.length, "x198 + x199 - x187 = 0"));
constraints.add(equationFromString(objective.length, "x200 - x188 = 0"));
constraints.add(equationFromString(objective.length, "x201 - x189 = 0"));
constraints.add(equationFromString(objective.length, "x202 - x190 = 0"));
constraints.add(equationFromString(objective.length, "x203 - x191 = 0"));
constraints.add(equationFromString(objective.length, "x204 - x192 = 0"));
SimplexSolver solver = new SimplexSolver();
RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true);
assertEquals(13366.0, solution.getValue(), .0000001);
//assertEquals(7518.0, solution.getValue(), .0000001);
}
/**
* Converts a test string to a {@link LinearConstraint}.
* Ex: x0 + x1 + x2 + x3 - x12 = 0
*/
private LinearConstraint equationFromString(int numCoefficients, String s) {
Relationship relationship;
if (s.contains(">=")) {
relationship = Relationship.GEQ;
} else if (s.contains("<=")) {
relationship = Relationship.LEQ;
} else if (s.contains("=")) {
relationship = Relationship.EQ;
} else {
throw new IllegalArgumentException();
}
String[] equationParts = s.split("[>|<]?=");
double rhs = Double.parseDouble(equationParts[1].trim());
RealVector lhs = new RealVectorImpl(numCoefficients);
String left = equationParts[0].replaceAll(" ?x", "");
String[] coefficients = left.split(" ");
for (String coefficient : coefficients) {
double value = coefficient.charAt(0) == '-' ? -1 : 1;
int index = Integer.parseInt(coefficient.replaceFirst("[+|-]", "").trim());
lhs.setEntry(index, value);
}
return new LinearConstraint(lhs, relationship, rhs);
}
}

View File

@ -0,0 +1,98 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.linear;
import java.util.ArrayList;
import java.util.Collection;
import org.apache.commons.math.optimization.GoalType;
import junit.framework.TestCase;
public class SimplexTableauTest extends TestCase {
public void testInitialization() {
LinearObjectiveFunction f = createFunction();
Collection<LinearConstraint> constraints = createConstraints();
SimplexTableau tableau =
new SimplexTableau(f, constraints, GoalType.MAXIMIZE, false);
double[][] expectedInitialTableau = {
{-1, 0, -1, -1, 2, 0, 0, 0, -4},
{ 0, 1, -15, -10, 25, 0, 0, 0, 0},
{ 0, 0, 1, 0, -1, 1, 0, 0, 2},
{ 0, 0, 0, 1, -1, 0, 1, 0, 3},
{ 0, 0, 1, 1, -2, 0, 0, 1, 4}
};
assertMatrixEquals(expectedInitialTableau, tableau.getData());
}
public void testdiscardArtificialVariables() {
LinearObjectiveFunction f = createFunction();
Collection<LinearConstraint> constraints = createConstraints();
SimplexTableau tableau =
new SimplexTableau(f, constraints, GoalType.MAXIMIZE, false);
double[][] expectedTableau = {
{ 1, -15, -10, 25, 0, 0, 0},
{ 0, 1, 0, -1, 1, 0, 2},
{ 0, 0, 1, -1, 0, 1, 3},
{ 0, 1, 1, -2, 0, 0, 4}
};
tableau.discardArtificialVariables();
assertMatrixEquals(expectedTableau, tableau.getData());
}
public void testTableauWithNoArtificialVars() {
LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {15, 10}, 0);
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] {1, 0}, Relationship.LEQ, 2));
constraints.add(new LinearConstraint(new double[] {0, 1}, Relationship.LEQ, 3));
constraints.add(new LinearConstraint(new double[] {1, 1}, Relationship.LEQ, 4));
SimplexTableau tableau =
new SimplexTableau(f, constraints, GoalType.MAXIMIZE, false);
double[][] initialTableau = {
{1, -15, -10, 25, 0, 0, 0, 0},
{0, 1, 0, -1, 1, 0, 0, 2},
{0, 0, 1, -1, 0, 1, 0, 3},
{0, 1, 1, -2, 0, 0, 1, 4}
};
assertMatrixEquals(initialTableau, tableau.getData());
}
private LinearObjectiveFunction createFunction() {
return new LinearObjectiveFunction(new double[] {15, 10}, 0);
}
private Collection<LinearConstraint> createConstraints() {
Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
constraints.add(new LinearConstraint(new double[] {1, 0}, Relationship.LEQ, 2));
constraints.add(new LinearConstraint(new double[] {0, 1}, Relationship.LEQ, 3));
constraints.add(new LinearConstraint(new double[] {1, 1}, Relationship.EQ, 4));
return constraints;
}
private void assertMatrixEquals(double[][] expected, double[][] result) {
assertEquals("Wrong number of rows.", expected.length, result.length);
for (int i = 0; i < expected.length; i++) {
assertEquals("Wrong number of columns.", expected[i].length, result[i].length);
for (int j = 0; j < expected[i].length; j++) {
assertEquals("Wrong value at position [" + i + "," + j + "]", expected[i][j], result[i][j]);
}
}
}
}